ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/Development/ray/src/common/interp2d.c

Comparing ray/src/common/interp2d.c (file contents):
Revision 2.9 by greg, Tue Feb 12 18:41:39 2013 UTC vs.
Revision 2.12 by greg, Fri Feb 15 19:15:16 2013 UTC

#	Line 23 | Line 23 | static const char RCSid[] = "$Id$";
23		* to reduce the influence of distant neighbors.  This yields a
24		* smooth interpolation regardless of how the sample points are
25		* initially distributed.  Evaluation is accelerated by use of
26	<	* a fast approximation to the atan2(y,x) function.
26	>	* a fast approximation to the atan2(y,x) function and an array
27	>	* of flags indicating where weights are (nearly) zero.
28		****************************************************************/
29
30		#include <stdio.h>
#	Line 149 | Line 150 | encode_diameter(const INTERP2 *ip, double d)
150		return(ed);
151		}
152
153	+	/* Compute unnormalized weight for a position relative to a sample */
154	+	double
155	+	interp2_wti(INTERP2 *ip, const int i, double x, double y)
156	+	{
157	+	double          dir, rd, r2, d2;
158	+	int             ri;
159	+	/* get relative direction */
160	+	x -= ip->spt[i][0];
161	+	y -= ip->spt[i][1];
162	+	dir = atan2a(y, x);
163	+	dir += 2.*PI*(dir < 0);
164	+	/* linear radius interpolation */
165	+	rd = dir * (NI2DIR/2./PI);
166	+	ri = (int)rd;
167	+	rd -= (double)ri;
168	+	rd = (1.-rd)*ip->da[i].dia[ri] + rd*ip->da[i].dia[(ri+1)%NI2DIR];
169	+	rd = ip->smf * DECODE_DIA(ip, rd);
170	+	r2 = 2.*rd*rd;
171	+	d2 = x*x + y*y;
172	+	if (d2 > 21.*r2)        /* result would be < 1e-9 */
173	+	return(.0);
174	+	/* Gaussian times harmonic weighting */
175	+	return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
176	+	}
177	+
178	+	/* private call to get grid flag index */
179	+	static int
180	+	interp2_flagpos(int fgi[2], INTERP2 *ip, double x, double y)
181	+	{
182	+	int     inbounds = 0;
183	+
184	+	if (ip == NULL)         /* paranoia */
185	+	return(-1);
186	+	/* need to compute interpolant? */
187	+	if (ip->da == NULL && !interp2_analyze(ip))
188	+	return(-1);
189	+	/* get x & y grid positions */
190	+	fgi[0] = (x - ip->smin[0]) * NI2DIM / (ip->smax[0] - ip->smin[0]);
191	+
192	+	if (fgi[0] >= NI2DIM)
193	+	fgi[0] = NI2DIM-1;
194	+	else if (fgi[0] < 0)
195	+	fgi[0] = 0;
196	+	else
197	+	++inbounds;
198	+
199	+	fgi[1] = (y - ip->smin[1]) * NI2DIM / (ip->smax[1] - ip->smin[1]);
200	+
201	+	if (fgi[1] >= NI2DIM)
202	+	fgi[1] = NI2DIM-1;
203	+	else if (fgi[1] < 0)
204	+	fgi[1] = 0;
205	+	else
206	+	++inbounds;
207	+
208	+	return(inbounds == 2);
209	+	}
210	+
211	+	#define setflg(fl,xi,yi)        ((fl)[yi] |= 1<<(xi))
212	+
213	+	#define chkflg(fl,xi,yi)        ((fl)[yi]>>(xi) & 1)
214	+
215	+	/* private flood function to determine sample influence */
216	+	static void
217	+	influence_flood(INTERP2 *ip, const int i, unsigned short visited[NI2DIM],
218	+	int xfi, int yfi)
219	+	{
220	+	double  gx = (xfi+.5)*(1./NI2DIM)*(ip->smax[0] - ip->smin[0]) +
221	+	ip->smin[0];
222	+	double  gy = (yfi+.5)*(1./NI2DIM)*(ip->smax[1] - ip->smin[1]) +
223	+	ip->smin[1];
224	+	double  dx = gx - ip->spt[i][0];
225	+	double  dy = gy - ip->spt[i][1];
226	+
227	+	setflg(visited, xfi, yfi);
228	+
229	+	if (dx*dx + dy*dy > 2.*ip->grid2 && interp2_wti(ip, i, gx, gy) <= 1e-7)
230	+	return;
231	+
232	+	setflg(ip->da[i].infl, xfi, yfi);
233	+
234	+	if (xfi > 0 && !chkflg(visited, xfi-1, yfi))
235	+	influence_flood(ip, i, visited, xfi-1, yfi);
236	+
237	+	if (yfi > 0 && !chkflg(visited, xfi, yfi-1))
238	+	influence_flood(ip, i, visited, xfi, yfi-1);
239	+
240	+	if (xfi < NI2DIM-1 && !chkflg(visited, xfi+1, yfi))
241	+	influence_flood(ip, i, visited, xfi+1, yfi);
242	+
243	+	if (yfi < NI2DIM-1 && !chkflg(visited, xfi, yfi+1))
244	+	influence_flood(ip, i, visited, xfi, yfi+1);
245	+	}
246	+
247		/* (Re)compute anisotropic basis function interpolant (normally automatic) */
248		int
249		interp2_analyze(INTERP2 *ip)
250		{
251		SAMPORD *sortord;
252		int     *rightrndx, *leftrndx, *endrndx;
253	<	int     bd;
253	>	int     i, j, bd;
254		/* sanity checks */
255	<	if (ip == NULL || (ip->ns <= 1) | (ip->dmin <= 0))
255	>	if (ip == NULL)
256		return(0);
257	<	/* need to allocate? */
258	<	if (ip->da == NULL) {
259	<	ip->da = (unsigned short (*)[NI2DIR])malloc(
165	<	sizeof(unsigned short)*NI2DIR*ip->ns);
166	<	if (ip->da == NULL)
167	<	return(0);
257	>	if (ip->da != NULL) {           /* free previous data if any */
258	>	free(ip->da);
259	>	ip->da = NULL;
260		}
261	<	/* get temporary arrays */
261	>	if ((ip->ns <= 1) | (ip->dmin <= 0))
262	>	return(0);
263	>	/* compute sample domain */
264	>	ip->smin[0] = ip->smin[1] = FHUGE;
265	>	ip->smax[0] = ip->smax[1] = -FHUGE;
266	>	for (i = ip->ns; i--; ) {
267	>	if (ip->spt[i][0] < ip->smin[0]) ip->smin[0] = ip->spt[i][0];
268	>	if (ip->spt[i][0] > ip->smax[0]) ip->smax[0] = ip->spt[i][0];
269	>	if (ip->spt[i][1] < ip->smin[1]) ip->smin[1] = ip->spt[i][1];
270	>	if (ip->spt[i][1] > ip->smax[1]) ip->smax[1] = ip->spt[i][1];
271	>	}
272	>	ip->grid2 = ((ip->smax[0]-ip->smin[0])*(ip->smax[0]-ip->smin[0]) +
273	>	(ip->smax[1]-ip->smin[1])*(ip->smax[1]-ip->smin[1])) *
274	>	(1./NI2DIM/NI2DIM);
275	>	if (ip->grid2 <= FTINY*ip->dmin*ip->dmin)
276	>	return(0);
277	>	/* allocate analysis data */
278	>	ip->da = (struct interp2_samp *)calloc( ip->ns,
279	>	sizeof(struct interp2_samp) );
280	>	if (ip->da == NULL)
281	>	return(0);
282	>	/* allocate temporary arrays */
283		sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
284		rightrndx = (int *)malloc(sizeof(int)*ip->ns);
285		leftrndx = (int *)malloc(sizeof(int)*ip->ns);
#	Line 178 | Line 291 | interp2_analyze(INTERP2 *ip)
291		for (bd = 0; bd < NI2DIR/2; bd++) {
292		const double        ang = 2.*PI/NI2DIR*bd;
293		int                 *sptr;
181	–	int                 i;
294		/* create right reverse index */
295		if (bd) {                   /* re-use from previous iteration? */
296		sptr = rightrndx;
#	Line 207 | Line 319 | interp2_analyze(INTERP2 *ip)
319		/* find nearest neighbors each side */
320		for (i = ip->ns; i--; ) {
321		const int       ii = sortord[i].si;
210	–	int             j;
322		/* preload with large radii */
323	<	ip->da[ii][bd] = ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
324	<	.5*(sortord[ip->ns-1].dm - sortord[0].dm));
323	>	ip->da[ii].dia[bd] =
324	>	ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
325	>	.5*(sortord[ip->ns-1].dm - sortord[0].dm));
326		for (j = i; ++j < ip->ns; )     /* nearest above */
327		if (rightrndx[sortord[j].si] > rightrndx[ii] &&
328		leftrndx[sortord[j].si] < leftrndx[ii]) {
329	<	ip->da[ii][bd] = encode_diameter(ip,
329	>	ip->da[ii].dia[bd] = encode_diameter(ip,
330		sortord[j].dm - sortord[i].dm);
331		break;
332		}
333		for (j = i; j-- > 0; )          /* nearest below */
334		if (rightrndx[sortord[j].si] < rightrndx[ii] &&
335		leftrndx[sortord[j].si] > leftrndx[ii]) {
336	<	ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
336	>	ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
337		sortord[i].dm - sortord[j].dm);
338		break;
339		}
340		}
341		}
342	<	free(sortord);                  /* clean up */
342	>	free(sortord);                  /* release temp arrays */
343		free(rightrndx);
344		free(leftrndx);
345		free(endrndx);
346	<	return(1);
347	<	}
346	>	/* fill influence maps */
347	>	for (i = ip->ns; i--; ) {
348	>	unsigned short  visited[NI2DIM];
349	>	int             fgi[2];
350
351	<	/* private call returns raw weight for a particular sample */
352	<	static double
353	<	get_wt(const INTERP2 *ip, const int i, double x, double y)
354	<	{
355	<	double  dir, rd, r2, d2;
242	<	int     ri;
243	<	/* get relative direction */
244	<	x -= ip->spt[i][0];
245	<	y -= ip->spt[i][1];
246	<	dir = atan2a(y, x);
247	<	dir += 2.*PI*(dir < 0);
248	<	/* linear radius interpolation */
249	<	rd = dir * (NI2DIR/2./PI);
250	<	ri = (int)rd;
251	<	rd -= (double)ri;
252	<	rd = (1.-rd)*ip->da[i][ri] + rd*ip->da[i][(ri+1)%NI2DIR];
253	<	rd = ip->smf * DECODE_DIA(ip, rd);
254	<	r2 = 2.*rd*rd;
255	<	d2 = x*x + y*y;
256	<	if (d2 > 21.*r2)        /* result would be < 1e-9 */
257	<	return(.0);
258	<	/* Gaussian times harmonic weighting */
259	<	return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
351	>	for (j = NI2DIM; j--; ) visited[j] = 0;
352	>	interp2_flagpos(fgi, ip, ip->spt[i][0], ip->spt[i][1]);
353	>	influence_flood(ip, i, visited, fgi[0], fgi[1]);
354	>	}
355	>	return(1);                      /* all done */
356		}
357
358		/* Assign full set of normalized weights to interpolate the given position */
#	Line 264 | Line 360 | int
360		interp2_weights(float wtv[], INTERP2 *ip, double x, double y)
361		{
362		double  wnorm;
363	+	int     fgi[2];
364		int     i;
365
366	<	if ((wtv == NULL) | (ip == NULL))
366	>	if (wtv == NULL)
367		return(0);
368	<	/* need to compute interpolant? */
369	<	if (ip->da == NULL && !interp2_analyze(ip))
368	>	/* get flag position */
369	>	if (interp2_flagpos(fgi, ip, x, y) < 0)
370		return(0);
371
372		wnorm = 0;                      /* compute raw weights */
373	<	for (i = ip->ns; i--; ) {
374	<	double  wt = get_wt(ip, i, x, y);
373	>	for (i = ip->ns; i--; )
374	>	if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
375	>	double  wt = interp2_wti(ip, i, x, y);
376		wtv[i] = wt;
377		wnorm += wt;
378	<	}
378	>	} else
379	>	wtv[i] = 0;
380		if (wnorm <= 0)                 /* too far from all our samples! */
381		return(0);
382		wnorm = 1./wnorm;               /* normalize weights */
#	Line 292 | Line 391 | int
391		interp2_topsamp(float wt[], int si[], const int n, INTERP2 *ip, double x, double y)
392		{
393		int     nn = 0;
394	+	int     fgi[2];
395		double  wnorm;
396		int     i, j;
397
398	<	if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
398	>	if ((n <= 0) | (wt == NULL) | (si == NULL))
399		return(0);
400	<	/* need to compute interpolant? */
401	<	if (ip->da == NULL && !interp2_analyze(ip))
400	>	/* get flag position */
401	>	if (interp2_flagpos(fgi, ip, x, y) < 0)
402		return(0);
403		/* identify top n weights */
404	<	for (i = ip->ns; i--; ) {
405	<	const double    wti = get_wt(ip, i, x, y);
406	<	if (wti <= 1e-9)
307	<	continue;
404	>	for (i = ip->ns; i--; )
405	>	if (chkflg(ip->da[i].infl, fgi[0], fgi[1])) {
406	>	const double    wti = interp2_wti(ip, i, x, y);
407		for (j = nn; j > 0; j--) {
408		if (wt[j-1] >= wti)
409		break;
#	Line 318 | Line 417 | interp2_topsamp(float wt[], int si[], const int n, INT
417		si[j] = i;
418		nn += (nn < n);
419		}
420	<	}
420	>	}
421		wnorm = 0;                      /* normalize sample weights */
422		for (j = nn; j--; )
423		wnorm += wt[j];

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)