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.1 by greg, Sat Feb 9 00:55:40 2013 UTC vs.
Revision 2.3 by greg, Mon Feb 11 20:01:15 2013 UTC

#	Line 18 | Line 18 | static const char RCSid[] = "$Id$";
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
21	<	* for a sample, we can use it in a Gaussian weighting scheme
22	<	* with anisotropic surround.  This gives us a fairly smooth
23	<	* interpolation however the sample points may be initially
24	<	* distributed.  Evaluation is accelerated by use of a fast
25	<	* approximation to the atan2(y,x) function.
21	>	* for a sample, we can use it in a modified Gaussian weighting
22	>	* scheme with anisotropic surround.  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		**************************************************************/
28
29		#include <stdio.h>
#	Line 39 | Line 40 | typedef struct {
40		float   dm;             /* distance measure in this direction */
41		} SAMPORD;
42
43	<	/* Allocate a new set of interpolation samples */
43	>	/* Allocate a new set of interpolation samples (caller assigns spt[] array) */
44		INTERP2 *
45		interp2_alloc(int nsamps)
46		{
#	Line 60 | Line 61 | interp2_alloc(int nsamps)
61		return(nip);
62		}
63
64	+	/* Resize interpolation array (caller must assign any new values) */
65	+	INTERP2 *
66	+	interp2_realloc(INTERP2 *ip, int nsamps)
67	+	{
68	+	if (ip == NULL)
69	+	return(interp2_alloc(nsamps));
70	+	if (nsamps <= 1) {
71	+	interp2_free(ip);
72	+	return(NULL);
73	+	}
74	+	if (nsamps == ip->ns);
75	+	return(ip);
76	+	if (ip->ra != NULL) {   /* will need to recompute distribution */
77	+	free(ip->ra);
78	+	ip->ra = NULL;
79	+	}
80	+	ip = (INTERP2 *)realloc(ip, sizeof(INTERP2)+sizeof(float)*2*(nsamps-1));
81	+	if (ip == NULL)
82	+	return(NULL);
83	+	ip->ns = nsamps;
84	+	return(ip);
85	+	}
86	+
87		/* private call-back to sort position index */
88		static int
89		cmp_spos(const void *p1, const void *p2)
#	Line 74 | Line 98 | cmp_spos(const void *p1, const void *p2)
98		return 0;
99		}
100
101	+	/* private routine to order samples in a particular direction */
102	+	static void
103	+	sort_samples(SAMPORD *sord, const INTERP2 *ip, double ang)
104	+	{
105	+	const double    cosd = cos(ang);
106	+	const double    sind = sin(ang);
107	+	int             i;
108	+
109	+	for (i = ip->ns; i--; ) {
110	+	sord[i].si = i;
111	+	sord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
112	+	}
113	+	qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
114	+	}
115	+
116		/* private routine to encode radius with range checks */
117		static int
118		encode_radius(const INTERP2 *ip, double r)
#	Line 87 | Line 126 | encode_radius(const INTERP2 *ip, double r)
126		return(er);
127		}
128
129	<	/* Compute anisotropic Gaussian basis function interpolant */
130	<	static int
131	<	interp2_compute(INTERP2 *ip)
129	>	/* (Re)compute anisotropic basis function interpolant (normally automatic) */
130	>	int
131	>	interp2_analyze(INTERP2 *ip)
132		{
133		SAMPORD *sortord;
134	<	int     *rightrndx, *leftrndx;
134	>	int     *rightrndx, *leftrndx, *endrndx;
135		int     bd;
136		/* sanity checks */
137		if (ip == NULL || (ip->ns <= 1) | (ip->rmin <= 0))
#	Line 108 | Line 147 | interp2_compute(INTERP2 *ip)
147		sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
148		rightrndx = (int *)malloc(sizeof(int)*ip->ns);
149		leftrndx = (int *)malloc(sizeof(int)*ip->ns);
150	<	if ((sortord == NULL) | (rightrndx == NULL) | (leftrndx == NULL))
150	>	endrndx = (int *)malloc(sizeof(int)*ip->ns);
151	>	if ((sortord == NULL) | (rightrndx == NULL) |
152	>	(leftrndx == NULL) | (endrndx == NULL))
153		return(0);
154		/* run through bidirections */
155		for (bd = 0; bd < NI2DIR/2; bd++) {
156		const double        ang = 2.*PI/NI2DIR*bd;
157	<	double              cosd, sind;
157	>	int                 *sptr;
158		int                 i;
159		/* create right reverse index */
160	<	if (bd) {                   /* re-use from prev. iteration? */
161	<	int     *sptr = rightrndx;
160	>	if (bd) {                   /* re-use from previous iteration? */
161	>	sptr = rightrndx;
162		rightrndx = leftrndx;
163		leftrndx = sptr;
164	<	} else {                    /* else compute it */
165	<	cosd = cos(ang + (PI/2. - PI/NI2DIR));
166	<	sind = sin(ang + (PI/2. - PI/NI2DIR));
126	<	for (i = 0; i < ip->ns; i++) {
127	<	sortord[i].si = i;
128	<	sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
129	<	}
130	<	qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
131	<	for (i = 0; i < ip->ns; i++)
164	>	} else {                    /* else sort first half-plane */
165	>	sort_samples(sortord, ip, PI/2. - PI/NI2DIR);
166	>	for (i = ip->ns; i--; )
167		rightrndx[sortord[i].si] = i;
168	+	/* & store reverse order for later */
169	+	for (i = ip->ns; i--; )
170	+	endrndx[sortord[i].si] = ip->ns-1 - i;
171		}
172		/* create new left reverse index */
173	<	cosd = cos(ang + (PI/2. + PI/NI2DIR));
174	<	sind = sin(ang + (PI/2. + PI/NI2DIR));
175	<	for (i = 0; i < ip->ns; i++) {
176	<	sortord[i].si = i;
177	<	sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
178	<	}
179	<	qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
142	<	for (i = 0; i < ip->ns; i++)
173	>	if (bd == NI2DIR/2 - 1) {   /* use order from first iteration? */
174	>	sptr = leftrndx;
175	>	leftrndx = endrndx;
176	>	endrndx = sptr;
177	>	} else {                    /* else compute new half-plane */
178	>	sort_samples(sortord, ip, ang + (PI/2. + PI/NI2DIR));
179	>	for (i = ip->ns; i--; )
180		leftrndx[sortord[i].si] = i;
144	–	/* sort grid values in this direction */
145	–	cosd = cos(ang);
146	–	sind = sin(ang);
147	–	for (i = 0; i < ip->ns; i++) {
148	–	sortord[i].si = i;
149	–	sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
181		}
182	<	qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
182	>	/* sort grid values in this direction */
183	>	sort_samples(sortord, ip, ang);
184		/* find nearest neighbors each side */
185	<	for (i = 0; i < ip->ns; i++) {
186	<	const int       rpos = rightrndx[sortord[i].si];
155	<	const int       lpos = leftrndx[sortord[i].si];
185	>	for (i = ip->ns; i--; ) {
186	>	const int       ii = sortord[i].si;
187		int             j;
188	<	/* preload with large radius */
189	<	ip->ra[i][bd] = ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
188	>	/* preload with large radii */
189	>	ip->ra[ii][bd] = ip->ra[ii][bd+NI2DIR/2] = encode_radius(ip,
190		.25*(sortord[ip->ns-1].dm - sortord[0].dm));
191		for (j = i; ++j < ip->ns; )     /* nearest above */
192	<	if (rightrndx[sortord[j].si] > rpos &&
193	<	leftrndx[sortord[j].si] < lpos) {
194	<	ip->ra[i][bd] = encode_radius(ip,
192	>	if (rightrndx[sortord[j].si] > rightrndx[ii] &&
193	>	leftrndx[sortord[j].si] < leftrndx[ii]) {
194	>	ip->ra[ii][bd] = encode_radius(ip,
195		.5*(sortord[j].dm - sortord[i].dm));
196		break;
197		}
198		for (j = i; j-- > 0; )          /* nearest below */
199	<	if (rightrndx[sortord[j].si] < rpos &&
200	<	leftrndx[sortord[j].si] > lpos) {
201	<	ip->ra[i][bd+NI2DIR/2] = encode_radius(ip,
199	>	if (rightrndx[sortord[j].si] < rightrndx[ii] &&
200	>	leftrndx[sortord[j].si] > leftrndx[ii]) {
201	>	ip->ra[ii][bd+NI2DIR/2] = encode_radius(ip,
202		.5*(sortord[i].dm - sortord[j].dm));
203		break;
204		}
#	Line 176 | Line 207 | interp2_compute(INTERP2 *ip)
207		free(sortord);                  /* clean up */
208		free(rightrndx);
209		free(leftrndx);
210	+	free(endrndx);
211		return(1);
212		}
213
214	<	/* private call returns log of raw weight for a particular sample */
214	>	/* private call returns raw weight for a particular sample */
215		static double
216	<	get_ln_wt(const INTERP2 *ip, const int i, double x, double y)
216	>	get_wt(const INTERP2 *ip, const int i, double x, double y)
217		{
218	<	double  dir, rd;
218	>	double  dir, rd, d2;
219		int     ri;
220		/* get relative direction */
221		x -= ip->spt[i][0];
#	Line 196 | Line 228 | get_ln_wt(const INTERP2 *ip, const int i, double x, do
228		rd -= (double)ri;
229		rd = (1.-rd)*ip->ra[i][ri] + rd*ip->ra[i][(ri+1)%NI2DIR];
230		rd = ip->smf * DECODE_RAD(ip, rd);
231	<	/* return log of Gaussian weight */
232	<	return( (x*x + y*y) / (-2.*rd*rd) );
231	>	d2 = x*x + y*y;
232	>	/* Gaussian times harmonic weighting */
233	>	return( exp(d2/(-2.*rd*rd)) * ip->rmin/(ip->rmin + sqrt(d2)) );
234		}
235
236		/* Assign full set of normalized weights to interpolate the given position */
#	Line 210 | Line 243 | interp2_weights(float wtv[], INTERP2 *ip, double x, do
243		if ((wtv == NULL) | (ip == NULL))
244		return(0);
245		/* need to compute interpolant? */
246	<	if (ip->ra == NULL && !interp2_compute(ip))
246	>	if (ip->ra == NULL && !interp2_analyze(ip))
247		return(0);
248
249		wnorm = 0;                      /* compute raw weights */
250		for (i = ip->ns; i--; ) {
251	<	double  wt = get_ln_wt(ip, i, x, y);
219	<	if (wt < -21.) {
220	<	wtv[i] = 0;     /* ignore weights < 1e-9 */
221	<	continue;
222	<	}
223	<	wt = exp(wt);           /* Gaussian weight */
251	>	double  wt = get_wt(ip, i, x, y);
252		wtv[i] = wt;
253		wnorm += wt;
254		}
#	Line 244 | Line 272 | interp2_topsamp(float wt[], int si[], const int n, INT
272		if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
273		return(0);
274		/* need to compute interpolant? */
275	<	if (ip->ra == NULL && !interp2_compute(ip))
275	>	if (ip->ra == NULL && !interp2_analyze(ip))
276		return(0);
277		/* identify top n weights */
278		for (i = ip->ns; i--; ) {
279	<	const double    lnwt = get_ln_wt(ip, i, x, y);
279	>	const double    wti = get_wt(ip, i, x, y);
280		for (j = nn; j > 0; j--) {
281	<	if (wt[j-1] >= lnwt)
281	>	if (wt[j-1] >= wti)
282		break;
283		if (j < n) {
284		wt[j] = wt[j-1];
#	Line 258 | Line 286 | interp2_topsamp(float wt[], int si[], const int n, INT
286		}
287		}
288		if (j < n) {            /* add/insert sample */
289	<	wt[j] = lnwt;
289	>	wt[j] = wti;
290		si[j] = i;
291		nn += (nn < n);
292		}
293		}
294	<	wnorm = 0;                      /* exponentiate and normalize */
295	<	for (j = nn; j--; ) {
296	<	double  dwt = exp(wt[j]);
269	<	wt[j] = dwt;
270	<	wnorm += dwt;
271	<	}
294	>	wnorm = 0;                      /* normalize sample weights */
295	>	for (j = nn; j--; )
296	>	wnorm += wt[j];
297		if (wnorm <= 0)
298		return(0);
299		wnorm = 1./wnorm;

Diff Legend

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