gb_lisa.w
上传用户:rrhhcc
上传日期:2015-12-11
资源大小:54129k
文件大小:27k
源码类别:

通讯编程

开发平台:

Visual C++

  1. % This file is part of the Stanford GraphBase (c) Stanford University 1993
  2. @i boilerplate.w %<< legal stuff: PLEASE READ IT BEFORE MAKING ANY CHANGES!
  3. @i gb_types.w
  4. deftitle{GB_,LISA}
  5. prerequisites{GB_,GRAPH}{GB_,IO}
  6. @* Introduction. This GraphBase module contains the |lisa| subroutine,
  7. which creates rectangular matrices of data based on Leonardo da Vinci's
  8. @^Vinci, Leonardo da@>
  9. {sl Gioconda/} (aka Mona Lisa). It also contains the |plane_lisa|
  10. subroutine, which constructs undirected planar graphs based on |lisa|,
  11. and the |bi_lisa| subroutine, which constructs undirected bipartite graphs.
  12. Another example of the use of |lisa| can be
  13. found in the demo program {sc ASSIGN_LISA}.
  14. @d plane_lisa p_lisa /* abbreviation for Procrustean external linkage */
  15. @(gb_lisa.h@>=
  16. #define plane_lisa p_lisa
  17. extern long* lisa();
  18. extern Graph *plane_lisa();
  19. extern Graph *bi_lisa();
  20. @ The subroutine call |lisa(m,n,d,m0,m1,n0,n1,d0,d1,area)|
  21. constructs an $mtimes n$ matrix of integers in the range
  22. $[0,.,.,dmskip1mu]$,
  23. based on the information in .{lisa.dat}. Storage space for the matrix is
  24. allocated in the memory area called |area|, using the normal GraphBase
  25. conventions explained in {sc GB_,GRAPH}.
  26. The entries of the matrix can be regarded as pixel data, with
  27. 0~representing black and $d$~representing white, and with intermediate
  28. values representing shades of gray.
  29. The data in .{lisa.dat} has 360 rows and 250 columns. The rows are numbered
  30. 0 to 359 from top to bottom, and the columns are numbered 0 to 249 from left
  31. to right. The output of |lisa| is generated from a rectangular section
  32. of the picture consisting of |m1-m0| rows and |n1-n0| columns; more
  33. precisely, |lisa| uses the data in positions $(k,l)$ for
  34. |m0<=k<m1| and |n0<=l<n1|.
  35. One way to understand the process of mapping |M=m1-m0| rows and |N=n1-n0|
  36. columns of input into |m|~rows and |n|~columns of output is to imagine
  37. a giant matrix of $mM$ rows and $nN$ columns in which the original input
  38. data has been replicated as an $Mtimes N$ array of submatrices of
  39. size $mtimes n$; each of the submatrices contains $mn$ identical pixel
  40. values. We can also regard the giant matrix as an $mtimes n$ array of
  41. submatrices of size $Mtimes N$. The pixel values to be output are obtained
  42. by averaging the $M_{}N$ pixel values in the submatrices of this second
  43. interpretation.
  44. More precisely, the output pixel value in a given row and column is obtained
  45. in two steps. First we sum the $M_{}N$ entries in the corresponding submatrix
  46. of the giant matrix, obtaining a value $D$ between 0 and~$255M_{}N$. Then we
  47. scale the value~$D$ linearly into the desired final range
  48. $[0,.,.,dmskip1mu]$ by
  49. setting the result to~0 if |D<d0|, to~$d$ if |D>=d1|, and to
  50. $lfloor d(D-|d0|)/(|d1|-|d0|)rfloor$ if |d0<=D<d1|.
  51. @d MAX_M 360 /* the total number of rows of input data */
  52. @d MAX_N 250 /* the total number of columns of input data */
  53. @d MAX_D 255 /* maximum pixel value in the input data */
  54. @ Default parameter values are automatically substituted when |m|, |n|, |d|,
  55. |m1|, |n1|, and/or |d1| are given as~0: If |m1=0| or |m1>360|,
  56. |m1|~is changed to 360; if |n1=0| or |n1>250|, |n1|~is
  57. changed to~250. Then if |m| is zero, it is changed
  58. to~|m1-m0|; if |n| is zero, it is changed to~|n1-n0|.
  59. If |d| is zero, it is changed to~255.
  60. If |d1| is zero, it is changed to |255(m1-m0)(n1-n0)|.
  61. After these substitutions have been made, the parameters must satisfy
  62. $$hbox{|m0<m1|, qquad|n0<n1|, qquad and qquad |d0<d1|.}$$
  63. Examples: The call |@t\{lisa_pix}@>=lisa(0,0,0,0,0,0,0,0,0,area)|
  64. is equivalent to the call
  65. |@t\{lisa_pix}@>=lisa(360,250,255,0,360,0,250,0,255*360*250,area)|;
  66. this special case delivers the original .{lisa.dat} data as a
  67. $360times250$ array of integers in the range $[0,.,.,255]$. You
  68. can access the pixel in row~$k$ and column~$l$ by writing
  69. $$hbox{|*(@[@t\{lisa_pix}@>@]+n*k+l)|},,$$ where |n| in this case is
  70. 250. A square array extracted from the top part of the picture,
  71. leaving out Mona's hands at the bottom, can be obtained by calling
  72. |lisa(250,250,255,0,250,0,250,0,0,area)|.
  73. The call |lisa(36,25,25500,0,0,0,0,0,0,area)| gives a $36times25$ array
  74. of pixel values in the range $[0,.,.,25500]$, obtained by summing
  75. $10times10$ subsquares of the original data.
  76. The call |lisa(100,100,100,0,0,0,0,0,0,area)| gives a $100times100$ array
  77. of pixel values in the range $[0,.,.,100]$; in this case the original
  78. data is effectively broken into subpixels and averaged appropriately.
  79. Notice that each output pixel in this example comes from 3.6 input
  80. rows and 2.5 input columns; therefore the image is being distorted
  81. (compressed vertically). However, our GraphBase applications are generally
  82. interested more in combinatorial test data, not in images per~se.
  83. If |(m1-m0)/m=(n1-n0)/n|, the output of |lisa| will represent ``square
  84. pixels.'' But if |(m1-m0)/m<(n1-n0)/n|, a halftone generated from the
  85. output will be compressed in the horizontal dimension; if
  86. |(m1-m0)/m>(n1-n0)/n|, it will be compressed in the vertical dimension.
  87. If you want to reduce the original image to binary data, with the value~0
  88. wherever the original pixels are less than some threshold value~|t|
  89. and the value~1 whenever they are |t| or more, call
  90. |lisa(m,n,1,m0,m1,n0,n1,@t}penalty0{@>0,t*(m1-m0)*(n1-n0),area)|.
  91. The subroutine call |lisa(1000,1000,255,0,250,0,250,0,0,area)| produces a
  92. million pixels from the upper part of the original image. This matrix
  93. contains more entries than the original data in .{lisa.dat}, but of course
  94. it is not any more accurate; it has simply been obtained by linear
  95. interpolation---in fact, by replicating the original
  96. data in $4times4$ subarrays.
  97. Mona Lisa's famous smile appears in the $16times32$ subarray defined by
  98. |m0=94|, |m1=110|, |n0=97|, |n1=129|. The |smile| macro makes this
  99. easily accessible. (See also |eyes|.)
  100. A string |lisa_id| is constructed, showing the actual parameter values
  101. used by |lisa| after defaults have been supplied.
  102. The |area| parameter is omitted from this string.
  103. @<gb_lisa.h@>=
  104. #define smile @tquad@> m0=94,m1=110,n0=97,n1=129 /* $16times32$ */
  105. #define eyes @tquad@> m0=61,m1=80,n0=91,n1=140 /* $20times50$ */
  106. extern char lisa_id[];
  107. @ @<Global variables@>=
  108. char lisa_id[]=
  109.   "lisa(360,250,9999999999,359,360,249,250,9999999999,9999999999)";
  110. @ If the |lisa| routine encounters a problem, it returns |NULL|
  111. (.{NULL}), after putting a nonzero number into the external variable
  112. |panic_code|. This code number identifies the type of failure.
  113. Otherwise |lisa| returns a pointer to the newly created array. (The
  114. external variable |panic_code| is defined in {sc GB_,GRAPH}.)
  115. @d panic(c) @+{@+panic_code=c;@+gb_trouble_code=0;@+return NULL;@+}
  116. @ The CEE/ file .{gb_lisa.c} begins as follows. (Other subroutines
  117. come later.)
  118. @p
  119. #include "gb_io.h" /* we will use the {sc GB_,IO} routines for input */
  120. #include "gb_graph.h" /* we will use the {sc GB_,GRAPH} data structures */
  121. @h@#
  122. @<Global variables@>@;
  123. @<Private variables@>@;
  124. @<Private subroutines@>@;
  125. @#
  126. long *lisa(m,n,d,m0,m1,n0,n1,d0,d1,area)
  127.   unsigned long m,n; /* number of rows and columns desired */
  128.   unsigned long d; /* maximum pixel value desired */
  129.   unsigned long m0,m1; /* input will be from rows $[|m0|,.,.,|m1|)$ */
  130.   unsigned long n0,n1; /* and from columns $[|n0|,.,.,|n1|)$ */
  131.   unsigned long d0,d1; /* lower and upper threshold of raw pixel scores */
  132.   Area area; /* where to allocate the matrix that will be output */
  133. {@+@<Local variables for |lisa|@>@;@#
  134.   @<Check the parameters and adjust them for defaults@>;
  135.   @<Allocate the matrix@>;
  136.   @<Read .{lisa.dat} and map it to the desired output form@>;
  137.   return matx;
  138. }
  139. @ @<Local variables for |lisa|@>=
  140. long *matx=NULL; /* the matrix constructed by |lisa| */
  141. register long k,l; /* the current row and column of output */
  142. register long i,j; /* all-purpose indices */
  143. long cap_M,cap_N; /* |m1-m0| and |n1-n0|, dimensions of the input */
  144. long cap_D; /* |d1-d0|, scale factor */
  145. @ @<Check the param...@>=
  146. if (m1==0 || m1>MAX_M) m1=MAX_M;
  147. if (m1<=m0) panic(bad_specs+1); /* |m0| must be less than |m1| */
  148. if (n1==0 || n1>MAX_N) n1=MAX_N;
  149. if (n1<=n0) panic(bad_specs+2); /* |n0| must be less than |n1| */
  150. cap_M=m1-m0;@+cap_N=n1-n0;
  151. if (m==0) m=cap_M;
  152. if (n==0) n=cap_N;
  153. if (d==0) d=MAX_D;
  154. if (d1==0) d1=MAX_D*cap_M*cap_N;
  155. if (d1<=d0) panic(bad_specs+3); /* |d0| must be less than |d1| */
  156. if (d1>=0x80000000) panic(bad_specs+4); /* |d1| must be less than $2^{31}$ */
  157. cap_D=d1-d0;
  158. sprintf(lisa_id,"lisa(%lu,%lu,%lu,%lu,%lu,%lu,%lu,%lu,%lu)",
  159.    m,n,d,m0,m1,n0,n1,d0,d1);
  160. @ @<Allocate the matrix@>=
  161. matx=gb_typed_alloc(m*n,long,area);
  162. if (gb_trouble_code) panic(no_room+1); /* no room for the output data */
  163. @ @<Read .{lisa.dat} and map it to the desired output form@>=
  164. @<Open the data file, skipping unwanted rows at the beginning@>;
  165. @<Generate the $m$ rows of output@>;
  166. @<Close the data file, skipping unwanted rows at the end@>;
  167. @* Elementary image processing.
  168. As mentioned in the introduction, we can visualize the input as a giant
  169. $mMtimes nN$ matrix, into which an $Mtimes N$ image is placed by replication
  170. of pixel values, and from which an $mtimes n$ image is derived by summation
  171. of pixel values and subsequent scaling. Here |M=m1-m0| and |N=n1-n0|.
  172. Let $(kappa,lambda)$ be a position in the giant matrix, where $0lekappa<mM$
  173. and $0lelambda<nN$. The corresponding indices of the input image are
  174. then $bigl(|m0|+lfloorkappa/mrfloor, |n0|+lfloorlambda/nrfloorbigr)$,
  175. and the corresponding indices of the output image are
  176. $bigl(lfloorkappa/Mrfloor,lfloorlambda/Nrfloorbigr)$. Our main job
  177. is to compute the sum of all pixel values that lie in each given row~|k|
  178. and column~|l| of the output image. Many elements are repeated in
  179. the sum, so we want to use multiplication instead of simple addition whenever
  180. possible.
  181. For example, let's consider the inner loop first, the loop on $l$ and
  182. $lambda$.  Suppose $n=3$, and suppose the input pixels in the current
  183. row of interest are $langle a_0,ldots,a_{N-1}rangle$. Then if $N=3$,
  184. we want to compute the output pixels $langle3a_0,3a_1,3a_2rangle$;
  185. if $N=4$, we want to compute
  186. $langle3a_0+a_1,2a_1+2a_2,a_2+3a_3rangle$; if $N=2$, we want to
  187. compute $langle2a_0,a_0+a_1,2a_1rangle$. The logic for doing this
  188. computation with the proper timing can be expressed conveniently in
  189. terms of four local variables:
  190. @<Local variables for |lisa|@>=
  191. long *cur_pix; /* current position within |in_row| */
  192. long lambda; /* right boundary in giant for the input pixel in |cur_pix| */
  193. long lam; /* the first giant column not yet used in the current row */
  194. long next_lam; /* right boundary in giant for the output pixel in column~|l| */
  195. @ @<Process one row of pixel sums, multiplying them by~|f|@>=
  196. lambda=n;@+cur_pix=in_row+n0;
  197. for (l=lam=0; l<n; l++) {@+register long sum=0;
  198.   next_lam=lam+cap_N;
  199.   do@+{@+register long nl; /* giant column where something new might happen */
  200.     if (lam>=lambda) cur_pix++,lambda+=n;
  201.     if (lambda<next_lam) nl=lambda;
  202.     else nl=next_lam;
  203.     sum+=(nl-lam)*(*cur_pix);
  204.     lam=nl;
  205.   }@+while (lam<next_lam);
  206.   *(out_row+l)+=f*sum;
  207. }
  208. @ The outer loop (on $k$ and $kappa$) is similar but slightly more
  209. complicated, because it deals with a vector of sums instead of a single
  210. sum and because it must invoke the input routine when we're done
  211. with a row of input data.
  212. %Generate them rows...
  213. @<Generate the $m$ rows of output@>=
  214. kappa=0;
  215. out_row=matx;
  216. for (k=kap=0; k<m;k++) {
  217.   for (l=0;l<n;l++) *(out_row+l)=0; /* clear the vector of sums */
  218.   next_kap=kap+cap_M;
  219.   do@+{@+register long nk; /* giant row where something new might happen */
  220.     if (kap>=kappa) {
  221.       @<Read a row of input into |in_row|@>;
  222.       kappa+=m;
  223.     }
  224.     if (kappa<next_kap) nk=kappa;
  225.     else nk=next_kap;
  226.     f=nk-kap;
  227.     @<Process one...@>;
  228.     kap=nk;
  229.   }@+while (kap<next_kap);
  230.   for (l=0; l<n; l++,out_row++) /* note that |out_row| will advance by~|n| */
  231.     @<Scale the sum found in |*out_row|@>;
  232. }
  233. @ @<Local variables for |lisa|@>=
  234. long kappa; /* bottom boundary in giant for the input pixels in |in_row| */
  235. long kap; /* the first giant row not yet used */
  236. long next_kap; /* bottom boundary in giant for the output pixel in row~|k| */
  237. long f; /* factor by which current input sums should be replicated */
  238. long *out_row; /* current position in |matx| */
  239. @* Integer scaling.
  240. Here's a general-purpose routine to compute $lfloor na/brfloor$ exactly
  241. without risking integer overflow, given integers $nge0$ and $0<ale b$.
  242. The idea is to solve the problem first for $n/2$, if $n$ is too large.
  243. We are careful to precompute values so that integer overflow cannot
  244. occur when $b$ is very large.
  245. @d el_gordo 0x7fffffff /* $2^{31}-1$, the largest single-precision |long| */
  246. @<Private sub...@>=
  247. static long na_over_b(n,a,b)
  248.   long n,a,b;
  249. {@+long nmax=el_gordo/a; /* the largest $n$ such that $na$ doesn't overflow */
  250.   register long r,k,q,br;
  251.   long a_thresh, b_thresh;
  252.   if (n<=nmax) return (n*a)/b;
  253.   a_thresh=b-a;
  254.   b_thresh=(b+1)>>1; /* $lceil b/2rceil$ */
  255.   k=0;
  256.   do@+{@+bit[k]=n&1; /* save the least significant bit of $n$ */
  257.     n>>=1; /* and shift it out */
  258.     k++;
  259.   }@+while (n>nmax);
  260.   r=n*a;@+ q=r/b;@+ r=r-q*b;
  261.   @<Maintain quotient |q| and remainder |r| while increasing |n|
  262.     back to its original value $2^kn+(|bit|[k-1]ldots |bit|[0])_2$@>;
  263.   return q;
  264. }
  265. @ @<Private var...@>=
  266. static long bit[30]; /* bits shifted out of |n| */
  267. @ @<Maintain quotient...@>=
  268. do@+{@+k--;@+ q<<=1;
  269.   if (r<b_thresh) r<<=1;
  270.   else q++,br=(b-r)<<1,r=b-br;
  271.   if (bit[k]) {
  272.     if (r<a_thresh) r+=a;
  273.     else q++,r-=a_thresh;
  274.   }
  275. }@+while (k);
  276. @ @<Scale the sum found in |*out_row|@>=
  277. if (*out_row<=d0) *out_row=0;
  278. else if (*out_row>=d1) *out_row=d;
  279. else *out_row=na_over_b(d,*out_row-d0,cap_D);
  280. @* Input data format.
  281. The file .{lisa.dat} contains 360 rows of pixel data, and each row
  282. appears on five consecutive lines of the file. The first four lines contain
  283. the data for 60 pixels; each sequence of four pixels is represented by five
  284. radix-85 digits, using the |icode| mapping of {sc GB_,IO}.
  285. The fifth and final line of each row contains $4+4+2=10$ more pixels,
  286. represented as $5+5+3$ radix-85 digits.
  287. @<Open the data file, skipping unwanted rows at the beginning@>=
  288. if (gb_open("lisa.dat")!=0)
  289.   panic(early_data_fault); /* couldn't open the file; |io_errors| tells why */
  290. for (i=0;i<m0;i++)
  291.   for (j=0;j<5;j++) gb_newline(); /* ignore one row of data */
  292. @ @<Close the data file, skipping unwanted rows at the end@>=
  293. for (i=m1;i<MAX_M;i++)
  294.   for (j=0;j<5;j++) gb_newline(); /* ignore one row of data */
  295. if (gb_close()!=0)
  296.   panic(late_data_fault);
  297.    /* checksum or other failure in data file; see |io_errors| */
  298. @ @<Read a row of input into |in_row|@>=
  299. {@+register long dd;
  300.   for (j=15,cur_pix=&in_row[0];;cur_pix+=4) {
  301.     dd=gb_digit(85);@+dd=dd*85+gb_digit(85);@+dd=dd*85+gb_digit(85);
  302.     if (cur_pix==&in_row[MAX_N-2]) break;
  303.     dd=dd*85+gb_digit(85);@+dd=dd*85+gb_digit(85);
  304.     *(cur_pix+3)=dd&0xff;@+dd=(dd>>8)&0xffffff;
  305.     *(cur_pix+2)=dd&0xff;@+dd>>=8;
  306.     *(cur_pix+1)=dd&0xff;@+*cur_pix=dd>>8;
  307.     if (--j==0) gb_newline(),j=15;
  308.   }
  309.   *(cur_pix+1)=dd&0xff;@+*cur_pix=dd>>8;@+gb_newline();
  310. }
  311. @ @<Private var...@>=
  312. static long in_row[MAX_N];
  313. @* Planar graphs. We can obtain a large family of planar graphs based on
  314. digitizations of Mona Lisa by using the following simple scheme: Each matrix
  315. of pixels defines a set of connected regions containing pixels of the same
  316. value. (Two pixels are considered adjacent if they share an edge.)
  317. These connected regions are taken to be vertices of an undirected graph;
  318. two vertices are adjacent if the corresponding regions have at least
  319. one pixel edge in common.
  320. We can also state the construction another way. If we take any planar
  321. graph and collapse two adjacent vertices, we obtain another planar
  322. graph. Suppose we start with the planar graph having $mn$ vertices
  323. $[k,l]$ for $0le k<m$ and $0le l<n$, where $[k,l]$ is adjacent to
  324. $[k,l-1]$ when $l>0$ and to $[k-1,l]$ when $k>0$. Then we can attach
  325. pixel values to each vertex, after which we can repeatedly collapse
  326. adjacent vertices whose pixel values are equal. The resulting planar
  327. graph is the same as the graph of connected regions that was described
  328. in the previous paragraph.
  329. The subroutine call |plane_lisa(m,n,d,m0,m1,n0,n1,d0,d1)| constructs
  330. the planar graph associated with the digitization produced by |lisa|.
  331. The description of |lisa|, given earlier, explains the significance of
  332. parameters |m|, |n|, |d|, |m0|, |m1|, |n0|, |n1|, |d0|, and |d1|. There will
  333. be at most $mn$ vertices, and the graph will be simply an $mtimes n$
  334. grid unless |d| is small enough to permit adjacent pixels to have
  335. equal values. The graph will also become rather trivial if |d| is
  336. too small.
  337. Utility fields |first_pixel| and |last_pixel| give, for each vertex,
  338. numbers of the form $k*n+l$, identifying the topmost/leftmost
  339. and bottommost/rightmost positions $[k,l]$ in the region corresponding
  340. to that vertex. Utility fields |matrix_rows| and |matrix_cols| in
  341. the |Graph| record contain the values of |m| and~|n|; thus, in particular,
  342. the value of |n| needed to decompose |first_pixel| and |last_pixel| into
  343. individual coordinates can be found in |g->matrix_cols|.
  344. The original pixel value of a vertex is placed into its |pixel_value|
  345. utility field.
  346. @d pixel_value x.I
  347. @d first_pixel y.I
  348. @d last_pixel z.I
  349. @d matrix_rows uu.I
  350. @d matrix_cols vv.I
  351. @p Graph *plane_lisa(m,n,d,m0,m1,n0,n1,d0,d1)
  352.   unsigned long m,n; /* number of rows and columns desired */
  353.   unsigned long d; /* maximum value desired */
  354.   unsigned long m0,m1; /* input will be from rows $[|m0|,.,.,|m1|)$ */
  355.   unsigned long n0,n1; /* and from columns $[|n0|,.,.,|n1|)$ */
  356.   unsigned long d0,d1; /* lower and upper threshold of raw pixel scores */
  357. {@+@<Local variables for |plane_lisa|@>@;@#
  358.   init_area(working_storage);
  359.   @<Figure out the number of connected regions, |regs|@>;
  360.   @<Set up a graph with |regs| vertices@>;
  361.   @<Put the appropriate edges into the graph@>;
  362. trouble: gb_free(working_storage);
  363.   if (gb_trouble_code) {
  364.     gb_recycle(new_graph);
  365.     panic(alloc_fault); /* oops, we ran out of memory somewhere back there */
  366.   }
  367.   return new_graph;
  368. }
  369. @ @<Local variables for |plane_lisa|@>=
  370. Graph *new_graph; /* the graph constructed by |plane_lisa| */
  371. register long j,k,l; /* all-purpose indices */
  372. Area working_storage; /* tables needed while |plane_lisa| does its thinking */
  373. long *a; /* the matrix constructed by |lisa| */
  374. long regs=0; /* number of vertices generated so far */
  375. @ @<gb_lisa.h@>=
  376. #define pixel_value @tquad@> x.I /* definitions for the header file */
  377. #define first_pixel @tquad@> y.I
  378. #define last_pixel @tquad@> z.I
  379. #define matrix_rows @tquad@> uu.I
  380. #define matrix_cols @tquad@> vv.I
  381. @ The following algorithm for counting the connected regions considers
  382. the array elements |a[k,l]| to be linearly ordered as they appear
  383. in memory. Thus we can speak of the $n$ elements preceding a given
  384. element |a[k,l]|, if $k>0$; these are the elements |a[k,l-1]|, dots,
  385. |a[k,0]|, |a[k-1,n-1]|, dots, |a[k-1,l]|. These $n$ elements appear
  386. in $n$ different columns.
  387. During the algorithm, we move through the array from bottom right
  388. to top left, maintaining an auxiliary table $langle f[0],ldots,f[n-1]
  389. rangle$ with the following significance: Whenever two of the
  390. $n$ elements preceding our current position $[k,l]$ are connected to
  391. each other by a sequence of pixels with equal value, where the connecting
  392. links do not involve pixels more than $n$ steps before our current
  393. position, those elements will be linked together in the $f$ array.
  394. More precisely, we will have $f[c_1]=c_2$, dots, $f[c_{j-1}]=c_j$,
  395. and $f[c_j]=c_j$, when there are $j$ equivalent elements in columns
  396. $c_1$, dots,~$c_j$. Here $c_1$ will be the ``last'' column and
  397. $c_j$ the ``first,'' in wraparound order; each element with $f[c]ne c$
  398. points to an earlier element.
  399. The main function of the |f| table is to identify the topmost/leftmost
  400. pixel of a region. If we are at position |[k,l]| and if we find $f[l]=l$
  401. while $a[k-1,l]ne a[k,l]$, there is no way to connect |[k,l]| to
  402. earlier positions, so we create a new vertex for it.
  403. We also change the |a| matrix, to facilitate another algorithm
  404. below. If position |[k,l]| is the topmost/leftmost pixel of a region,
  405. we set |a[k,l]=-1-a[k,l]|; otherwise we set |a[k,l]=f[l]|, the column of
  406. a preceding element belonging to the same region.
  407. @<Figure out the number...@>=
  408. a=lisa(m,n,d,m0,m1,n0,n1,d0,d1,working_storage);
  409. if (a==NULL) return NULL; /* |panic_code| has been set by |lisa| */
  410. sscanf(lisa_id,"lisa(%lu,%lu,",&m,&n); /* adjust for defaults */
  411. f=gb_typed_alloc(n,unsigned long,working_storage);
  412. if (f==NULL) {
  413.   gb_free(working_storage); /* recycle the |a| matrix */
  414.   panic(no_room+2); /* there's no room for the |f| vector */
  415. }
  416. @<Pass over the |a| matrix from bottom right to top left, looking
  417.   for the beginnings of connected regions@>;
  418. @ @<Local variables for |plane_lisa|@>=
  419. unsigned long *f; /* beginning of array |f|;
  420.                      $f[j]$ is the column of an equivalent element */
  421. long *apos; /* the location of |a[k,l]| */
  422. @ We maintain a pointer |apos| equal to |&a[k,l]|, so that
  423. |*(apos-1)=a[k,l-1]| and |*(apos-n)=a[k-1,l]| when $l>0$ and $k>0$.
  424. The loop that replaces $f[j]$ by $j$ can cause this algorithm to
  425. take time $mn^2$. We could improve the worst case by using path
  426. compression, but the extra complication is rarely worth the trouble.
  427. @<Pass over the |a| matrix from bottom right to top left, looking
  428.   for the beginnings of connected regions@>=
  429. for (k=m, apos=a+n*(m+1)-1; k>=0; k--)
  430.   for (l=n-1; l>=0; l--,apos--) {
  431.     if (k<m) {
  432.       if (k>0&&*(apos-n)==*apos) {
  433.         for (j=l; f[j]!=j; j=f[j]) ; /* find the first element */
  434.         f[j]=l; /* link it to the new first element */
  435.         *apos=l;
  436.       }@+else if (f[l]==l) *apos=-1-*apos,regs++; /* new region found */
  437.         else *apos=f[l];
  438.     }
  439.     if (k>0&&l<n-1&&*(apos-n)==*(apos-n+1)) f[l+1]=l;
  440.     f[l]=l;
  441.   }
  442. @ @<Set up a graph with |regs| vertices@>=
  443. new_graph=gb_new_graph(regs);
  444. if (new_graph==NULL)
  445.   panic(no_room); /* out of memory before we're even started */
  446. sprintf(new_graph->id,"plane_%s",lisa_id);
  447. strcpy(new_graph->util_types,"ZZZIIIZZIIZZZZ");
  448. new_graph->matrix_rows=m;
  449. new_graph->matrix_cols=n;
  450. @ Now we make another pass over the matrix, this time from top left
  451. to bottom right. An auxiliary vector of length |n| is once again
  452. sufficient to tell us when one region is adjacent to a previous one.
  453. In this case the vector is called |u|, and it contains pointers to
  454. the vertices in the $n$ positions before our current position.
  455. We assume that a pointer to a |Vertex| takes the same amount of
  456. memory as an |unsigned long|, hence |u| can share the space formerly
  457. occupied by~|f|; if this is not the case, a system-dependent
  458. change should be made here.
  459. @^system dependencies@>
  460. The vertex names are simply integers, starting with 0.
  461. @<Put the appropriate edges into the graph@>=
  462. regs=0;
  463. u=(Vertex**)f;
  464. for (l=0;l<n;l++) u[l]=NULL;
  465. for (k=0,apos=a,aloc=0;k<m;k++)
  466.   for (l=0;l<n;l++,apos++,aloc++) {
  467.     w=u[l];
  468.     if (*apos<0) {
  469.       sprintf(str_buf,"%ld",regs);
  470.       v=new_graph->vertices+regs;
  471.       v->name=gb_save_string(str_buf);
  472.       v->pixel_value=-*apos-1;
  473.       v->first_pixel=aloc;
  474.       regs++;
  475.     }@+else v=u[*apos];
  476.     u[l]=v;
  477.     v->last_pixel=aloc;
  478.     if (gb_trouble_code) goto trouble;
  479.     if (k>0 && v!=w) adjac(v,w);
  480.     if (l>0 && v!=u[l-1]) adjac(v,u[l-1]);
  481.   }
  482. @ @<Local variables for |pl...@>=
  483. Vertex **u; /* table of vertices for previous $n$ pixels */
  484. Vertex *v; /* vertex corresponding to position |[k,l]| */
  485. Vertex *w; /* vertex corresponding to position |[k-1,l]| */
  486. long aloc; /* $k*n+l$ */
  487. @ The |adjac| routine makes two vertices adjacent, if they aren't already.
  488. A faster way to recognize duplicates would probably speed things up.
  489. @<Private sub...@>=
  490. static void adjac(u,v)
  491.   Vertex *u,*v;
  492. {@+Arc *a;
  493.   for (a=u->arcs;a;a=a->next)
  494.     if (a->tip==v) return;
  495.   gb_new_edge(u,v,1L);
  496. }
  497. @* Bipartite graphs. An even simpler class of Mona-Lisa-based graphs
  498. is obtained by considering the |m| rows and |n| columns to be individual
  499. vertices, with a row adjacent to a column if the associated pixel value
  500. is sufficiently large or sufficiently small. All edges have length~1.
  501. The subroutine call |bi_lisa(m,n,m0,m1,n0,n1,thresh,c)| constructs
  502. the bipartite graph corresponding to the $mtimes n$
  503. digitization produced by |lisa|, using parameters |(m0,m1,n0,n1)| to
  504. define a rectangular subpicture as described earlier.
  505. The threshold parameter |thresh| should be between 0 and~65535.
  506. If the pixel value in row |k| and column |l| is at least |thresh/65535| of
  507. its maximum, vertices |k| and~|l| will be adjacent.
  508. If |c!=0|, however, the convention is reversed; vertices are then
  509. adjacent when the corresponding pixel value is {sl smaller/} than
  510. |thresh/65535|. Thus adjacencies come from ``light'' areas of
  511. da Vinci's painting when |c=0| and from ``dark'' areas when |c!=0|. There
  512. are |m+n| vertices and up to $mtimes n$ edges.
  513. The actual pixel value is recorded in utility field |b.I| of each arc,
  514. and scaled to be in the range $[0,65535]$.
  515. @p Graph *bi_lisa(m,n,m0,m1,n0,n1,thresh,c)
  516.   unsigned long m,n; /* number of rows and columns desired */
  517.   unsigned long m0,m1; /* input will be from rows $[|m0|,.,.,|m1|)$ */
  518.   unsigned long n0,n1; /* and from columns $[|n0|,.,.,|n1|)$ */
  519.   unsigned long thresh; /* threshold defining adjacency */
  520.   long c; /* should we prefer dark pixels to light pixels? */
  521. {@+@<Local variables for |bi_lisa|@>@;@#
  522.   init_area(working_storage);
  523.   @<Set up a bipartite graph with |m+n| vertices@>;
  524.   @<Put the appropriate edges into the bigraph@>;
  525.   gb_free(working_storage);
  526.   if (gb_trouble_code) {
  527.     gb_recycle(new_graph);
  528.     panic(alloc_fault); /* oops, we ran out of memory somewhere back there */
  529.   }
  530.   return new_graph;
  531. }
  532. @ @<Local variables for |bi_lisa|@>=
  533. Graph *new_graph; /* the graph constructed by |bi_lisa| */
  534. register long k,l; /* all-purpose indices */
  535. Area working_storage; /* tables needed while |bi_lisa| does its thinking */
  536. long *a; /* the matrix constructed by |lisa| */
  537. long *apos; /* the location of |a[k,l]| */
  538. register Vertex *u,*v; /* current vertices of interest */
  539. @ @<Set up a bipartite graph...@>=
  540. a=lisa(m,n,65535L,m0,m1,n0,n1,0L,0L,working_storage);
  541. if (a==NULL) return NULL; /* |panic_code| has been set by |lisa| */
  542. sscanf(lisa_id,"lisa(%lu,%lu,65535,%lu,%lu,%lu,%lu",&m,&n,&m0,&m1,&n0,&n1);
  543. new_graph=gb_new_graph(m+n);
  544. if (new_graph==NULL)
  545.   panic(no_room); /* out of memory before we're even started */
  546. sprintf(new_graph->id,"bi_lisa(%lu,%lu,%lu,%lu,%lu,%lu,%lu,%c)",
  547.    m,n,m0,m1,n0,n1,thresh,c?'1':'0');
  548. new_graph->util_types[7]='I'; /* enable field |b.I| */
  549. mark_bipartite(new_graph,m);
  550. for (k=0,v=new_graph->vertices;k<m;k++,v++) {
  551.   sprintf(str_buf,"r%ld",k); /* row vertices are called |"r0"|, |"r1"|, etc. */
  552.   v->name=gb_save_string(str_buf);
  553. }
  554. for (l=0;l<n;l++,v++) {
  555.   sprintf(str_buf,"c%ld",l); /* column vertices are called |"c0"|,
  556.                                             |"c1"|, etc. */
  557.   v->name=gb_save_string(str_buf);
  558. }
  559. @ Since we've called |lisa| with |d=65535|, the determination of
  560. adjacency is simple.
  561. @<Put the appropriate edges into the bigraph@>=
  562. for (u=new_graph->vertices,apos=a;u<new_graph->vertices+m;u++)
  563.   for (v=new_graph->vertices+m;v<new_graph->vertices+m+n;apos++,v++) {
  564.     if (c?*apos<thresh:*apos>=thresh) {
  565.       gb_new_edge(u,v,1L);
  566.       u->arcs->b.I=v->arcs->b.I=*apos;
  567.     }
  568.   }
  569. @* Index. As usual, we close with an index that
  570. shows where the identifiers of \{gb_lisa} are defined and used.