apclust.c

00001 /*
00002 
00003 $Id: apclust.c,v 1.1 2005/09/13 13:25:26 jim Exp $
00004 
00005 */
00006 
00007 #include <stdio.h>
00008 #include <cpl.h>
00009 
00010 #include "imcore.h"
00011 #include "util.h"
00012 
00013 static void minmax_xy(int, plstruct *, int *, int *, int *, int *);
00014 
00015 void apclust(ap_t *ap, int np, plstruct *plstr) {
00016 
00017     int i,i1,loop,ik;
00018     int is;    /* parent name for image in this slice */
00019     int ip;    /* parent name for image on last line */
00020     int ib;    /* data block name */
00021     int k,j,ix1,ix2,iy1,iy2,nwork,nx,ny,kk;
00022     float i2compare,icompare;
00023     short int *work;
00024  
00025     /* A couple of useful things */
00026 
00027     i2compare = ap->thresh;
00028     icompare = i2compare * ap->multiply;
00029 
00030     /* Get the min and max positions. Create a raster with the IDs of the
00031        pixels in the pixel list (algorithm prefers data to be in a raster) */
00032 
00033     minmax_xy(np,plstr,&ix1,&ix2,&iy1,&iy2);
00034     nx = ix2 - ix1 + 1;
00035     ny = iy2 - iy1 + 1;
00036     nwork = nx*ny;
00037     work = cpl_malloc(nwork*sizeof(*work));
00038     for (i = 0; i < nwork; i++)
00039         work[i] = -1;
00040     for (k = 0; k < np; k++) {
00041         i = plstr[k].x - 1;
00042         j = plstr[k].y - 1;
00043         kk = (j - iy1)*nx + i - ix1;
00044         work[kk] = k;
00045     }
00046 
00047     /* Now do the job */
00048             
00049     for (j = iy1; j <= iy2; j++) {
00050         for (i = ix1; i <= ix2; i++) {
00051             kk = (j - iy1)*nx + i - ix1;
00052             k = work[kk];
00053             if (k < 0) {
00054                 ap->lastline[i + 1] = 0;
00055             } else {
00056                 if (plstr[k].zsm > icompare) {
00057 
00058                     /* Pixel is above threshold, find which parent it belongs to. */
00059 
00060                     is = ap->lastline[i];       /* Parent last pixel this line */
00061                     ip = ap->lastline[i + 1];   /* Guess belongs to above line */
00062                     if (ip == 0) {
00063 
00064                         /* New parent, or, horizontal slice: */
00065 
00066                         if (is == 0) {
00067 
00068                             /* Ah - new parent. */
00069 
00070                             if (ap->ipstack > ap->maxpa*3/4) {
00071                                 for (ik = 0; ik < ap->maxpa*3/8; ik++)
00072                                     apfu(ap);
00073                             }
00074                             ip = ap->pstack[ap->ipstack++];
00075                             ap->parent[ip].first = ap->bstack[ap->ibstack];
00076                             ap->parent[ip].pnop = 0;
00077                             ap->parent[ip].pnbp = 0;
00078                             ap->parent[ip].growing = 0;
00079                             if (j == 0)
00080 
00081                                 /* It touches first line: */
00082 
00083                                 ap->parent[ip].touch = 1;
00084                             else
00085                                 ap->parent[ip].touch = 0;
00086 
00087                             /* For hunt thru list for terminates: */
00088 
00089                             if (ip > ap->maxip)
00090                                 ap->maxip = ip;
00091                         } else {
00092 
00093                             /* Slice with no vertical join: */
00094 
00095                             ip = is;
00096                         }
00097                     } else if ((ip > 0 && is > 0) && (ip != is)) {
00098 
00099                         /* merge: Join linked lists: */
00100 
00101                         ap->blink[ap->parent[ip].last] = ap->parent[is].first;
00102 
00103                         /* Copy `last block': */
00104 
00105                         ap->parent[ip].last = ap->parent[is].last;
00106                         ap->parent[ip].pnop += ap->parent[is].pnop;
00107                         ap->parent[ip].pnbp += ap->parent[is].pnbp;
00108 
00109                         /* Fix `lastline' correlator array: */
00110 
00111                         ib = ap->parent[is].first;
00112                         loop = 1;
00113                         while (loop) {
00114                             i1 = ap->plessey[ib].x;
00115                             if (ap->lastline[i1 + 1] == is)
00116                                 ap->lastline[i1 + 1] = ip;
00117                             if (ap->parent[is].last == ib)
00118                                 loop = 0;
00119                             else
00120                                 ib = ap->blink[ib];
00121                         }
00122 
00123                         /* Mark parent inactive: */
00124 
00125                         ap->parent[is].pnop = -1;
00126                         ap->parent[is].pnbp = -1;
00127 
00128                         /* return name to stack: */
00129 
00130                         ap->pstack[--ap->ipstack] = is;
00131                     }
00132 
00133                     /* Add in pixel to linked list: */
00134 
00135                     ib = ap->bstack[ap->ibstack++];
00136 
00137                     /* Patch forward link into last data block: */
00138 
00139                     if (ap->parent[ip].pnop > 0)
00140                         ap->blink[ap->parent[ip].last] = ib;
00141 
00142                     /* Remember last block in chain: */
00143 
00144                     ap->parent[ip].last = ib;
00145 
00146                     /* Store the data: */
00147 
00148                     ap->plessey[ib].x = i;
00149                     ap->plessey[ib].y = j;
00150                     ap->plessey[ib].z = plstr[k].z;
00151                     ap->plessey[ib].zsm = plstr[k].zsm;
00152 
00153                     /* increment active count: */
00154 
00155                     ap->parent[ip].pnop++;
00156 
00157                     /* remember which parent this pixel was for next line: */
00158 
00159                     ap->lastline[i + 1] = ip;
00160 
00161                 } else {
00162 
00163                     /* Pixel was below threshold, mark lastline: */
00164 
00165                     ap->lastline[i + 1] = 0;
00166                 }
00167             }
00168         }
00169     }
00170 
00171     /* Check for images touching left & right edges:
00172        OR the touch flag with 2 for left, 4 for right: */
00173 
00174     if(ap->lastline[1] > 0 )
00175         ap->parent[ap->lastline[1]].touch |= 2;
00176     if(ap->lastline[ap->lsiz] > 0)
00177         ap->parent[ap->lastline[ap->lsiz]].touch |= 4;
00178     cpl_free(work);
00179 }
00180 
00181 static void minmax_xy(int np, plstruct *plstr, int *ix1, int *ix2, int *iy1,
00182                       int *iy2) {
00183     int i;
00184 
00185     /* Get the minmax of the positions of the pixels in a plstruct.  Take
00186        1 away from each position so that it runs from 0 rather than 1 */
00187 
00188     *ix1 = plstr[0].x - 1;
00189     *ix2 = plstr[0].x - 1;
00190     *iy1 = plstr[0].y - 1;
00191     *iy2 = plstr[0].y - 1; 
00192     for (i = 1; i < np; i++) {
00193         *ix1 = MIN(*ix1,plstr[i].x - 1);
00194         *ix2 = MAX(*ix2,plstr[i].x - 1);
00195         *iy1 = MIN(*iy1,plstr[i].y - 1);
00196         *iy2 = MAX(*iy2,plstr[i].y - 1);
00197     }
00198 }
00199 
00200 
00201 /*
00202 
00203 $Log: apclust.c,v $
00204 Revision 1.1  2005/09/13 13:25:26  jim
00205 Initial entry after modifications to make cpl compliant
00206 
00207 
00208 */

Generated on Sat Apr 6 04:03:06 2013 for VIRCAM Pipeline by  doxygen 1.5.1