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.2 by greg, Sat Feb 9 17:39:21 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 182 | Line 183 | interp2_analyze(INTERP2 *ip)
183		sort_samples(sortord, ip, ang);
184		/* find nearest neighbors each side */
185		for (i = ip->ns; i--; ) {
186	<	const int       rpos = rightrndx[sortord[i].si];
186	<	const int       lpos = leftrndx[sortord[i].si];
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 211 | Line 211 | interp2_analyze(INTERP2 *ip)
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 228 | 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 247 | Line 248 | interp2_weights(float wtv[], INTERP2 *ip, double x, do
248
249		wnorm = 0;                      /* compute raw weights */
250		for (i = ip->ns; i--; ) {
251	<	double  wt = get_ln_wt(ip, i, x, y);
251	<	if (wt < -21.) {
252	<	wtv[i] = 0;     /* ignore weights < 1e-9 */
253	<	continue;
254	<	}
255	<	wt = exp(wt);           /* Gaussian weight */
251	>	double  wt = get_wt(ip, i, x, y);
252		wtv[i] = wt;
253		wnorm += wt;
254		}
#	Line 280 | Line 276 | interp2_topsamp(float wt[], int si[], const int n, INT
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 290 | 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]);
301	<	wt[j] = dwt;
302	<	wnorm += dwt;
303	<	}
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)