~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina5/tina-libs/tina/image/imgGen_get_interp.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  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 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.