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.10 by greg, Thu Feb 14 19:57:10 2013 UTC vs.
Revision 2.11 by greg, Fri Feb 15 01:26:47 2013 UTC

#	Line 168 | Line 168 | interp2_analyze(INTERP2 *ip)
168		return(0);
169		/* compute sample domain */
170		ip->smin[0] = ip->smin[1] = FHUGE;
171	<	ip->smul[0] = ip->smul[1] = -FHUGE;
171	>	ip->smax[0] = ip->smax[1] = -FHUGE;
172		for (i = ip->ns; i--; ) {
173		if (ip->spt[i][0] < ip->smin[0])
174		ip->smin[0] = ip->spt[i][0];
175	<	if (ip->spt[i][0] > ip->smul[0])
176	<	ip->smul[0] = ip->spt[i][0];
175	>	if (ip->spt[i][0] > ip->smax[0])
176	>	ip->smax[0] = ip->spt[i][0];
177		if (ip->spt[i][1] < ip->smin[1])
178		ip->smin[1] = ip->spt[i][1];
179	<	if (ip->spt[i][1] > ip->smul[1])
180	<	ip->smul[1] = ip->spt[i][1];
179	>	if (ip->spt[i][1] > ip->smax[1])
180	>	ip->smax[1] = ip->spt[i][1];
181		}
182	<	ip->smul[0] -= ip->smin[0];
183	<	ip->smul[1] -= ip->smin[1];
184	<	ip->grid2 = (ip->smul[0]*ip->smul[0] + ip->smul[1]*ip->smul[1]) *
185	<	(4./NI2DIM/NI2DIM);
182	>	ip->grid2 = ((ip->smax[0]-ip->smin[0])*(ip->smax[0]-ip->smin[0]) +
183	>	(ip->smax[1]-ip->smin[1])*(ip->smax[1]-ip->smin[1])) *
184	>	(1./NI2DIM/NI2DIM);
185		if (ip->grid2 <= FTINY*ip->dmin*ip->dmin)
186		return(0);
188	–	if (ip->smul[0] > FTINY)
189	–	ip->smul[0] = NI2DIM / ip->smul[0];
190	–	if (ip->smul[1] > FTINY)
191	–	ip->smul[1] = NI2DIM / ip->smul[1];
187		/* allocate analysis data */
188		ip->da = (struct interp2_samp *)calloc( ip->ns,
189		sizeof(struct interp2_samp) );
#	Line 266 | Line 261 | interp2_analyze(INTERP2 *ip)
261		double
262		interp2_wti(INTERP2 *ip, const int i, double x, double y)
263		{
264	<	int     xfi, yfi;
265	<	double  dir, rd, r2, d2;
266	<	int     ri;
267	<	/* need to compute interpolant? */
273	<	if (ip->da == NULL && !interp2_analyze(ip))
274	<	return(0);
275	<	/* get grid position */
276	<	xfi = (x - ip->smin[0]) * ip->smul[0];
277	<	if (xfi >= NI2DIM)
278	<	xfi = NI2DIM-1;
279	<	else
280	<	xfi *= (xfi >= 0);
281	<	yfi = (y - ip->smin[1]) * ip->smul[1];
282	<	if (yfi >= NI2DIM)
283	<	yfi = NI2DIM-1;
284	<	else
285	<	yfi *= (yfi >= 0);
286	<	x -= ip->spt[i][0];     /* check distance */
264	>	double          dir, rd, r2, d2;
265	>	int             ri;
266	>	/* get relative direction */
267	>	x -= ip->spt[i][0];
268		y -= ip->spt[i][1];
269	<	d2 = x*x + y*y;
289	<	/* zero weight this zone? */
290	<	if (d2 > ip->grid2 && ip->da[i].blkflg[yfi] & 1<<xfi)
291	<	return(.0);
292	<
293	<	dir = atan2a(y, x);     /* get relative direction */
269	>	dir = atan2a(y, x);
270		dir += 2.*PI*(dir < 0);
271		/* linear radius interpolation */
272		rd = dir * (NI2DIR/2./PI);
#	Line 299 | Line 275 | interp2_wti(INTERP2 *ip, const int i, double x, double
275		rd = (1.-rd)*ip->da[i].dia[ri] + rd*ip->da[i].dia[(ri+1)%NI2DIR];
276		rd = ip->smf * DECODE_DIA(ip, rd);
277		r2 = 2.*rd*rd;
278	<	if (d2 > 21.*r2) {      /* result would be < 1e-9 */
279	<	ip->da[i].blkflg[yfi] |= 1<<xfi;
278	>	d2 = x*x + y*y;
279	>	if (d2 > 21.*r2)        /* result would be < 1e-9 */
280		return(.0);
305	–	}
281		/* Gaussian times harmonic weighting */
282		return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
283		}
284
285	+	/* private call to get grid flag index */
286	+	static int
287	+	interp2_flagpos(int fgi[2], INTERP2 *ip, double x, double y)
288	+	{
289	+	int     ingrid = 1;
290	+
291	+	if (ip == NULL)         /* paranoia */
292	+	return(-1);
293	+	/* need to compute interpolant? */
294	+	if (ip->da == NULL && !interp2_analyze(ip))
295	+	return(-1);
296	+	/* get grid position */
297	+	fgi[0] = (x - ip->smin[0]) * NI2DIM / (ip->smax[0] - ip->smin[0]);
298	+	if (fgi[0] >= NI2DIM) {
299	+	fgi[0] = NI2DIM-1;
300	+	ingrid = 0;
301	+	} else if (fgi[0] < 0) {
302	+	fgi[0] = 0;
303	+	ingrid = 0;
304	+	}
305	+	fgi[1] = (y - ip->smin[1]) * NI2DIM / (ip->smax[1] - ip->smin[1]);
306	+	if (fgi[1] >= NI2DIM) {
307	+	fgi[1] = NI2DIM-1;
308	+	ingrid = 0;
309	+	} else if (fgi[1] < 0) {
310	+	fgi[1] = 0;
311	+	ingrid = 0;
312	+	}
313	+	return(ingrid);
314	+	}
315	+
316	+	/* private call to set black flag if not too close to the given sample */
317	+	static void
318	+	setblk(INTERP2 *ip, const int i, const int gi[2])
319	+	{
320	+	double  dx = (gi[0]+.5)*(1./NI2DIM)*(ip->smax[0] - ip->smin[0]) +
321	+	ip->smin[0] - ip->spt[i][0];
322	+	double  dy = (gi[1]+.5)*(1./NI2DIM)*(ip->smax[1] - ip->smin[1]) +
323	+	ip->smin[1] - ip->spt[i][1];
324	+
325	+	if (dx*dx + dy*dy > 2.*ip->grid2)
326	+	ip->da[i].blkflg[gi[1]] |= 1<<gi[0];
327	+	}
328	+
329	+	#define chkblk(ip,i,gi) ((ip)->da[i].blkflg[(gi)[1]]>>(gi)[0] & 1)
330	+
331		/* Assign full set of normalized weights to interpolate the given position */
332		int
333		interp2_weights(float wtv[], INTERP2 *ip, double x, double y)
334		{
335		double  wnorm;
336	+	int     fgi[2];
337	+	int     ingrid;
338		int     i;
339
340	<	if ((wtv == NULL) | (ip == NULL))
340	>	if (wtv == NULL)
341		return(0);
342	+	/* get flag position */
343	+	if ((ingrid = interp2_flagpos(fgi, ip, x, y)) < 0)
344	+	return(0);
345
346		wnorm = 0;                      /* compute raw weights */
347	<	for (i = ip->ns; i--; ) {
347	>	for (i = ip->ns; i--; )
348	>	if (chkblk(ip, i, fgi)) {
349	>	wtv[i] = 0;
350	>	} else {
351		double  wt = interp2_wti(ip, i, x, y);
352		wtv[i] = wt;
353		wnorm += wt;
354	<	}
354	>	if (wt <= 1e-9 && ingrid)
355	>	setblk(ip, i, fgi);
356	>	}
357		if (wnorm <= 0)                 /* too far from all our samples! */
358		return(0);
359		wnorm = 1./wnorm;               /* normalize weights */
#	Line 337 | Line 368 | int
368		interp2_topsamp(float wt[], int si[], const int n, INTERP2 *ip, double x, double y)
369		{
370		int     nn = 0;
371	+	int     fgi[2];
372	+	int     ingrid;
373		double  wnorm;
374		int     i, j;
375
376	<	if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
376	>	if ((n <= 0) | (wt == NULL) | (si == NULL))
377		return(0);
378	+	/* get flag position */
379	+	if ((ingrid = interp2_flagpos(fgi, ip, x, y)) < 0)
380	+	return(0);
381		/* identify top n weights */
382	<	for (i = ip->ns; i--; ) {
382	>	for (i = ip->ns; i--; )
383	>	if (!chkblk(ip, i, fgi)) {
384		const double    wti = interp2_wti(ip, i, x, y);
385	<	if (wti <= 1e-9)
385	>	if (wti <= 1e-9) {
386	>	if (ingrid)
387	>	setblk(ip, i, fgi);
388		continue;
389	+	}
390		for (j = nn; j > 0; j--) {
391		if (wt[j-1] >= wti)
392		break;
#	Line 360 | Line 400 | interp2_topsamp(float wt[], int si[], const int n, INT
400		si[j] = i;
401		nn += (nn < n);
402		}
403	<	}
403	>	}
404		wnorm = 0;                      /* normalize sample weights */
405		for (j = nn; j--; )
406		wnorm += wt[j];

Diff Legend

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