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.2 by greg, Sat Feb 9 17:39:21 2013 UTC

#	Line 39 | Line 39 | typedef struct {
39		float   dm;             /* distance measure in this direction */
40		} SAMPORD;
41
42	<	/* Allocate a new set of interpolation samples */
42	>	/* Allocate a new set of interpolation samples (caller assigns spt[] array) */
43		INTERP2 *
44		interp2_alloc(int nsamps)
45		{
#	Line 60 | Line 60 | interp2_alloc(int nsamps)
60		return(nip);
61		}
62
63	+	/* Resize interpolation array (caller must assign any new values) */
64	+	INTERP2 *
65	+	interp2_realloc(INTERP2 *ip, int nsamps)
66	+	{
67	+	if (ip == NULL)
68	+	return(interp2_alloc(nsamps));
69	+	if (nsamps <= 1) {
70	+	interp2_free(ip);
71	+	return(NULL);
72	+	}
73	+	if (nsamps == ip->ns);
74	+	return(ip);
75	+	if (ip->ra != NULL) {   /* will need to recompute distribution */
76	+	free(ip->ra);
77	+	ip->ra = NULL;
78	+	}
79	+	ip = (INTERP2 *)realloc(ip, sizeof(INTERP2)+sizeof(float)*2*(nsamps-1));
80	+	if (ip == NULL)
81	+	return(NULL);
82	+	ip->ns = nsamps;
83	+	return(ip);
84	+	}
85	+
86		/* private call-back to sort position index */
87		static int
88		cmp_spos(const void *p1, const void *p2)
#	Line 74 | Line 97 | cmp_spos(const void *p1, const void *p2)
97		return 0;
98		}
99
100	+	/* private routine to order samples in a particular direction */
101	+	static void
102	+	sort_samples(SAMPORD *sord, const INTERP2 *ip, double ang)
103	+	{
104	+	const double    cosd = cos(ang);
105	+	const double    sind = sin(ang);
106	+	int             i;
107	+
108	+	for (i = ip->ns; i--; ) {
109	+	sord[i].si = i;
110	+	sord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
111	+	}
112	+	qsort(sord, ip->ns, sizeof(SAMPORD), &cmp_spos);
113	+	}
114	+
115		/* private routine to encode radius with range checks */
116		static int
117		encode_radius(const INTERP2 *ip, double r)
#	Line 87 | Line 125 | encode_radius(const INTERP2 *ip, double r)
125		return(er);
126		}
127
128	<	/* Compute anisotropic Gaussian basis function interpolant */
129	<	static int
130	<	interp2_compute(INTERP2 *ip)
128	>	/* (Re)compute anisotropic basis function interpolant (normally automatic) */
129	>	int
130	>	interp2_analyze(INTERP2 *ip)
131		{
132		SAMPORD *sortord;
133	<	int     *rightrndx, *leftrndx;
133	>	int     *rightrndx, *leftrndx, *endrndx;
134		int     bd;
135		/* sanity checks */
136		if (ip == NULL || (ip->ns <= 1) | (ip->rmin <= 0))
#	Line 108 | Line 146 | interp2_compute(INTERP2 *ip)
146		sortord = (SAMPORD *)malloc(sizeof(SAMPORD)*ip->ns);
147		rightrndx = (int *)malloc(sizeof(int)*ip->ns);
148		leftrndx = (int *)malloc(sizeof(int)*ip->ns);
149	<	if ((sortord == NULL) | (rightrndx == NULL) | (leftrndx == NULL))
149	>	endrndx = (int *)malloc(sizeof(int)*ip->ns);
150	>	if ((sortord == NULL) | (rightrndx == NULL) |
151	>	(leftrndx == NULL) | (endrndx == NULL))
152		return(0);
153		/* run through bidirections */
154		for (bd = 0; bd < NI2DIR/2; bd++) {
155		const double        ang = 2.*PI/NI2DIR*bd;
156	<	double              cosd, sind;
156	>	int                 *sptr;
157		int                 i;
158		/* create right reverse index */
159	<	if (bd) {                   /* re-use from prev. iteration? */
160	<	int     *sptr = rightrndx;
159	>	if (bd) {                   /* re-use from previous iteration? */
160	>	sptr = rightrndx;
161		rightrndx = leftrndx;
162		leftrndx = sptr;
163	<	} else {                    /* else compute it */
164	<	cosd = cos(ang + (PI/2. - PI/NI2DIR));
165	<	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++)
163	>	} else {                    /* else sort first half-plane */
164	>	sort_samples(sortord, ip, PI/2. - PI/NI2DIR);
165	>	for (i = ip->ns; i--; )
166		rightrndx[sortord[i].si] = i;
167	+	/* & store reverse order for later */
168	+	for (i = ip->ns; i--; )
169	+	endrndx[sortord[i].si] = ip->ns-1 - i;
170		}
171		/* create new left reverse index */
172	<	cosd = cos(ang + (PI/2. + PI/NI2DIR));
173	<	sind = sin(ang + (PI/2. + PI/NI2DIR));
174	<	for (i = 0; i < ip->ns; i++) {
175	<	sortord[i].si = i;
176	<	sortord[i].dm = cosd*ip->spt[i][0] + sind*ip->spt[i][1];
177	<	}
178	<	qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
142	<	for (i = 0; i < ip->ns; i++)
172	>	if (bd == NI2DIR/2 - 1) {   /* use order from first iteration? */
173	>	sptr = leftrndx;
174	>	leftrndx = endrndx;
175	>	endrndx = sptr;
176	>	} else {                    /* else compute new half-plane */
177	>	sort_samples(sortord, ip, ang + (PI/2. + PI/NI2DIR));
178	>	for (i = ip->ns; i--; )
179		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];
180		}
181	<	qsort(sortord, ip->ns, sizeof(SAMPORD), &cmp_spos);
181	>	/* sort grid values in this direction */
182	>	sort_samples(sortord, ip, ang);
183		/* find nearest neighbors each side */
184	<	for (i = 0; i < ip->ns; i++) {
184	>	for (i = ip->ns; i--; ) {
185		const int       rpos = rightrndx[sortord[i].si];
186		const int       lpos = leftrndx[sortord[i].si];
187		int             j;
#	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
#	Line 210 | Line 242 | interp2_weights(float wtv[], INTERP2 *ip, double x, do
242		if ((wtv == NULL) | (ip == NULL))
243		return(0);
244		/* need to compute interpolant? */
245	<	if (ip->ra == NULL && !interp2_compute(ip))
245	>	if (ip->ra == NULL && !interp2_analyze(ip))
246		return(0);
247
248		wnorm = 0;                      /* compute raw weights */
#	Line 244 | Line 276 | interp2_topsamp(float wt[], int si[], const int n, INT
276		if ((n <= 0) | (wt == NULL) | (si == NULL) | (ip == NULL))
277		return(0);
278		/* need to compute interpolant? */
279	<	if (ip->ra == NULL && !interp2_compute(ip))
279	>	if (ip->ra == NULL && !interp2_analyze(ip))
280		return(0);
281		/* identify top n weights */
282		for (i = ip->ns; i--; ) {

Diff Legend

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