gb_lisa.w

模拟器提供了一个简单易用的平台

% This file is part of the Stanford GraphBase (c) Stanford University 1993
@i boilerplate.w %<< legal stuff: PLEASE READ IT BEFORE MAKING ANY CHANGES!
@i gb_types.w

\def\title{GB\_\,LISA}

\prerequisites{GB\_\,GRAPH}{GB\_\,IO}
@* Introduction. This GraphBase module contains the |lisa| subroutine,
which creates rectangular matrices of data based on Leonardo da Vinci's
@^Vinci, Leonardo da@>
{\sl Gioconda\/} (aka Mona Lisa). It also contains the |plane_lisa|
subroutine, which constructs undirected planar graphs based on |lisa|,
and the |bi_lisa| subroutine, which constructs undirected bipartite graphs.
Another example of the use of |lisa| can be
found in the demo program {\sc ASSIGN\_LISA}.

@d plane_lisa p_lisa /* abbreviation for Procrustean external linkage */

@(gb_lisa.h@>=
#define plane_lisa p_lisa
extern long* lisa();
extern Graph *plane_lisa();
extern Graph *bi_lisa();

@ The subroutine call |lisa(m,n,d,m0,m1,n0,n1,d0,d1,area)|
constructs an $m\times n$ matrix of integers in the range
$[0\,.\,.\,d\mskip1mu]$,
based on the information in \.{lisa.dat}. Storage space for the matrix is
allocated in the memory area called |area|, using the normal GraphBase
conventions explained in {\sc GB\_\,GRAPH}.
The entries of the matrix can be regarded as pixel data, with
0~representing black and $d$~representing white, and with intermediate
values representing shades of gray.

The data in \.{lisa.dat} has 360 rows and 250 columns. The rows are numbered
0 to 359 from top to bottom, and the columns are numbered 0 to 249 from left
to right. The output of |lisa| is generated from a rectangular section
of the picture consisting of |m1-m0| rows and |n1-n0| columns; more
precisely, |lisa| uses the data in positions $(k,l)$ for
|m0<=k<m1| and |n0<=l<n1|.

One way to understand the process of mapping |M=m1-m0| rows and |N=n1-n0|
columns of input into |m|~rows and |n|~columns of output is to imagine
a giant matrix of $mM$ rows and $nN$ columns in which the original input
data has been replicated as an $M\times N$ array of submatrices of
size $m\times n$; each of the submatrices contains $mn$ identical pixel
values. We can also regard the giant matrix as an $m\times n$ array of
submatrices of size $M\times N$. The pixel values to be output are obtained
by averaging the $M_{}N$ pixel values in the submatrices of this second
interpretation.

More precisely, the output pixel value in a given row and column is obtained
in two steps. First we sum the $M_{}N$ entries in the corresponding submatrix
of the giant matrix, obtaining a value $D$ between 0 and~$255M_{}N$. Then we
scale the value~$D$ linearly into the desired final range
$[0\,.\,.\,d\mskip1mu]$ by
setting the result to~0 if |D<d0|, to~$d$ if |D>=d1|, and to
$\lfloor d(D-|d0|)/(|d1|-|d0|)\rfloor$ if |d0<=D<d1|.

@d MAX_M 360 /* the total number of rows of input data */
@d MAX_N 250 /* the total number of columns of input data */
@d MAX_D 255 /* maximum pixel value in the input data */

@ Default parameter values are automatically substituted when |m|, |n|, |d|,
|m1|, |n1|, and/or |d1| are given as~0: If |m1=0| or |m1>360|,
|m1|~is changed to 360; if |n1=0| or |n1>250|, |n1|~is
changed to~250. Then if |m| is zero, it is changed
to~|m1-m0|; if |n| is zero, it is changed to~|n1-n0|.
If |d| is zero, it is changed to~255.
If |d1| is zero, it is changed to |255(m1-m0)(n1-n0)|.
After these substitutions have been made, the parameters must satisfy
$$\hbox{|m0<m1|, \qquad|n0<n1|, \qquad and \qquad |d0<d1|.}$$

Examples: The call |@t\\{lisa\_pix}@>=lisa(0,0,0,0,0,0,0,0,0,area)|
is equivalent to the call
|@t\\{lisa\_pix}@>=lisa(360,250,255,0,360,0,250,0,255*360*250,area)|;
this special case delivers the original \.{lisa.dat} data as a
$360\times250$ array of integers in the range $[0\,.\,.\,255]$. You
can access the pixel in row~$k$ and column~$l$ by writing
$$\hbox{|*(@[@t\\{lisa\_pix}@>@]+n*k+l)|}\,,$$ where |n| in this case is
250. A square array extracted from the top part of the picture,
leaving out Mona's hands at the bottom, can be obtained by calling
|lisa(250,250,255,0,250,0,250,0,0,area)|.

The call |lisa(36,25,25500,0,0,0,0,0,0,area)| gives a $36\times25$ array
of pixel values in the range $[0\,.\,.\,25500]$, obtained by summing
$10\times10$ subsquares of the original data.

The call |lisa(100,100,100,0,0,0,0,0,0,area)| gives a $100\times100$ array
of pixel values in the range $[0\,.\,.\,100]$; in this case the original
data is effectively broken into subpixels and averaged appropriately.
Notice that each output pixel in this example comes from 3.6 input
rows and 2.5 input columns; therefore the image is being distorted
(compressed vertically). However, our GraphBase applications are generally
interested more in combinatorial test data, not in images per~se.
If |(m1-m0)/m=(n1-n0)/n|, the output of |lisa| will represent ``square
pixels.'' But if |(m1-m0)/m<(n1-n0)/n|, a halftone generated from the
output will be compressed in the horizontal dimension; if
|(m1-m0)/m>(n1-n0)/n|, it will be compressed in the vertical dimension.

If you want to reduce the original image to binary data, with the value~0
wherever the original pixels are less than some threshold value~|t|
and the value~1 whenever they are |t| or more, call
|lisa(m,n,1,m0,m1,n0,n1,@t}\penalty0{@>0,t*(m1-m0)*(n1-n0),area)|.

The subroutine call |lisa(1000,1000,255,0,250,0,250,0,0,area)| produces a
million pixels from the upper part of the original image. This matrix
contains more entries than the original data in \.{lisa.dat}, but of course
it is not any more accurate; it has simply been obtained by linear
interpolation---in fact, by replicating the original
data in $4\times4$ subarrays.

Mona Lisa's famous smile appears in the $16\times32$ subarray defined by
|m0=94|, |m1=110|, |n0=97|, |n1=129|. The |smile| macro makes this
easily accessible. (See also |eyes|.)

A string |lisa_id| is constructed, showing the actual parameter values
used by |lisa| after defaults have been supplied.
The |area| parameter is omitted from this string.

@<gb_lisa.h@>=
#define smile @t\quad@> m0=94,m1=110,n0=97,n1=129 /* $16\times32$ */
#define eyes @t\quad@> m0=61,m1=80,n0=91,n1=140 /* $20\times50$ */
extern char lisa_id[];

@ @<Global variables@>=
char lisa_id[]=
  "lisa(360,250,9999999999,359,360,249,250,9999999999,9999999999)";

@ If the |lisa| routine encounters a problem, it returns |NULL|
(\.{NULL}), after putting a nonzero number into the external variable
|panic_code|. This code number identifies the type of failure.
Otherwise |lisa| returns a pointer to the newly created array. (The
external variable |panic_code| is defined in {\sc GB\_\,GRAPH}.)

@d panic(c) @+{@+panic_code=c;@+gb_trouble_code=0;@+return NULL;@+}

@ The \CEE/ file \.{gb\_lisa.c} begins as follows. (Other subroutines
come later.)

@p
#include "gb_io.h" /* we will use the {\sc GB\_\,IO} routines for input */
#include "gb_graph.h" /* we will use the {\sc GB\_\,GRAPH} data structures */
@h@#
@<Global variables@>@;
@<Private variables@>@;
@<Private subroutines@>@;
@#
long *lisa(m,n,d,m0,m1,n0,n1,d0,d1,area)
  unsigned long m,n; /* number of rows and columns desired */
  unsigned long d; /* maximum pixel value desired */
  unsigned long m0,m1; /* input will be from rows $[|m0|\,.\,.\,|m1|)$ */
  unsigned long n0,n1; /* and from columns $[|n0|\,.\,.\,|n1|)$ */
  unsigned long d0,d1; /* lower and upper threshold of raw pixel scores */
  Area area; /* where to allocate the matrix that will be output */
{@+@<Local variables for |lisa|@>@;@#
  @<Check the parameters and adjust them for defaults@>;
  @<Allocate the matrix@>;
  @<Read \.{lisa.dat} and map it to the desired output form@>;
  return matx;
}

@ @<Local variables for |lisa|@>=
long *matx=NULL; /* the matrix constructed by |lisa| */
register long k,l; /* the current row and column of output */
register long i,j; /* all-purpose indices */
long cap_M,cap_N; /* |m1-m0| and |n1-n0|, dimensions of the input */
long cap_D; /* |d1-d0|, scale factor */

@ @<Check the param...@>=
if (m1==0 || m1>MAX_M) m1=MAX_M;
if (m1<=m0) panic(bad_specs+1); /* |m0| must be less than |m1| */
if (n1==0 || n1>MAX_N) n1=MAX_N;
if (n1<=n0) panic(bad_specs+2); /* |n0| must be less than |n1| */
cap_M=m1-m0;@+cap_N=n1-n0;
if (m==0) m=cap_M;
if (n==0) n=cap_N;
if (d==0) d=MAX_D;
if (d1==0) d1=MAX_D*cap_M*cap_N;
if (d1<=d0) panic(bad_specs+3); /* |d0| must be less than |d1| */
if (d1>=0x80000000) panic(bad_specs+4); /* |d1| must be less than $2^{31}$ */
cap_D=d1-d0;
sprintf(lisa_id,"lisa(%lu,%lu,%lu,%lu,%lu,%lu,%lu,%lu,%lu)",
   m,n,d,m0,m1,n0,n1,d0,d1);

@ @<Allocate the matrix@>=
matx=gb_typed_alloc(m*n,long,area);
if (gb_trouble_code) panic(no_room+1); /* no room for the output data */

@ @<Read \.{lisa.dat} and map it to the desired output form@>=
@<Open the data file, skipping unwanted rows at the beginning@>;
@<Generate the $m$ rows of output@>;
@<Close the data file, skipping unwanted rows at the end@>;

@* Elementary image processing.
As mentioned in the introduction, we can visualize the input as a giant
$mM\times nN$ matrix, into which an $M\times N$ image is placed by replication
of pixel values, and from which an $m\times n$ image is derived by summation
of pixel values and subsequent scaling. Here |M=m1-m0| and |N=n1-n0|.

Let $(\kappa,\lambda)$ be a position in the giant matrix, where $0\le\kappa<mM$
and $0\le\lambda<nN$. The corresponding indices of the input image are
then $\bigl(|m0|+\lfloor\kappa/m\rfloor, |n0|+\lfloor\lambda/n\rfloor\bigr)$,
and the corresponding indices of the output image are
$\bigl(\lfloor\kappa/M\rfloor,\lfloor\lambda/N\rfloor\bigr)$. Our main job
is to compute the sum of all pixel values that lie in each given row~|k|
and column~|l| of the output image. Many elements are repeated in
the sum, so we want to use multiplication instead of simple addition whenever
possible.

For example, let's consider the inner loop first, the loop on $l$ and
$\lambda$.  Suppose $n=3$, and suppose the input pixels in the current
row of interest are $\langle a_0,\ldots,a_{N-1}\rangle$. Then if $N=3$,
we want to compute the output pixels $\langle3a_0,3a_1,3a_2\rangle$;
if $N=4$, we want to compute
$\langle3a_0+a_1,2a_1+2a_2,a_2+3a_3\rangle$; if $N=2$, we want to
compute $\langle2a_0,a_0+a_1,2a_1\rangle$. The logic for doing this
computation with the proper timing can be expressed conveniently in
terms of four local variables:

@<Local variables for |lisa|@>=
long *cur_pix; /* current position within |in_row| */
long lambda; /* right boundary in giant for the input pixel in |cur_pix| */
long lam; /* the first giant column not yet used in the current row */
long next_lam; /* right boundary in giant for the output pixel in column~|l| */

@ @<Process one row of pixel sums, multiplying them by~|f|@>=
lambda=n;@+cur_pix=in_row+n0;
for (l=lam=0; l<n; l++) {@+register long sum=0;
  next_lam=lam+cap_N;
  do@+{@+register long nl; /* giant column where something new might happen */
    if (lam>=lambda) cur_pix++,lambda+=n;
    if (lambda<next_lam) nl=lambda;
    else nl=next_lam;
    sum+=(nl-lam)*(*cur_pix);
    lam=nl;
  }@+while (lam<next_lam);
  *(out_row+l)+=f*sum;
}

@ The outer loop (on $k$ and $\kappa$) is similar but slightly more
complicated, because it deals with a vector of sums instead of a single
sum and because it must invoke the input routine when we're done
with a row of input data.

%Generate them rows...
@<Generate the $m$ rows of output@>=
kappa=0;
out_row=matx;
for (k=kap=0; k<m;k++) {
  for (l=0;l<n;l++) *(out_row+l)=0; /* clear the vector of sums */
  next_kap=kap+cap_M;
  do@+{@+register long nk; /* giant row where something new might happen */
    if (kap>=kappa) {
      @<Read a row of input into |in_row|@>;
      kappa+=m;
    }
    if (kappa<next_kap) nk=kappa;
    else nk=next_kap;
    f=nk-kap;
    @<Process one...@>;
    kap=nk;
  }@+while (kap<next_kap);
  for (l=0; l<n; l++,out_row++) /* note that |out_row| will advance by~|n| */
    @<Scale the sum found in |*out_row|@>;
}

@ @<Local variables for |lisa|@>=
long kappa; /* bottom boundary in giant for the input pixels in |in_row| */
long kap; /* the first giant row not yet used */
long next_kap; /* bottom boundary in giant for the output pixel in row~|k| */
long f; /* factor by which current input sums should be replicated */
long *out_row; /* current position in |matx| */

@* Integer scaling.
Here's a general-purpose routine to compute $\lfloor na/b\rfloor$ exactly
without risking integer overflow, given integers $n\ge0$ and $0<a\le b$.
The idea is to solve the problem first for $n/2$, if $n$ is too large.

We are careful to precompute values so that integer overflow cannot
occur when $b$ is very large.

@d el_gordo 0x7fffffff /* $2^{31}-1$, the largest single-precision |long| */

@<Private sub...@>=
static long na_over_b(n,a,b)
  long n,a,b;
{@+long nmax=el_gordo/a; /* the largest $n$ such that $na$ doesn't overflow */
  register long r,k,q,br;
  long a_thresh, b_thresh;
  if (n<=nmax) return (n*a)/b;
  a_thresh=b-a;
  b_thresh=(b+1)>>1; /* $\lceil b/2\rceil$ */
  k=0;
  do@+{@+bit[k]=n&1; /* save the least significant bit of $n$ */
    n>>=1; /* and shift it out */
    k++;
  }@+while (n>nmax);
  r=n*a;@+ q=r/b;@+ r=r-q*b;
  @<Maintain quotient |q| and remainder |r| while increasing |n|
    back to its original value $2^kn+(|bit|[k-1]\ldots |bit|[0])_2$@>;
  return q;
}

@ @<Private var...@>=
static long bit[30]; /* bits shifted out of |n| */

@ @<Maintain quotient...@>=
do@+{@+k--;@+ q<<=1;
  if (r<b_thresh) r<<=1;
  else q++,br=(b-r)<<1,r=b-br;
  if (bit[k]) {
    if (r<a_thresh) r+=a;
    else q++,r-=a_thresh;
  }
}@+while (k);

@ @<Scale the sum found in |*out_row|@>=
if (*out_row<=d0) *out_row=0;
else if (*out_row>=d1) *out_row=d;
else *out_row=na_over_b(d,*out_row-d0,cap_D);

@* Input data format.
The file \.{lisa.dat} contains 360 rows of pixel data, and each row
appears on five consecutive lines of the file. The first four lines contain
the data for 60 pixels; each sequence of four pixels is represented by five
radix-85 digits, using the |icode| mapping of {\sc GB\_\,IO}.