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.5 by greg, Mon Feb 11 23:33:35 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	>	* direction.  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.
27	>	****************************************************************/
28
29		#include <stdio.h>
30		#include <stdlib.h>
31		#include "rtmath.h"
32		#include "interp2d.h"
33
34	<	#define DECODE_RAD(ip,er)       ((ip)->rmin*(1. + .5*(er)))
35	<	#define ENCODE_RAD(ip,r)        ((int)(2.*(r)/(ip)->rmin) - 2)
34	>	#define DECODE_DIA(ip,ed)       ((ip)->dmin*(1. + .5*(ed)))
35	>	#define ENCODE_DIA(ip,d)        ((int)(2.*(d)/(ip)->dmin) - 2)
36
37		/* Sample order (private) */
38		typedef struct {
#	Line 54 | Line 54 | interp2_alloc(int nsamps)
54		return(NULL);
55
56		nip->ns = nsamps;
57	<	nip->rmin = .5;         /* default radius minimum */
57	>	nip->dmin = 1;          /* default minimum diameter */
58		nip->smf = NI2DSMF;     /* default smoothing factor */
59	<	nip->ra = NULL;
59	>	nip->da = NULL;
60		/* caller must assign spt[] array */
61		return(nip);
62		}
#	Line 73 | Line 73 | interp2_realloc(INTERP2 *ip, int nsamps)
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;
76	>	if (ip->da != NULL) {   /* will need to recompute distribution */
77	>	free(ip->da);
78	>	ip->da = NULL;
79		}
80		ip = (INTERP2 *)realloc(ip, sizeof(INTERP2)+sizeof(float)*2*(nsamps-1));
81		if (ip == NULL)
#	Line 84 | Line 84 | interp2_realloc(INTERP2 *ip, int nsamps)
84		return(ip);
85		}
86
87	+	/* Set minimum distance under which samples will start to merge */
88	+	void
89	+	interp2_spacing(INTERP2 *ip, double mind)
90	+	{
91	+	if (mind <= 0)
92	+	return;
93	+	if ((.998*ip->dmin <= mind) && (mind <= 1.002*ip->dmin))
94	+	return;
95	+	if (ip->da != NULL) {   /* will need to recompute distribution */
96	+	free(ip->da);
97	+	ip->da = NULL;
98	+	}
99	+	ip->dmin = mind;
100	+	}
101	+
102	+	/* Modify smoothing parameter by the given factor */
103	+	void
104	+	interp2_smooth(INTERP2 *ip, double sf)
105	+	{
106	+	if ((ip->smf *= sf) < NI2DSMF)
107	+	ip->smf = NI2DSMF;
108	+	}
109	+
110		/* private call-back to sort position index */
111		static int
112		cmp_spos(const void *p1, const void *p2)
#	Line 113 | Line 136 | sort_samples(SAMPORD *sord, const INTERP2 *ip, double
136		qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
137		}
138
139	<	/* private routine to encode radius with range checks */
139	>	/* private routine to encode sample diameter with range checks */
140		static int
141	<	encode_radius(const INTERP2 *ip, double r)
141	>	encode_diameter(const INTERP2 *ip, double d)
142		{
143	<	const int       er = ENCODE_RAD(ip, r);
143	>	const int       ed = ENCODE_DIA(ip, d);
144
145	<	if (er <= 0)
145	>	if (ed <= 0)
146		return(0);
147	<	if (er >= 0xffff)
147	>	if (ed >= 0xffff)
148		return(0xffff);
149	<	return(er);
149	>	return(ed);
150		}
151
152		/* (Re)compute anisotropic basis function interpolant (normally automatic) */
#	Line 134 | Line 157 | interp2_analyze(INTERP2 *ip)
157		int     *rightrndx, *leftrndx, *endrndx;
158		int     bd;
159		/* sanity checks */
160	<	if (ip == NULL || (ip->ns <= 1) | (ip->rmin <= 0))
160	>	if (ip == NULL || (ip->ns <= 1) | (ip->dmin <= 0))
161		return(0);
162		/* need to allocate? */
163	<	if (ip->ra == NULL) {
164	<	ip->ra = (unsigned short (*)[NI2DIR])malloc(
163	>	if (ip->da == NULL) {
164	>	ip->da = (unsigned short (*)[NI2DIR])malloc(
165		sizeof(unsigned short)*NI2DIR*ip->ns);
166	<	if (ip->ra == NULL)
166	>	if (ip->da == NULL)
167		return(0);
168		}
169		/* get temporary arrays */
#	Line 186 | Line 209 | interp2_analyze(INTERP2 *ip)
209		const int       ii = sortord[i].si;
210		int             j;
211		/* preload with large radii */
212	<	ip->ra[ii][bd] = ip->ra[ii][bd+NI2DIR/2] = encode_radius(ip,
213	<	.25*(sortord[ip->ns-1].dm - sortord[0].dm));
212	>	ip->da[ii][bd] = ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
213	>	.5*(sortord[ip->ns-1].dm - sortord[0].dm));
214		for (j = i; ++j < ip->ns; )     /* nearest above */
215		if (rightrndx[sortord[j].si] > rightrndx[ii] &&
216		leftrndx[sortord[j].si] < leftrndx[ii]) {
217	<	ip->ra[ii][bd] = encode_radius(ip,
218	<	.5*(sortord[j].dm - sortord[i].dm));
217	>	ip->da[ii][bd] = encode_diameter(ip,
218	>	sortord[j].dm - sortord[i].dm);
219		break;
220		}
221		for (j = i; j-- > 0; )          /* nearest below */
222		if (rightrndx[sortord[j].si] < rightrndx[ii] &&
223		leftrndx[sortord[j].si] > leftrndx[ii]) {
224	<	ip->ra[ii][bd+NI2DIR/2] = encode_radius(ip,
225	<	.5*(sortord[i].dm - sortord[j].dm));
224	>	ip->da[ii][bd+NI2DIR/2] = encode_diameter(ip,
225	>	sortord[i].dm - sortord[j].dm);
226		break;
227		}
228		}
#	Line 226 | Line 249 | get_wt(const INTERP2 *ip, const int i, double x, doubl
249		rd = dir * (NI2DIR/2./PI);
250		ri = (int)rd;
251		rd -= (double)ri;
252	<	rd = (1.-rd)*ip->ra[i][ri] + rd*ip->ra[i][(ri+1)%NI2DIR];
253	<	rd = ip->smf * DECODE_RAD(ip, rd);
252	>	rd = (1.-rd)*ip->da[i][ri] + rd*ip->da[i][(ri+1)%NI2DIR];
253	>	rd = ip->smf * DECODE_DIA(ip, rd);
254		d2 = x*x + y*y;
255		/* Gaussian times harmonic weighting */
256	<	return( exp(d2/(-2.*rd*rd)) * ip->rmin/(ip->rmin + sqrt(d2)) );
256	>	return( exp(d2/(-2.*rd*rd)) * ip->dmin/(ip->dmin + sqrt(d2)) );
257		}
258
259		/* Assign full set of normalized weights to interpolate the given position */
#	Line 243 | Line 266 | interp2_weights(float wtv[], INTERP2 *ip, double x, do
266		if ((wtv == NULL) | (ip == NULL))
267		return(0);
268		/* need to compute interpolant? */
269	<	if (ip->ra == NULL && !interp2_analyze(ip))
269	>	if (ip->da == NULL && !interp2_analyze(ip))
270		return(0);
271
272		wnorm = 0;                      /* compute raw weights */
#	Line 272 | Line 295 | interp2_topsamp(float wt[], int si[], const int n, INT
295		if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
296		return(0);
297		/* need to compute interpolant? */
298	<	if (ip->ra == NULL && !interp2_analyze(ip))
298	>	if (ip->da == NULL && !interp2_analyze(ip))
299		return(0);
300		/* identify top n weights */
301		for (i = ip->ns; i--; ) {
#	Line 308 | Line 331 | interp2_free(INTERP2 *ip)
331		{
332		if (ip == NULL)
333		return;
334	<	if (ip->ra != NULL)
335	<	free(ip->ra);
334	>	if (ip->da != NULL)
335	>	free(ip->da);
336		free(ip);
337		}

Diff Legend

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