1 /**********
2 *
3 * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
4 * University of Manchester, UK. All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without modification,
7 * are permitted provided that the following conditions are met:
8 *
9 * . Redistributions of source code must retain the above copyright notice,
10 * this list of conditions and the following disclaimer.
11 *
12 * . Redistributions in binary form must reproduce the above copyright notice,
13 * this list of conditions and the following disclaimer in the documentation
14 * and/or other materials provided with the distribution.
15 *
16 * . Neither the name of the University of Manchester nor the names of its
17 * contributors may be used to endorse or promote products derived from this
18 * software without specific prior written permission.
19 *
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
25 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 * POSSIBILITY OF SUCH DAMAGE.
32 *
33 **********
34 *
35 * Program : TINA
36 * File : $Source: /home/tina/cvs/tina-libs/tina/image/imgGen_get_interp.c,v $
37 * Date : $Date: 2003/09/30 16:53:49 $
38 * Version : $Revision: 1.7 $
39 * CVS Id : $Id: imgGen_get_interp.c,v 1.7 2003/09/30 16:53:49 ipoole Exp $
40 */
41
42 /**
43 * @file
44 * @brief Image pixel access by bi-quadratic interpolation.
45 *
46 * There are two versions - one gives continuous
47 * interpolation at tissue boundaries, while Sinc interpolation eliminates
48 * first order truncation error. Sinc is the theoretical optimum for data acquired
49 * in the Fourier domain and gives better results than bi-quadratic which lacks differential
50 * continuity.
51 *
52 */
53
54 #include "imgGen_get_interp.h"
55
56 #if HAVE_CONFIG_H
57 #include <config.h>
58 #endif
59
60 #include <stdio.h>
61 #include <math.h>
62 #include <tina/sys/sysDef.h>
63 #include <tina/sys/sysPro.h>
64 #include <tina/math/mathDef.h>
65 #include <tina/math/mathPro.h>
66 #include <tina/image/img_GenDef.h>
67 #include <tina/image/imgGen_get.h>
68
69
70 double im_get_quadmaxf(Imrect * image, float x, float y, float *px,
71 float *py)
72 /* image of type float */
73 /*
74 * fits simple quadratic function to image and returns interpolated value and
75 * max/min of local area NAT 10/10/90
76 */
77 {
78 Imregion *region = image->region;
79 double pixval[3][3];
80 double a, b, c, d, e, f;
81 double inter;
82 double temp;
83 int i, j, n, m;
84 double xs, ys;
85
86 i = tina_int(x - 1);
87 j = tina_int(y - 1);
88
89 if (j < region->ly + 1 || j > region->uy - 1
90 || i < region->lx + 1 || i > region->ux - 1)
91 return (0);
92
93 for (n = 0; n < 3; n++)
94 {
95 for (m = 0; m < 3; m++)
96 {
97 pixval[n][m] = im_get_pixf(image, j + n, i + m);
98 }
99 }
100
101 a = pixval[1][1];
102 b = (pixval[0][2] - pixval[0][0]
103 + pixval[1][2] - pixval[1][0] + pixval[2][2] - pixval[2][0]) / 6.0;
104 c = (pixval[2][0] - pixval[0][0]
105 + pixval[2][1] - pixval[0][1] + pixval[2][2] - pixval[0][2]) / 6.0;
106 d = (pixval[0][0] - 2.0 * pixval[0][1] + pixval[0][2]
107 + 3.0 * pixval[1][0] - 6.0 * pixval[1][1] + 3.0 * pixval[1][2]
108 + pixval[2][0] - 2.0 * pixval[2][1] + pixval[2][2]) / 10.0;
109 e = (pixval[0][0] - pixval[2][0] + pixval[2][2] - pixval[0][2]) / 4.0;
110 f = (pixval[0][0] + 3.0 * pixval[0][1] + pixval[0][2]
111 - 2.0 * pixval[1][0] - 6.0 * pixval[1][1] - 2.0 * pixval[1][2]
112 + pixval[2][0] + 3.0 * pixval[2][1] + pixval[2][2]) / 10.0;
113
114 temp = 4.0 * d * f - e * e;
115 *px = 1.5 + (double) i + (e * c - 2 * f * b) / temp;
116 *py = 1.5 + (double) j + (e * b - 2 * d * c) / temp;
117
118 xs = (e * c - 2 * f * b) / temp;
119 ys = (e * b - 2 * d * c) / temp;
120
121 inter = a + b * xs + c * ys + d * xs * xs + e * xs * ys + f * ys * ys;
122
123 return (inter);
124 }
125
126 double im_get_quadinterpf(Imrect * image, float x, float y, float *pdx,
127 float *pdy)
128 /* image of type float */
129 /*
130 * fits simple quadratic function to image and returns interpolated value NAT
131 * 10/10/90
132 */
133 {
134 Imregion *region = image->region;
135 double pixval[3][3];
136 double a, b, c, d, e, f;
137 double inter;
138 int i, j, n, m;
139 double xs, ys;
140
141 i = tina_int(x - 1);
142 j = tina_int(y - 1);
143
144 /* assume that pixel contents are exaluated at the centre */
145 xs = x - (float) i - 1.5;
146 ys = y - (float) j - 1.5;
147
148 if (j < region->ly + 1 || j > region->uy - 1
149 || i < region->lx + 1 || i > region->ux - 1)
150 return (0);
151
152 for (n = 0; n < 3; n++)
153 {
154 for (m = 0; m < 3; m++)
155 {
156 pixval[n][m] = im_get_pixf(image, j + n, i + m);
157 }
158 }
159
160 a = pixval[1][1];
161 b = (pixval[0][2] - pixval[0][0]
162 + pixval[1][2] - pixval[1][0] + pixval[2][2] - pixval[2][0]) / 6.0;
163 c = (pixval[2][0] - pixval[0][0]
164 + pixval[2][1] - pixval[0][1] + pixval[2][2] - pixval[0][2]) / 6.0;
165 d = (pixval[0][0] - 2.0 * pixval[0][1] + pixval[0][2]
166 + 3.0 * pixval[1][0] - 6.0 * pixval[1][1] + 3.0 * pixval[1][2]
167 + pixval[2][0] - 2.0 * pixval[2][1] + pixval[2][2]) / 10.0;
168 e = (pixval[0][0] - pixval[2][0] + pixval[2][2] - pixval[0][2]) / 4.0;
169 f = (pixval[0][0] + 3.0 * pixval[0][1] + pixval[0][2]
170 - 2.0 * pixval[1][0] - 6.0 * pixval[1][1] - 2.0 * pixval[1][2]
171 + pixval[2][0] + 3.0 * pixval[2][1] + pixval[2][2]) / 10.0;
172
173 *pdx = b + 2.0 * d * xs + e * ys;
174 *pdy = c + e * xs + 2.0 * f * ys;
175
176 inter = a + b * xs + c * ys + d * xs * xs + e * xs * ys + f * ys * ys;
177
178 return (inter);
179 }
180
181 /* sinc image interpolation using a 5*5 kernel */
182
183 #define KERN_RANGE 2 /* 1 = 3*3 kernel, 2 = 5*5 kernel, 3 = 7*7
184 * kernel */
185
186 static float han_sinc5(float delta)
187 {
188 double sinc, phase, hanning;
189
190 phase = PI * delta;
191 if (fabs(phase) > 0.000001)
192 sinc = sin(phase) / phase;
193 else
194 sinc = 1.0;
195 hanning = (0.5 + 0.5 * cos(phase / (float) (1.0 + KERN_RANGE)));
196 return (hanning * sinc);
197 }
198
199 double im_get_sinc5interpf(Imrect * image, float x, float y,
200 float *pdx, float *pdy)
201 {
202
203 static float *hx = NULL, *hy, *hxstart; /* static data! */
204 float hxtot = 0.0, hytot = 0.0;
205 float ypos, pfac, newpix = 0.0;
206 float *thispix, *hxptr;
207 float xoffset, yoffset;
208 int pos, rowoffset, xint, yint;
209
210 xint = tina_int(x);
211 yint = tina_int(y);
212
213 if (yint < image->region->ly + KERN_RANGE ||
214 yint > image->region->uy - KERN_RANGE - 1 ||
215 xint < image->region->lx + KERN_RANGE ||
216 xint > image->region->ux - KERN_RANGE - 1)
217 return (0.0);
218
219 xoffset = x - (float) xint - 0.5;
220 yoffset = y - (float) yint - 0.5;
221
222 if (hx == NULL)
223 {
224 hx = fvector_alloc(-KERN_RANGE, KERN_RANGE + 1);
225 hy = fvector_alloc(-KERN_RANGE, KERN_RANGE + 1);
226 hxstart = &hx[-KERN_RANGE];
227 }
228 for (pos = -KERN_RANGE; pos <= KERN_RANGE; pos++)
229 {
230 hy[pos] = han_sinc5(yoffset - (float) pos);
231 hytot += hy[pos];
232 hx[pos] = han_sinc5(xoffset - (float) pos);
233 hxtot += hx[pos];
234 }
235
236 rowoffset = tina_int(x - (float) KERN_RANGE);
237
238 for (pos = -KERN_RANGE; pos <= KERN_RANGE; pos++)
239 {
240 ypos = y + (float) pos;
241 thispix = &((float **) image->data)[tina_int(ypos)][rowoffset];
242 hxptr = hxstart;
243 pfac = 0.0;
244 switch (KERN_RANGE)
245 {
246 /* note: no break statements at case ends */
247 case 6:
248 pfac += *(thispix++) * *(hxptr++);
249 pfac += *(thispix++) * *(hxptr++);
250 case 5:
251 pfac += *(thispix++) * *(hxptr++);
252 pfac += *(thispix++) * *(hxptr++);
253 case 4:
254 pfac += *(thispix++) * *(hxptr++);
255 pfac += *(thispix++) * *(hxptr++);
256 case 3:
257 pfac += *(thispix++) * *(hxptr++);
258 pfac += *(thispix++) * *(hxptr++);
259 case 2:
260 pfac += *(thispix++) * *(hxptr++);
261 pfac += *(thispix++) * *(hxptr++);
262 case 1:
263 pfac += *(thispix++) * *(hxptr++);
264 pfac += *(thispix++) * *(hxptr++);
265 pfac += *(thispix++) * *(hxptr++);
266 }
267 newpix += hy[pos] * pfac;
268 }
269
270 return (newpix / (hxtot * hytot));
271 }
272
273 #undef KERN_RANGE
274
275 /* sinc image interpolation using a 3*3 kernel */
276
277 #define KERN_RANGE 1 /* 1 = 3*3 kernel, 2 = 5*5 kernel, 3 = 7*7
278 * kernel */
279 /*
280 * UNUSED
281 static float han_sinc3(float delta)
282 {
283 double sinc, phase, hanning;
284
285 phase = PI * delta;
286 if (fabs(phase) > 0.000001)
287 sinc = sin(phase) / phase;
288 else
289 sinc = 1.0;
290 hanning = (0.5 + 0.5 * cos(phase / (float) (1.0 + KERN_RANGE)));
291 return (hanning * sinc);
292 }
293 */
294
295 double im_get_sinc3interpf(Imrect * image, float x, float y,
296 float *pdx, float *pdy)
297 {
298
299 static float *hx = NULL, *hy, *hxstart; /* static data! */
300 float hxtot = 0.0, hytot = 0.0;
301 float ypos, pfac, newpix = 0.0;
302 float *thispix, *hxptr;
303 float xoffset, yoffset;
304 int pos, rowoffset, xint, yint;
305
306 xint = tina_int(x);
307 yint = tina_int(y);
308
309 if (yint < image->region->ly + KERN_RANGE ||
310 yint > image->region->uy - KERN_RANGE - 1 ||
311 xint < image->region->lx + KERN_RANGE ||
312 xint > image->region->ux - KERN_RANGE - 1)
313 return (0.0);
314
315 xoffset = x - (float) xint - 0.5;
316 yoffset = y - (float) yint - 0.5;
317
318 if (hx == NULL)
319 {
320 hx = fvector_alloc(-KERN_RANGE, KERN_RANGE + 1);
321 hy = fvector_alloc(-KERN_RANGE, KERN_RANGE + 1);
322 hxstart = &hx[-KERN_RANGE];
323 }
324 for (pos = -KERN_RANGE; pos <= KERN_RANGE; pos++)
325 {
326 hy[pos] = han_sinc5(yoffset - (float) pos);
327 hytot += hy[pos];
328 hx[pos] = han_sinc5(xoffset - (float) pos);
329 hxtot += hx[pos];
330 }
331
332 rowoffset = tina_int(x - (float) KERN_RANGE);
333
334 for (pos = -KERN_RANGE; pos <= KERN_RANGE; pos++)
335 {
336 ypos = y + (float) pos;
337 thispix = &((float **) image->data)[tina_int(ypos)][rowoffset];
338 hxptr = hxstart;
339 pfac = 0.0;
340 switch (KERN_RANGE)
341 {
342 /* note: no break statements at case ends */
343 case 6:
344 pfac += *(thispix++) * *(hxptr++);
345 pfac += *(thispix++) * *(hxptr++);
346 case 5:
347 pfac += *(thispix++) * *(hxptr++);
348 pfac += *(thispix++) * *(hxptr++);
349 case 4:
350 pfac += *(thispix++) * *(hxptr++);
351 pfac += *(thispix++) * *(hxptr++);
352 case 3:
353 pfac += *(thispix++) * *(hxptr++);
354 pfac += *(thispix++) * *(hxptr++);
355 case 2:
356 pfac += *(thispix++) * *(hxptr++);
357 pfac += *(thispix++) * *(hxptr++);
358 case 1:
359 pfac += *(thispix++) * *(hxptr++);
360 pfac += *(thispix++) * *(hxptr++);
361 pfac += *(thispix++) * *(hxptr++);
362 }
363 newpix += hy[pos] * pfac;
364 }
365
366 return (newpix / (hxtot * hytot));
367 }
368
369 #undef KERN_RANGE
370
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.