home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_6.2 / XVORA / KNNS.C < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  17.4 KB  |  687 lines
  1. /*
  2.  * k-N(earest)N(eighbor) S(hells)
  3.  *
  4.  * construct kNN shells on the basis of Voronoi diagram output;
  5.  *
  6.  */
  7. #include <stdio.h>
  8. #include <stdlib.h>
  9. #include "vora.h"
  10.  
  11. #undef    DEBUG
  12. #undef    KNN_DEBUG
  13. #undef MN_DEBUG
  14. #undef    KNN_ARRAY_VER
  15.  
  16. /* globals */
  17. extern int KNN_MAX;
  18. int DISPL_KNNS = 0;             /* when set, display full KNN shells */
  19.  
  20.  
  21. /*
  22.  * for a given Site, evaluate the average, mn, of the nof NN of its nnn NN
  23.  */
  24. float
  25. eval_mn (struct vSite *vS, struct vSite *vsa)
  26. {
  27.   float mn;
  28.   int inn;
  29.  
  30.   mn = (float) 0.0;
  31.   for (inn = 0; inn < vS->nnn; inn++) {
  32. #ifdef MN_DEBUG
  33.     printf ("\t %d (with %d NNs)\n",
  34.             vS->nnsi[inn], (vsa + (vS->nnsi[inn]))->nnn);
  35. #endif
  36.     mn += (float) (vsa + (vS->nnsi[inn]))->nnn;
  37.   }
  38.   mn /= (float) vS->nnn;
  39.  
  40.   return (mn);
  41. }
  42.  
  43.  
  44. /*
  45.  * evaluate, for each site, the average, m, of the nof NN of its NN;
  46.  * ignore sites which have been flagged (eFlag, aoiFlag) out-of-bounds
  47.  *
  48.  * --> Aboav-Weaire law for cellular structures 
  49.  * K.B.Lauritsen, C. Moukarzel & H.J.Herrmann, J. Physique I 3, 1941 (1993)
  50.  */
  51. int
  52. m_of_n (int *ndn, float **m_n, int ns, struct vSite *vsa)
  53. {
  54.   int is, nas;
  55.   int cnn;
  56.   float mn_av;
  57.   struct vSite *vS;
  58.  
  59.   nas = 0;
  60.   for (is = 0; is < ns; is++) {
  61.     vS = (vsa + is);
  62.     if ((vS->eFlag == UnSet) &&
  63.         (vS->aoiFlag == UnSet)) {
  64.       cnn = vS->nnn;
  65.       mn_av = *(m_n[cnn] + ndn[cnn]) = eval_mn (vS, vsa);
  66. #ifdef MN_DEBUG
  67.       printf ("\nM_OF_N: process site: %d\n", vS->sitenbr);
  68.       printf ("\tmn_av: %f, ndn[%d]: %d, *(mn[]+ndn[]): %f\n", mn_av, cnn, ndn[cnn], *(m_n[cnn] + ndn[cnn]));
  69. #endif
  70.       ndn[cnn]++;
  71.       nas++;
  72.     }
  73.   }
  74.   return (nas);
  75. }
  76.  
  77.  
  78. /*
  79.  * evaluate the average area, Ak, of domains in the kNN shells of
  80.  * the current Site
  81.  */
  82. void
  83. eval_aka (struct vSite *vsa, int kNN, double *aka, struct linklist *sha)
  84. {
  85.   int k;
  86.  
  87.   for (k = 0; k < kNN; k++) {   /* loop over k-th NN shells */
  88.     *(aka + k) = xak_kNN_list (vsa, sha + k, k);
  89.   }
  90. }
  91.  
  92.  
  93.  
  94.  
  95. /*
  96.  * evaluate, for the current Site, the product c(0)c(k)/nkNN, as a function
  97.  * of k: c denotes topological charge, n-6, of a given site, nkNN denotes
  98.  * the number of sites in the k-th NN shell; these products represent the
  99.  * correlation function of topological charge as a function of topological
  100.  * distance, k;
  101.  */
  102. int
  103. eval_tctc (struct vSite *vsa, int kNN, double *tctc, struct linklist *sha)
  104. {
  105.   int k;
  106.   double c0, ck_bar;
  107.  
  108.  
  109.   c0 = xc0 (vsa, sha + 0);      /* topol charge of ref Site site */
  110.   *(tctc + 0) = c0 * c0;
  111.   for (k = 1; k < kNN; k++) {   /* loop over k-th NN shells */
  112.     ck_bar = xck_kNN_list (vsa, sha + k, k);
  113.     *(tctc + k) = c0 * ck_bar;
  114.   }
  115.   return ((int) c0);
  116. }
  117.  
  118.  
  119. /*
  120.  * evaluate the fractal measure, c, of a set of kNN shells, associated
  121.  * with a given site, such that N_k = c*k**2, 0<=k<=kNN, N_k representing
  122.  * the nof sites contained within the k-th NN shell, i.e., 
  123.  * N_k = sum from {j=0 to j=k} N_j;
  124.  */
  125. double
  126. eval_frms (int kNN, struct linklist *sha)
  127. {
  128.   int k;
  129.   int nk;
  130.   double frms;
  131.  
  132.   nk = 1;
  133.   for (k = 1; k < kNN; k++) {   /* <= ??? */
  134.     nk += ll_length (sha + k);
  135.   }
  136.   frms = (double) nk / (double) ((kNN - 1) * (kNN - 1));
  137.   return (frms);
  138. }
  139.  
  140.  
  141.  
  142. /*
  143.  * compare list entries 
  144.  *   key: Site index
  145.  */
  146. int
  147. cmpSiteInd (struct nnSite *oldSite, struct nnSite *newSite)
  148. {
  149.   return (SIGN (oldSite->sitenbr - newSite->sitenbr));
  150. }
  151.  
  152. /*
  153.  * init (empty) kNN shell lists
  154.  */
  155. void
  156. init_sk (struct linklist *knnsll)
  157. {
  158.   llzero (knnsll);              /* init empty list */
  159.   llsetsize (sizeof (struct nnSite), knnsll);
  160.  
  161.   llsetmatch (cmpSiteInd, knnsll);
  162. }
  163.  
  164.  
  165.  
  166. /*
  167.  * scan list, comparing index of newSite with that of Sites already 
  168.  * in current list; add Site in order of increasing indices; 
  169.  * employs matching function, here cmpSiteInd();
  170.  * (see also: llsearch() in llist.c )
  171.  */
  172. Boolean
  173. scanSiteInd (struct nnSite *newSite, struct linklist *list)
  174. {
  175.   Boolean MatchStatus = False;
  176.  
  177.   llhead (list);
  178.   do {
  179.     if (((*list->match) (list->clp->item, newSite)) > 0)
  180.       MatchStatus = True;
  181.  
  182.   } while ((MatchStatus != True) && (llnext (list) == True));
  183.  
  184.   return (MatchStatus);
  185. }
  186.  
  187.  
  188. /*
  189.  * insert new site into current list, maintaining the list sorted
  190.  * in order of increasing nnSite indices
  191.  */
  192. void
  193. expand_kNN_list (struct nnSite *nnnS, struct linklist *list)
  194. {
  195.   Boolean MatchStatus;
  196.  
  197.   llhead (list);
  198.   if ((MatchStatus = scanSiteInd (nnnS, list)) == True)
  199.     lladdprev ((char *) nnnS, list);
  200.   else
  201.     lladd ((char *) nnnS, list);
  202. }
  203.  
  204.  
  205. /*
  206.  * remove duplicate entries from kNN list which is assumed to contain
  207.  * kNN site indices in increasing order
  208.  */
  209. void
  210. contract_kNN_list (struct linklist *list)
  211. {
  212.   struct nnSite *nnS, *nnSpp;
  213.  
  214.  
  215.   if (ll_length (list) < 2) {
  216. #ifdef KNN_DEBUG
  217.     printf ("...list contains less than two entries\n");
  218. #endif
  219.   }
  220.   else {
  221.     llhead (list);
  222.     nnS = (struct nnSite *) list->clp->item;
  223.     while (llnext (list) == True) {
  224.       nnSpp = (struct nnSite *) list->clp->item;
  225.       if (nnS->sitenbr == nnSpp->sitenbr) {
  226.         lldelete (list);
  227.         nnS = (struct nnSite *) list->clp->item;
  228.       }
  229.       else {
  230.         nnS = nnSpp;
  231.       }
  232.     }
  233.   }
  234. }
  235.  
  236.  
  237. /*
  238.  * display list
  239.  */
  240. void
  241. show_kNN_list (struct linklist *list, int k)
  242. {
  243.   int is;
  244.  
  245.   printf ("SHOW_KNN_LIST: %2d-NN shell:", k);
  246.   llhead (list);
  247.   if (ll_length (list) != 0) {  /* list not empty */
  248.     is = 0;
  249.     do {
  250.       printf (" %3d ",
  251.               ((struct nnSite *) list->clp->item)->sitenbr);
  252.       is++;
  253.  
  254.       if ((is + 1) % 10 == 0)
  255.         printf ("\n");
  256.     } while (llnext (list) == True);
  257.     printf ("\n");
  258.   }
  259.   else
  260.     printf ("  list empty\n");
  261. }
  262.  
  263.  
  264.  
  265.  
  266. /*
  267.  * extract average area of k-th NN shells from corresponding lists
  268.  */
  269. double
  270. xak_kNN_list (struct vSite *vsa, struct linklist *list, int k)
  271. {
  272.   int is, lk;
  273.   double ak = 0.0;
  274.  
  275.   llhead (list);
  276.   if ((lk = ll_length (list)) != 0) {  /* list not empty */
  277.     do {
  278.       is = ((struct nnSite *) list->clp->item)->sitenbr;
  279.       ak += (double) ((vsa + is)->area);
  280.  
  281.     } while (llnext (list) == True);
  282.     ak /= (double) lk;
  283.   }
  284.   else
  285.     printf ("  list for %d-th NN shell empty\n", k);
  286.  
  287.   return (ak);
  288. }
  289.  
  290.  
  291. /*
  292.  * extract topol charge of ref Site (stored at top of list sha+0)
  293.  */
  294. double
  295. xc0 (struct vSite *vsa, struct linklist *list)
  296. {
  297.   int is, lk;
  298.   double c0 = 100.0;
  299.  
  300.   llhead (list);
  301.   if ((lk = ll_length (list)) != 0) {  /* list not empty */
  302.     is = ((struct nnSite *) list->clp->item)->sitenbr;
  303.     c0 = (double) ((vsa + is)->nnn - 6);
  304.   }
  305.   else
  306.     printf ("  list empty, something wrong\n");
  307.  
  308.   return (c0);
  309. }
  310.  
  311.  
  312. /*
  313.  * extract average topological charge of k-th NN shells 
  314.  * from corresponding lists
  315.  */
  316. double
  317. xck_kNN_list (struct vSite *vsa, struct linklist *list, int k)
  318. {
  319.   int is, lk;
  320.   double tck = 0.0;
  321.  
  322.   llhead (list);
  323.   if ((lk = ll_length (list)) != 0) {  /* list not empty */
  324.     do {
  325.       is = ((struct nnSite *) list->clp->item)->sitenbr;
  326.       tck += (double) ((vsa + is)->nnn - 6);
  327.  
  328.     } while (llnext (list) == True);
  329.     tck /= (double) lk;
  330.   }
  331.   else
  332.     printf ("  list for %d-th NN shell empty\n", k);
  333.  
  334.   return (tck);
  335. }
  336.  
  337.  
  338. /*
  339.  * copy vSite information to root of kNN shell lists
  340.  */
  341. void
  342. init_sk0 (struct vSite *vsa, int sitenbr, struct linklist *list)
  343. {
  344.   struct nnSite nnS, *pnnS = &nnS;
  345.  
  346.   (vsa + sitenbr)->status = InActive;
  347.   pnnS->sitenbr = sitenbr;
  348.  
  349.   expand_kNN_list (pnnS, list);
  350. }
  351.  
  352.  
  353.  
  354. /*
  355.  * set Status of all Sites to Active
  356.  */
  357. void
  358. SetActive (struct vSite *vsa, int ns)
  359. {
  360.   int is;
  361.  
  362.   for (is = 0; is < ns; is++)
  363.     (vsa + is)->status = Active;
  364. }
  365.  
  366.  
  367. /*
  368.  * extract, from Voronoi diagram, kNN shells of a given site, 
  369.  * in order of increasing k: k=0 refers to the Site itself;
  370.  * the NN shells are constructed on the basis of an array of Voronoi sites
  371.  * which gives site coordinates, site index and NN site indices;
  372.  *
  373.  * the desired list of NN shells is implemented as a list of lists, for
  374.  * each site in the set; for each k, there is a list of indices giving the
  375.  * sites of the corresponding shell; for given ROOT site, each site is
  376.  * uniquely assigned to one shell: that is, as k increases, only those
  377.  * site indices area added to the current shell which have not been 
  378.  * assigned to any previous shell; once assigned, a site is given InActive
  379.  * status; the recursion comes to an end when either a maximum k is reached
  380.  * or else, if there are no more Active sites to be added to new shells
  381.  */
  382. struct linklist *
  383. kNNs (struct vSite *vsa, struct linklist *plist, struct linklist *clist, int k)
  384. {
  385.   struct nnSite *pnnS, nnS, *cnnS = &nnS;
  386.   struct vSite *vS;
  387.   int inn;
  388.   Boolean Descend_Status;
  389.  
  390.  
  391. #ifdef DEBUG
  392.   printf ("\nKNNS:\n");
  393.   printf ("...pointers:\n");
  394.   printf ("\tvsa: %p, plist: %p, clist: %p\n", vsa, plist, clist);
  395.   printf ("\textract_kNN: %p\n", kNNs);
  396.   printf ("...new parent list length: %3d\n", ll_length (plist));
  397.   show_kNN_list (plist, k);
  398. #endif
  399.  
  400.   llhead (plist);               /* parent list */
  401.   llhead (clist);               /* child list */
  402.   Descend_Status = False;
  403.   do {                          /* step through parent list */
  404.     pnnS = (struct nnSite *) plist->clp->item;
  405.     vS = vsa + (pnnS->sitenbr);
  406.  
  407.     for (inn = 0; inn < vS->nnn; inn++) {  /* loop over NN ind of vS */
  408.       if ((vsa + (vS->nnsi[inn]))->status == Active) {
  409.  
  410.         cnnS->sitenbr = (vS->nnsi[inn]);
  411.         expand_kNN_list (cnnS, clist);
  412.  
  413.         if (Descend_Status == False)
  414.           Descend_Status = True;
  415.  
  416.         (vsa + (vS->nnsi[inn]))->status = InActive;
  417.       }
  418.     }
  419.   } while (llnext (plist) == True);
  420.  
  421.   if (Descend_Status == True) {
  422.     contract_kNN_list (clist);
  423.     return (clist);
  424.   }
  425.   else {
  426. #ifdef DEBUG
  427.     printf ("\nKNNS: no more Active sites\n");
  428. #endif
  429.     return (clist = (struct linklist *) NULL);
  430.   }
  431. }
  432.  
  433.  
  434.  
  435. /* 
  436.  * code to construct, on the basiis of the Voronoi diagram, the k-NN shells 
  437.  * for each site in a point set 
  438.  *
  439.  * the function construct_kNNs() handles the kNN shells for a single
  440.  * site at a time, and evaluates quantities of interest before going
  441.  * on to the next site; that is, a single array of shells, sk, properly
  442.  * dimensioned to KNN_MAX, is repeatedly used to store the kNN 
  443.  * shells of successively examined sites
  444.  *
  445.  * a note on the data structures:
  446.  *   sk is an array of pointers to struct kNNShell; one such array is
  447.  *   associated with each site, is <= ns, in the set; each entry in sk, 
  448.  *   sk[is], contains the size, kNN, of the actual nof shells built for
  449.  *   a given ``root'' site, and quantities of interest evaluated from that
  450.  *   shell, such as the fractal measure; 
  451.  *   for each root Site, an array, sha, of (maximally) KNN_MAX shells is
  452.  *   constructed (and then recycled for the following root Site): the k-th 
  453.  *   element, sha[k], of that array represents the k-th NN shell of the root 
  454.  *   Site in the form of a list of structs nnSite, containing, among other
  455.  *   entries, the indices (into the Site array vsa) identifying the k-th NN;
  456.  *
  457.  *   if memory consumption still presents a concern, this version could
  458.  *   be modified to dispense with the array, sk, and just store information
  459.  *   for a single Site at a time: in that case, evaluation of quantities 
  460.  *   such as the fractal measure would have to integrated in kNN construction
  461.  *   
  462.  *   a memory-consuming, but conceptually appealing version:
  463.  *   -------------------------------------------------------
  464.  *   for the akNNs version, sk[is] contains, in addition, an array of 
  465.  *   kNN lists; the k-th element, (sk+is)->sha[k], in any such array is 
  466.  *   a list containing the set of indices comprising the k-NN shell of 
  467.  *   the ``root'' or reference site;
  468.  *   sites which have been flagged (eFlag, aoiFlag) out-of-bounds
  469.  *   are not processed: their kNN shells (k>0) contain 0 elements
  470.  */
  471. double
  472. construct_kNNs (int ns, struct vSite *vsa, struct linklist *sha, struct kNNShell *sk)
  473. {
  474.   int is, k;
  475.   int nas;
  476.   double frms;
  477.   struct vSite *vS;
  478.  
  479.   struct linklist *plist, *clist;  /* ptr to parent(k) list */
  480.  
  481.  
  482. /*
  483.  * construct kNN shells, one site at a time
  484.  */
  485.   frms = 0.0;
  486.   nas = 0;
  487.   printf ("\tSite:");
  488.   for (is = 0; is < ns; is++) { /* loop over sites */
  489.     (sk + is)->sitenbr = (vsa + is)->sitenbr;
  490.  
  491.     if (is % 50 == 0)
  492.       printf ("...%3d", is);
  493.     if ((is + 1) == ns)
  494.       printf ("...%3d\n", is);
  495.  
  496. /* set status of all Sites to Active (must be repeated) */
  497.     SetActive (vsa, ns);
  498.  
  499. /* set up empty list structs */
  500.     for (k = 0; k < KNN_MAX; k++)
  501.       init_sk (sha + k);
  502.  
  503. /* loop over successive k-NN shells of current ``root'' site */
  504.     init_sk0 (vsa, (vsa + is)->sitenbr, sha + 0);
  505.     if (DISPL_KNNS == 1)
  506.       show_kNN_list (sha + 0, 0);
  507.  
  508.     if (((vS = (vsa + is))->eFlag == UnSet) &&
  509.         (vS->aoiFlag == UnSet)) {  /* Site out-of-bounds ? */
  510.       k = 0;
  511.       plist = (struct linklist *) (sha + 0);
  512.       do {
  513.         if (++k >= KNN_MAX) {
  514. #ifdef KNN_DEBUG
  515.           printf ("...reached max NN shell: ");
  516.           printf (" k: %d\n", k);
  517. #endif
  518.           plist = (struct linklist *) NULL;
  519.         }
  520.         else {
  521.           clist = (struct linklist *) (&sha[k]);
  522.           plist = kNNs (vsa, plist, clist, k);
  523.  
  524. #ifdef KNN_DEBUG
  525.           printf ("...%d-th shell completed\n", k);
  526.           for (j = 0; j <= k; j++)
  527.             show_kNN_list (sha + j, j);
  528. #endif
  529.         }
  530.       } while (plist != (struct linklist *) NULL);
  531.       (sk + is)->kNN = k;
  532.  
  533.       if (DISPL_KNNS == 1) {
  534.         printf ("...%d-NN shells for Site %d\n",
  535.                 k, vS->sitenbr);
  536.         for (k = 0; k < (sk + is)->kNN; k++)
  537.           show_kNN_list (sha + k, k);
  538.       }
  539.  
  540. /*
  541.  * employ kNN shells of current root Site to evaluate further
  542.  * quantities of interest; soter these in sk[is]
  543.  *
  544.  * fractal measure:
  545.  *   c, in Nk = c*k**2, where Nk denotes the number 
  546.  *   of domains in the k-th NN shell of a given site
  547.  */
  548.       (sk + is)->frms = eval_frms ((sk + is)->kNN, sha);
  549. #ifdef KNN_DEBUG
  550.       printf ("\n...fractal measure (from %d shells): %lf\n",
  551.               (sk + is)->kNN, (sk + is)->frms);
  552. #endif
  553.       frms += (sk + is)->frms;
  554.       nas++;
  555.  
  556.  
  557. /*
  558.  * average area, <Ak>, of domains in k-th NN shell, as function of k
  559.  */
  560.       eval_aka (vsa, (sk + is)->kNN, (sk + is)->aka, sha);
  561.  
  562. /*
  563.  * correlation function, tccf(k) = c(0)c(k), of topological charge,
  564.  * c=n-6, as function of topological distance, k, for given site;
  565.  */
  566.       (sk + is)->tc = eval_tctc (vsa, (sk + is)->kNN, (sk + is)->tctc, sha);
  567.     }
  568.     else {                      /* Site skipped */
  569.       (sk + is)->kNN = 0;
  570.     }
  571.   }
  572.   frms /= (double) nas;
  573.   return (frms);
  574. }
  575.  
  576.  
  577.  
  578. /* 
  579.  * code to construct k-NN shells for each site in a point set 
  580.  * on the basis of Voronoi diagram:
  581.  * this version, construct_akNNs() handles an entire array of shells,
  582.  * one for each site in the set (see loop over sites); this has the 
  583.  * disadvantage of requiring large amounts of memory, but has the
  584.  * advantage of making available all requisite data for subsequent
  585.  * computation of topological quantities
  586.  *
  587.  * 3/94: not completely debugged, not maintained
  588.  *
  589.  * a note on the data structures:
  590.  *    sk is an array of pointers to struct kNNShell; one such array is
  591.  *    associated with each site, is <= ns, in the set; each entry in sk, 
  592.  *    sk[is], contains the size, kNN, of the actual nof shells built for
  593.  *    a given ``root'' site, and an array of kNN lists; the k-th element, 
  594.  *    (sk+is)->sha[k], in any such array is a list containing the set 
  595.  *    of indices comprising the k-NN shell of the ``root'' site
  596.  */
  597. #ifdef KNN_ARRAY_VER
  598. struct kNNShell *
  599. construct_akNNs (int ns, struct vSite *vsa)
  600. {
  601.   int is, j, k, kNN;
  602.  
  603.   struct kNNShell *sk;          /* array of struct kNNShell */
  604.   struct linklist *plist, *clist;  /* ptr to parent(k) list */
  605.   struct linklist *sksi;
  606.  
  607. /*
  608.  * mem alloc
  609.  */
  610.   if ((sk = (struct kNNShell *) calloc (ns, sizeof (struct kNNShell))) == NULL)
  611.       fail_alloc ("sk", 1);
  612.  
  613.   for (is = 0; is < ns; is++) { /* loop over sites */
  614.  
  615. /*
  616.  * handle shell structure mem alloc
  617.  */
  618.     if (((sk + is)->sha = (struct linklist *) calloc (KNN_MAX, sizeof (struct linklist))) == NULL)
  619.         fail_alloc ("(sk+is)->sha", 1);
  620.  
  621.     sksi = (sk + is)->sha;      /* array of kNN lists */
  622.  
  623. /*
  624.  * printf("\n...size of struct linklist: %d\n", 
  625.  * sizeof(struct linklist));
  626.  * printf("...size of sksi: %d\n", sizeof(sksi));
  627.  */
  628.  
  629. /*
  630.  * set status of all Sites to Active (must be repeated)
  631.  */
  632.     SetActive (vsa, ns);
  633.  
  634. /* set up empty list structs */
  635.     printf ("\n...initialize empty k-NN shell lists\n");
  636.     for (k = 0; k < KNN_MAX; k++)
  637.       init_sk (&((sk + is)->sha[k]));
  638.  
  639. /* loop over successive k-NN shells of current ``root'' site */
  640.     init_sk0 (vsa, (vsa + is)->sitenbr, &((sk + is)->sha[0]));
  641.     show_kNN_list (&((sk + is)->sha[0]), 0);
  642.  
  643.     k = 0;
  644.     plist = &((sk + is)->sha[0]);
  645.     do {
  646. /*                      clist = (struct linklist *)(sksi+k); */
  647. /*                      
  648.  * printf("\n...pointers:\n");
  649.  * printf("\tsksi: %p, sksi+k: %p, ++sksi: %p\n",
  650.  * sksi, sksi+k, ++sksi);
  651.  * printf("\tplist: %p, clist: %p\n", plist, clist);
  652.  * --sksi;
  653.  */
  654.  
  655.       if (++k >= KNN_MAX) {
  656.         printf ("...reached max NN shell, k: %d\n", k);
  657.         plist = (struct linklist *) NULL;
  658.       }
  659.       else {
  660.         printf ("...new k: %d\n", k);
  661.         clist = &((sk + is)->sha[k]);
  662.  
  663.  
  664.         plist = kNNs (vsa, plist, clist, k);
  665.  
  666. #ifdef KNN_DEBUG
  667.         printf ("...%d-th shell completed\n", k);
  668.         for (j = 0; j <= k; j++)
  669.           show_kNN_list (&((sk + is)->sha[j]), j);
  670. #endif
  671.       }
  672.     } while (plist != (struct linklist *) NULL);
  673.     (sk + is)->kNN = kNN = k;
  674.  
  675.     printf ("...%d-NN shells for Site %d\n", k, (vsa + is)->sitenbr);
  676.     for (k = 0; k < (sk + is)->kNN; k++)
  677.       show_kNN_list (&((sk + is)->sha[k]), k);
  678.  
  679.     if (((sk + is)->sha = (struct linklist *) realloc ((sk + is)->sha, kNN * sizeof (struct linklist))) == NULL)
  680.         fail_alloc ("realloc (sk+is)->sha", 1);
  681.   }
  682.   for (is = 0; is < ns; is++)
  683.     free ((sk + is));
  684.   return (sk);
  685. }
  686. #endif
  687.