--- ray/src/common/interp2d.c 2013/02/14 19:57:10 2.10 +++ ray/src/common/interp2d.c 2013/02/16 00:41:12 2.13 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: interp2d.c,v 2.10 2013/02/14 19:57:10 greg Exp $"; +static const char RCSid[] = "$Id: interp2d.c,v 2.13 2013/02/16 00:41:12 greg Exp $"; #endif /* * General interpolation method for unstructured values on 2-D plane. @@ -23,8 +23,8 @@ static const char RCSid[] = "$Id: interp2d.c,v 2.10 20 * to reduce the influence of distant neighbors. This yields a * smooth interpolation regardless of how the sample points are * initially distributed. Evaluation is accelerated by use of - * a fast approximation to the atan2(y,x) function and an array - * of flags indicating where weights are (nearly) zero. + * a fast approximation to the atan2(y,x) function and a low-res + * map indicating where sample weights are significant. ****************************************************************/ #include @@ -41,6 +41,19 @@ typedef struct { float dm; /* distance measure in this direction */ } SAMPORD; +/* private routine to encode sample diameter with range checks */ +static int +encode_diameter(const INTERP2 *ip, double d) +{ + const int ed = ENCODE_DIA(ip, d); + + if (ed <= 0) + return(0); + if (ed >= 0xffff) + return(0xffff); + return(ed); +} + /* Allocate a new set of interpolation samples (caller assigns spt[] array) */ INTERP2 * interp2_alloc(int nsamps) @@ -100,12 +113,131 @@ interp2_spacing(INTERP2 *ip, double mind) ip->dmin = mind; } +/* Compute unnormalized weight for a position relative to a sample */ +double +interp2_wti(INTERP2 *ip, const int i, double x, double y) +{ + double dir, rd, r2, d2; + int ri; + /* get relative direction */ + x -= ip->spt[i][0]; + y -= ip->spt[i][1]; + dir = atan2a(y, x); + dir += 2.*PI*(dir < 0); + /* linear radius interpolation */ + rd = dir * (NI2DIR/2./PI); + ri = (int)rd; + rd -= (double)ri; + rd = (1.-rd)*ip->da[i].dia[ri] + rd*ip->da[i].dia[(ri+1)%NI2DIR]; + rd = ip->smf * DECODE_DIA(ip, rd); + r2 = 2.*rd*rd; + d2 = x*x + y*y; + if (d2 > 21.*r2) /* result would be < 1e-9 */ + return(.0); + /* Gaussian times harmonic weighting */ + return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) ); +} + +/* private call to get grid flag index */ +static int +interp2_flagpos(int fgi[2], INTERP2 *ip, double x, double y) +{ + int inbounds = 0; + + if (ip == NULL) /* paranoia */ + return(-1); + /* need to compute interpolant? */ + if (ip->da == NULL && !interp2_analyze(ip)) + return(-1); + /* get x & y grid positions */ + fgi[0] = (x - ip->smin[0]) * NI2DIM / (ip->smax[0] - ip->smin[0]); + + if (fgi[0] >= NI2DIM) + fgi[0] = NI2DIM-1; + else if (fgi[0] < 0) + fgi[0] = 0; + else + ++inbounds; + + fgi[1] = (y - ip->smin[1]) * NI2DIM / (ip->smax[1] - ip->smin[1]); + + if (fgi[1] >= NI2DIM) + fgi[1] = NI2DIM-1; + else if (fgi[1] < 0) + fgi[1] = 0; + else + ++inbounds; + + return(inbounds == 2); +} + +#define setflg(fl,xi,yi) ((fl)[yi] |= 1<<(xi)) + +#define chkflg(fl,xi,yi) ((fl)[yi]>>(xi) & 1) + +/* private flood function to determine sample influence */ +static void +influence_flood(INTERP2 *ip, const int i, unsigned short visited[NI2DIM], + int xfi, int yfi) +{ + double gx = (xfi+.5)*(1./NI2DIM)*(ip->smax[0] - ip->smin[0]) + + ip->smin[0]; + double gy = (yfi+.5)*(1./NI2DIM)*(ip->smax[1] - ip->smin[1]) + + ip->smin[1]; + double dx = gx - ip->spt[i][0]; + double dy = gy - ip->spt[i][1]; + + setflg(visited, xfi, yfi); + + if (dx*dx + dy*dy > 2.*ip->grid2 && interp2_wti(ip, i, gx, gy) <= 1e-7) + return; + + setflg(ip->da[i].infl, xfi, yfi); + + if (xfi > 0 && !chkflg(visited, xfi-1, yfi)) + influence_flood(ip, i, visited, xfi-1, yfi); + + if (yfi > 0 && !chkflg(visited, xfi, yfi-1)) + influence_flood(ip, i, visited, xfi, yfi-1); + + if (xfi < NI2DIM-1 && !chkflg(visited, xfi+1, yfi)) + influence_flood(ip, i, visited, xfi+1, yfi); + + if (yfi < NI2DIM-1 && !chkflg(visited, xfi, yfi+1)) + influence_flood(ip, i, visited, xfi, yfi+1); +} + +/* private call to compute sample influence maps */ +static void +map_influence(INTERP2 *ip) +{ + unsigned short visited[NI2DIM]; + int fgi[2]; + int i, j; + + for (i = ip->ns; i--; ) { + for (j = NI2DIM; j--; ) { + ip->da[i].infl[j] = 0; + visited[j] = 0; + } + interp2_flagpos(fgi, ip, ip->spt[i][0], ip->spt[i][1]); + + influence_flood(ip, i, visited, fgi[0], fgi[1]); + } +} + /* Modify smoothing parameter by the given factor */ void interp2_smooth(INTERP2 *ip, double sf) { + float old_smf = ip->smf; + if ((ip->smf *= sf) < NI2DSMF) ip->smf = NI2DSMF; + /* need to recompute influence maps? */ + if (ip->da != NULL && (old_smf*.85 > ip->smf) | + (ip->smf > old_smf*1.15)) + map_influence(ip); } /* private call-back to sort position index */ @@ -137,19 +269,6 @@ sort_samples(SAMPORD *sord, const INTERP2 *ip, double qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos); } -/* private routine to encode sample diameter with range checks */ -static int -encode_diameter(const INTERP2 *ip, double d) -{ - const int ed = ENCODE_DIA(ip, d); - - if (ed <= 0) - return(0); - if (ed >= 0xffff) - return(0xffff); - return(ed); -} - /* (Re)compute anisotropic basis function interpolant (normally automatic) */ int interp2_analyze(INTERP2 *ip) @@ -168,33 +287,24 @@ interp2_analyze(INTERP2 *ip) return(0); /* compute sample domain */ ip->smin[0] = ip->smin[1] = FHUGE; - ip->smul[0] = ip->smul[1] = -FHUGE; + ip->smax[0] = ip->smax[1] = -FHUGE; for (i = ip->ns; i--; ) { - if (ip->spt[i][0] < ip->smin[0]) - ip->smin[0] = ip->spt[i][0]; - if (ip->spt[i][0] > ip->smul[0]) - ip->smul[0] = ip->spt[i][0]; - if (ip->spt[i][1] < ip->smin[1]) - ip->smin[1] = ip->spt[i][1]; - if (ip->spt[i][1] > ip->smul[1]) - ip->smul[1] = ip->spt[i][1]; + if (ip->spt[i][0] < ip->smin[0]) ip->smin[0] = ip->spt[i][0]; + if (ip->spt[i][0] > ip->smax[0]) ip->smax[0] = ip->spt[i][0]; + if (ip->spt[i][1] < ip->smin[1]) ip->smin[1] = ip->spt[i][1]; + if (ip->spt[i][1] > ip->smax[1]) ip->smax[1] = ip->spt[i][1]; } - ip->smul[0] -= ip->smin[0]; - ip->smul[1] -= ip->smin[1]; - ip->grid2 = (ip->smul[0]*ip->smul[0] + ip->smul[1]*ip->smul[1]) * - (4./NI2DIM/NI2DIM); + ip->grid2 = ((ip->smax[0]-ip->smin[0])*(ip->smax[0]-ip->smin[0]) + + (ip->smax[1]-ip->smin[1])*(ip->smax[1]-ip->smin[1])) * + (1./NI2DIM/NI2DIM); if (ip->grid2 <= FTINY*ip->dmin*ip->dmin) return(0); - if (ip->smul[0] > FTINY) - ip->smul[0] = NI2DIM / ip->smul[0]; - if (ip->smul[1] > FTINY) - ip->smul[1] = NI2DIM / ip->smul[1]; /* allocate analysis data */ - ip->da = (struct interp2_samp *)calloc( ip->ns, - sizeof(struct interp2_samp) ); + ip->da = (struct interp2_samp *)malloc( + sizeof(struct interp2_samp)*ip->ns ); if (ip->da == NULL) return(0); - /* get temporary arrays */ + /* allocate temporary arrays */ sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns); rightrndx = (int *)malloc(sizeof(int)*ip->ns); leftrndx = (int *)malloc(sizeof(int)*ip->ns); @@ -255,74 +365,37 @@ interp2_analyze(INTERP2 *ip) } } } - free(sortord); /* clean up */ + free(sortord); /* release temp arrays */ free(rightrndx); free(leftrndx); free(endrndx); - return(1); + /* map sample influence areas */ + map_influence(ip); + return(1); /* all done */ } -/* Compute unnormalized weight for a position relative to a sample */ -double -interp2_wti(INTERP2 *ip, const int i, double x, double y) -{ - int xfi, yfi; - double dir, rd, r2, d2; - int ri; - /* need to compute interpolant? */ - if (ip->da == NULL && !interp2_analyze(ip)) - return(0); - /* get grid position */ - xfi = (x - ip->smin[0]) * ip->smul[0]; - if (xfi >= NI2DIM) - xfi = NI2DIM-1; - else - xfi *= (xfi >= 0); - yfi = (y - ip->smin[1]) * ip->smul[1]; - if (yfi >= NI2DIM) - yfi = NI2DIM-1; - else - yfi *= (yfi >= 0); - x -= ip->spt[i][0]; /* check distance */ - y -= ip->spt[i][1]; - d2 = x*x + y*y; - /* zero weight this zone? */ - if (d2 > ip->grid2 && ip->da[i].blkflg[yfi] & 1<da[i].dia[ri] + rd*ip->da[i].dia[(ri+1)%NI2DIR]; - rd = ip->smf * DECODE_DIA(ip, rd); - r2 = 2.*rd*rd; - if (d2 > 21.*r2) { /* result would be < 1e-9 */ - ip->da[i].blkflg[yfi] |= 1<dmin/(ip->dmin + sqrt(d2)) ); -} - /* Assign full set of normalized weights to interpolate the given position */ int interp2_weights(float wtv[], INTERP2 *ip, double x, double y) { double wnorm; + int fgi[2]; int i; - if ((wtv == NULL) | (ip == NULL)) + if (wtv == NULL) return(0); + /* get flag position */ + if (interp2_flagpos(fgi, ip, x, y) < 0) + return(0); wnorm = 0; /* compute raw weights */ - for (i = ip->ns; i--; ) { + for (i = ip->ns; i--; ) + if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) { double wt = interp2_wti(ip, i, x, y); wtv[i] = wt; wnorm += wt; - } + } else + wtv[i] = 0; if (wnorm <= 0) /* too far from all our samples! */ return(0); wnorm = 1./wnorm; /* normalize weights */ @@ -337,16 +410,19 @@ int interp2_topsamp(float wt[], int si[], const int n, INTERP2 *ip, double x, double y) { int nn = 0; + int fgi[2]; double wnorm; int i, j; - if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL)) + if ((n <= 0) | (wt == NULL) | (si == NULL)) return(0); + /* get flag position */ + if (interp2_flagpos(fgi, ip, x, y) < 0) + return(0); /* identify top n weights */ - for (i = ip->ns; i--; ) { + for (i = ip->ns; i--; ) + if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) { const double wti = interp2_wti(ip, i, x, y); - if (wti <= 1e-9) - continue; for (j = nn; j > 0; j--) { if (wt[j-1] >= wti) break; @@ -360,7 +436,7 @@ interp2_topsamp(float wt[], int si[], const int n, INT si[j] = i; nn += (nn < n); } - } + } wnorm = 0; /* normalize sample weights */ for (j = nn; j--; ) wnorm += wt[j];