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.3 by greg, Mon Feb 11 20:01:15 2013 UTC vs.
Revision 2.10 by greg, Thu Feb 14 19:57:10 2013 UTC

#	Line 9 | Line 9 | static const char RCSid[] = "$Id$";
9
10		#include "copyright.h"
11
12	<	/*************************************************************
12	>	/***************************************************************
13		* This is a general method for 2-D interpolation similar to
14		* radial basis functions but allowing for a good deal of local
15		* anisotropy in the point distribution.  Each sample point
16		* is examined to determine the closest neighboring samples in
17		* each of NI2DIR surrounding directions.  To speed this
18	<	* calculation, we sort the data into 3 half-planes and
19	<	* perform simple tests to see which neighbor is closest in
20	<	* a each direction.  Once we have our approximate neighborhood
18	>	* calculation, we sort the data into half-planes and apply
19	>	* simple tests to see which neighbor is closest in each
20	>	* angular slice.  Once we have our approximate neighborhood
21		* for a sample, we can use it in a modified Gaussian weighting
22	<	* scheme with anisotropic surround.  Harmonic weighting is added
22	>	* with allowing local anisotropy.  Harmonic weighting is added
23		* to reduce the influence of distant neighbors.  This yields a
24		* smooth interpolation regardless of how the sample points are
25	<	* initiallydistributed.  Evaluation is accelerated by use of a
26	<	* fast approximation to the atan2(y,x) function.
27	<	**************************************************************/
25	>	* initially distributed.  Evaluation is accelerated by use of
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>
31		#include <stdlib.h>
32		#include "rtmath.h"
33		#include "interp2d.h"
34
35	<	#define DECODE_RAD(ip,er)       ((ip)->rmin*(1. + .5*(er)))
36	<	#define ENCODE_RAD(ip,r)        ((int)(2.*(r)/(ip)->rmin) - 2)
35	>	#define DECODE_DIA(ip,ed)       ((ip)->dmin*(1. + .5*(ed)))
36	>	#define ENCODE_DIA(ip,d)        ((int)(2.*(d)/(ip)->dmin) - 2)
37
38		/* Sample order (private) */
39		typedef struct {
#	Line 54 | Line 55 | interp2_alloc(int nsamps)
55		return(NULL);
56
57		nip->ns = nsamps;
58	<	nip->rmin = .5;         /* default radius minimum */
58	>	nip->dmin = 1;          /* default minimum diameter */
59		nip->smf = NI2DSMF;     /* default smoothing factor */
60	<	nip->ra = NULL;
60	>	nip->da = NULL;
61		/* caller must assign spt[] array */
62		return(nip);
63		}
#	Line 71 | Line 72 | interp2_realloc(INTERP2 *ip, int nsamps)
72		interp2_free(ip);
73		return(NULL);
74		}
75	<	if (nsamps == ip->ns);
75	>	if (nsamps == ip->ns)
76		return(ip);
77	<	if (ip->ra != NULL) {   /* will need to recompute distribution */
78	<	free(ip->ra);
79	<	ip->ra = NULL;
77	>	if (ip->da != NULL) {   /* will need to recompute distribution */
78	>	free(ip->da);
79	>	ip->da = NULL;
80		}
81		ip = (INTERP2 *)realloc(ip, sizeof(INTERP2)+sizeof(float)*2*(nsamps-1));
82		if (ip == NULL)
#	Line 84 | Line 85 | interp2_realloc(INTERP2 *ip, int nsamps)
85		return(ip);
86		}
87
88	+	/* Set minimum distance under which samples will start to merge */
89	+	void
90	+	interp2_spacing(INTERP2 *ip, double mind)
91	+	{
92	+	if (mind <= 0)
93	+	return;
94	+	if ((.998*ip->dmin <= mind) & (mind <= 1.002*ip->dmin))
95	+	return;
96	+	if (ip->da != NULL) {   /* will need to recompute distribution */
97	+	free(ip->da);
98	+	ip->da = NULL;
99	+	}
100	+	ip->dmin = mind;
101	+	}
102	+
103	+	/* Modify smoothing parameter by the given factor */
104	+	void
105	+	interp2_smooth(INTERP2 *ip, double sf)
106	+	{
107	+	if ((ip->smf *= sf) < NI2DSMF)
108	+	ip->smf = NI2DSMF;
109	+	}
110	+
111		/* private call-back to sort position index */
112		static int
113		cmp_spos(const void *p1, const void *p2)
#	Line 113 | Line 137 | sort_samples(SAMPORD *sord, const INTERP2 *ip, double
137		qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
138		}
139
140	<	/* private routine to encode radius with range checks */
140	>	/* private routine to encode sample diameter with range checks */
141		static int
142	<	encode_radius(const INTERP2 *ip, double r)
142	>	encode_diameter(const INTERP2 *ip, double d)
143		{
144	<	const int       er = ENCODE_RAD(ip, r);
144	>	const int       ed = ENCODE_DIA(ip, d);
145
146	<	if (er <= 0)
146	>	if (ed <= 0)
147		return(0);
148	<	if (er >= 0xffff)
148	>	if (ed >= 0xffff)
149		return(0xffff);
150	<	return(er);
150	>	return(ed);
151		}
152
153		/* (Re)compute anisotropic basis function interpolant (normally automatic) */
#	Line 132 | Line 156 | interp2_analyze(INTERP2 *ip)
156		{
157		SAMPORD *sortord;
158		int     *rightrndx, *leftrndx, *endrndx;
159	<	int     bd;
159	>	int     i, bd;
160		/* sanity checks */
161	<	if (ip == NULL || (ip->ns <= 1) | (ip->rmin <= 0))
161	>	if (ip == NULL)
162		return(0);
163	<	/* need to allocate? */
164	<	if (ip->ra == NULL) {
165	<	ip->ra = (unsigned short (*)[NI2DIR])malloc(
142	<	sizeof(unsigned short)*NI2DIR*ip->ns);
143	<	if (ip->ra == NULL)
144	<	return(0);
163	>	if (ip->da != NULL) {           /* free previous data if any */
164	>	free(ip->da);
165	>	ip->da = NULL;
166		}
167	+	if ((ip->ns <= 1) | (ip->dmin <= 0))
168	+	return(0);
169	+	/* compute sample domain */
170	+	ip->smin[0] = ip->smin[1] = FHUGE;
171	+	ip->smul[0] = ip->smul[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];
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];
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);
186	+	if (ip->grid2 <= FTINY*ip->dmin*ip->dmin)
187	+	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];
192	+	/* allocate analysis data */
193	+	ip->da = (struct interp2_samp *)calloc( ip->ns,
194	+	sizeof(struct interp2_samp) );
195	+	if (ip->da == NULL)
196	+	return(0);
197		/* get temporary arrays */
198		sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
199		rightrndx = (int *)malloc(sizeof(int)*ip->ns);
#	Line 155 | Line 206 | interp2_analyze(INTERP2 *ip)
206		for (bd = 0; bd < NI2DIR/2; bd++) {
207		const double        ang = 2.*PI/NI2DIR*bd;
208		int                 *sptr;
158	–	int                 i;
209		/* create right reverse index */
210		if (bd) {                   /* re-use from previous iteration? */
211		sptr = rightrndx;
#	Line 186 | Line 236 | interp2_analyze(INTERP2 *ip)
236		const int       ii = sortord[i].si;
237		int             j;
238		/* preload with large radii */
239	<	ip->ra[ii][bd] = ip->ra[ii][bd+NI2DIR/2] = encode_radius(ip,
240	<	.25*(sortord[ip->ns-1].dm - sortord[0].dm));
239	>	ip->da[ii].dia[bd] =
240	>	ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
241	>	.5*(sortord[ip->ns-1].dm - sortord[0].dm));
242		for (j = i; ++j < ip->ns; )     /* nearest above */
243		if (rightrndx[sortord[j].si] > rightrndx[ii] &&
244		leftrndx[sortord[j].si] < leftrndx[ii]) {
245	<	ip->ra[ii][bd] = encode_radius(ip,
246	<	.5*(sortord[j].dm - sortord[i].dm));
245	>	ip->da[ii].dia[bd] = encode_diameter(ip,
246	>	sortord[j].dm - sortord[i].dm);
247		break;
248		}
249		for (j = i; j-- > 0; )          /* nearest below */
250		if (rightrndx[sortord[j].si] < rightrndx[ii] &&
251		leftrndx[sortord[j].si] > leftrndx[ii]) {
252	<	ip->ra[ii][bd+NI2DIR/2] = encode_radius(ip,
253	<	.5*(sortord[i].dm - sortord[j].dm));
252	>	ip->da[ii].dia[bd+NI2DIR/2] = encode_diameter(ip,
253	>	sortord[i].dm - sortord[j].dm);
254		break;
255		}
256		}
#	Line 211 | Line 262 | interp2_analyze(INTERP2 *ip)
262		return(1);
263		}
264
265	<	/* private call returns raw weight for a particular sample */
266	<	static double
267	<	get_wt(const INTERP2 *ip, const int i, double x, double y)
265	>	/* Compute unnormalized weight for a position relative to a sample */
266	>	double
267	>	interp2_wti(INTERP2 *ip, const int i, double x, double y)
268		{
269	<	double  dir, rd, d2;
269	>	int     xfi, yfi;
270	>	double  dir, rd, r2, d2;
271		int     ri;
272	<	/* get relative direction */
273	<	x -= ip->spt[i][0];
272	>	/* 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 */
287		y -= ip->spt[i][1];
288	<	dir = atan2a(y, x);
288	>	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 */
294		dir += 2.*PI*(dir < 0);
295		/* linear radius interpolation */
296		rd = dir * (NI2DIR/2./PI);
297		ri = (int)rd;
298		rd -= (double)ri;
299	<	rd = (1.-rd)*ip->ra[i][ri] + rd*ip->ra[i][(ri+1)%NI2DIR];
300	<	rd = ip->smf * DECODE_RAD(ip, rd);
301	<	d2 = x*x + y*y;
299	>	rd = (1.-rd)*ip->da[i].dia[ri] + rd*ip->da[i].dia[(ri+1)%NI2DIR];
300	>	rd = ip->smf * DECODE_DIA(ip, rd);
301	>	r2 = 2.*rd*rd;
302	>	if (d2 > 21.*r2) {      /* result would be < 1e-9 */
303	>	ip->da[i].blkflg[yfi] |= 1<<xfi;
304	>	return(.0);
305	>	}
306		/* Gaussian times harmonic weighting */
307	<	return( exp(d2/(-2.*rd*rd)) * ip->rmin/(ip->rmin + sqrt(d2)) );
307	>	return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
308		}
309
310		/* Assign full set of normalized weights to interpolate the given position */
#	Line 242 | Line 316 | interp2_weights(float wtv[], INTERP2 *ip, double x, do
316
317		if ((wtv == NULL) | (ip == NULL))
318		return(0);
245	–	/* need to compute interpolant? */
246	–	if (ip->ra == NULL && !interp2_analyze(ip))
247	–	return(0);
319
320		wnorm = 0;                      /* compute raw weights */
321		for (i = ip->ns; i--; ) {
322	<	double  wt = get_wt(ip, i, x, y);
322	>	double  wt = interp2_wti(ip, i, x, y);
323		wtv[i] = wt;
324		wnorm += wt;
325		}
#	Line 271 | Line 342 | interp2_topsamp(float wt[], int si[], const int n, INT
342
343		if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
344		return(0);
274	–	/* need to compute interpolant? */
275	–	if (ip->ra == NULL && !interp2_analyze(ip))
276	–	return(0);
345		/* identify top n weights */
346		for (i = ip->ns; i--; ) {
347	<	const double    wti = get_wt(ip, i, x, y);
347	>	const double    wti = interp2_wti(ip, i, x, y);
348	>	if (wti <= 1e-9)
349	>	continue;
350		for (j = nn; j > 0; j--) {
351		if (wt[j-1] >= wti)
352		break;
#	Line 308 | Line 378 | interp2_free(INTERP2 *ip)
378		{
379		if (ip == NULL)
380		return;
381	<	if (ip->ra != NULL)
382	<	free(ip->ra);
381	>	if (ip->da != NULL)
382	>	free(ip->da);
383		free(ip);
384		}

Diff Legend

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