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

Linux Cross Reference
Tina6/tina-libs/tina/image/imgSeq_slice.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/imgSeq_slice.c,v $
 37  * Date    :  $Date: 2004/08/06 15:04:58 $
 38  * Version :  $Revision: 1.10 $
 39  * CVS Id  :  $Id: imgSeq_slice.c,v 1.10 2004/08/06 15:04:58 paul Exp $
 40  */
 41 
 42 /** 
 43  *  @file
 44  *  @brief Functions associated with sub-pixel interpolation of 3D volume data, such as
 45  *  that which can be stored in the sequence data structure. 
 46  *
 47  *  Uses various interpolation approaches including tri-linear and re-normalised Sinc interpolation.
 48  *  Routines can generate images at fixed orthogonal planes for display and other purposes.  
 49  *
 50  */
 51 
 52 
 53 #include "imgSeq_slice.h"
 54 
 55 #if HAVE_CONFIG_H
 56 #include <config.h>
 57 #endif
 58 
 59 #include <stdio.h>
 60 #include <math.h>
 61 #include <tina/sys/sysDef.h>
 62 #include <tina/sys/sysPro.h>
 63 #include <tina/math/mathDef.h>
 64 #include <tina/math/mathPro.h>
 65 #include <tina/image/imgDef.h>
 66 #include <tina/image/imgPro.h>
 67 
 68 
 69 static void ***imptrs=NULL;
 70 static int imptrlx, imptrux,imptrly, imptruy, imptrlz, imptruz;
 71 static int blx, bux, bly, buy, blz, buz;
 72 static float (*interp)(void ***rasptrs, Vec3 pos) = NULL;
 73 static float (*interp_mask)(void ***rasptrs, Vec3 pos, int *mask_pix) = NULL;
 74 static float (*seq_pixel)(void ***rasptrs, int z, int y, int x) = NULL;
 75 static double (*sinc_accum)(void ***rasptrs, int z, int y, int x, float *hxp) = NULL;
 76 static int barrier=1;
 77 static double outside = 0.0;
 78 
 79 
 80 void seq_set_outside(double out)
 81 {
 82     /* The outside value is returned by the interpolation functions when
 83        a voxel outside the limits of the current sequence is requested.
 84        It can be used as a flag to indicate that this has happened by 
 85        setting it to a value that does not occur in the sequence, and 
 86        then testing for voxels of e.g. more than a slightly smaller 
 87        value.  Remeber to set it back to zero after you have used it, 
 88        to prevent interference with the rest of TINA.*/
 89 
 90     outside = out;
 91 }
 92 
 93 void ***seq_limits(int *lz, int *uz, int *ly, int *uy, int *lx, int *ux)
 94 {
 95    *lz = imptrlz;
 96    *uz = imptruz;
 97    *ly = imptrly;
 98    *uy = imptruy;
 99    *lx = imptrlx;
100    *ux = imptrux;
101    return(imptrs);
102 }
103 
104 static float seq_cpixel(void ***rasptrs, int z, int y, int x)
105 {
106     char *row;
107     row = rasptrs[z][y];
108     return((float)row[x]);
109 }
110 
111 static float seq_ucpixel(void ***rasptrs, int z, int y, int x)
112 {
113     unsigned char *row;
114     row = rasptrs[z][y];
115     return((float)row[x]);
116 }
117 
118 static float seq_spixel(void ***rasptrs, int z, int y, int x)
119 {
120     short *row;
121     row = rasptrs[z][y];
122     return((float)row[x]);
123 }
124 
125 static float seq_uspixel(void ***rasptrs, int z, int y, int x)
126 {
127     unsigned short *row;
128     row = rasptrs[z][y];
129     return((float)row[x]);
130 }
131 
132 static float seq_ipixel(void ***rasptrs, int z, int y, int x)
133 {
134     int *row;
135     row = rasptrs[z][y];
136     return((float)row[x]);
137 }
138 
139 static float seq_fpixel(void ***rasptrs,int z, int y, int x)
140 {
141     float *row;
142     row = rasptrs[z][y];
143     return((float)row[x]);
144 }
145 
146 static double sinc_accumuc(void *** rasptrs, int z, int y, int x, float *hxp)
147 {
148      unsigned char *row;
149      unsigned char *thispix;
150      double pfac;
151  
152      row = rasptrs[z][y];
153      thispix = &row[x];
154      pfac = 0.0;
155      switch(barrier)
156      {
157         case 6:
158           pfac += *(thispix++)* *(hxp++);
159           pfac += *(thispix++)* *(hxp++);
160 
161         case 5:
162           pfac += *(thispix++)* *(hxp++);
163           pfac += *(thispix++)* *(hxp++);
164 
165         case 4:
166           pfac += *(thispix++)* *(hxp++);
167           pfac += *(thispix++)* *(hxp++);
168 
169         case 3:
170           pfac += *(thispix++)* *(hxp++);
171           pfac += *(thispix++)* *(hxp++);
172 
173         case 2:
174           pfac += *(thispix++)* *(hxp++);
175           pfac += *(thispix++)* *(hxp++);
176 
177        case 1:
178          pfac += *(thispix++)* *(hxp++);
179          pfac += *(thispix++)* *(hxp++);
180          pfac += *(thispix++)* *(hxp++);
181     }
182     return (pfac);
183 }
184 
185 static double sinc_accumc(void *** rasptrs, int z, int y, int x, float *hxp)
186 {
187      char *row;
188      char *thispix;
189      double pfac;
190 
191      row = rasptrs[z][y];
192      thispix = &row[x];
193      pfac = 0.0;
194      switch(barrier)
195      {
196         case 6:
197           pfac += *(thispix++)* *(hxp++);
198           pfac += *(thispix++)* *(hxp++);
199 
200         case 5:
201           pfac += *(thispix++)* *(hxp++);
202           pfac += *(thispix++)* *(hxp++);
203 
204         case 4:
205           pfac += *(thispix++)* *(hxp++);
206           pfac += *(thispix++)* *(hxp++);
207 
208         case 3:
209           pfac += *(thispix++)* *(hxp++);
210           pfac += *(thispix++)* *(hxp++);
211 
212         case 2:
213           pfac += *(thispix++)* *(hxp++);
214           pfac += *(thispix++)* *(hxp++);
215 
216        case 1:
217          pfac += *(thispix++)* *(hxp++);
218          pfac += *(thispix++)* *(hxp++);
219          pfac += *(thispix++)* *(hxp++);
220     }
221     return (pfac);
222 }
223 
224 static double sinc_accums(void *** rasptrs, int z, int y, int x, float *hxp)
225 {
226      short *row;
227      short *thispix;
228      double pfac;
229 
230      row = rasptrs[z][y];
231      thispix = &row[x];
232      pfac = 0.0;
233      switch(barrier)
234      {
235         case 6:
236           pfac += *(thispix++)* *(hxp++);
237           pfac += *(thispix++)* *(hxp++);
238 
239         case 5:
240           pfac += *(thispix++)* *(hxp++);
241           pfac += *(thispix++)* *(hxp++);
242 
243         case 4:
244           pfac += *(thispix++)* *(hxp++);
245           pfac += *(thispix++)* *(hxp++);
246 
247         case 3:
248           pfac += *(thispix++)* *(hxp++);
249           pfac += *(thispix++)* *(hxp++);
250 
251         case 2:
252           pfac += *(thispix++)* *(hxp++);
253           pfac += *(thispix++)* *(hxp++);
254 
255        case 1:
256          pfac += *(thispix++)* *(hxp++);
257          pfac += *(thispix++)* *(hxp++);
258          pfac += *(thispix++)* *(hxp++);
259     }
260     return (pfac);
261 }
262 
263 static double sinc_accumus(void *** rasptrs, int z, int y, int x, float *hxp)
264 {
265      unsigned short *row;
266      unsigned short *thispix;
267      double pfac;
268 
269      row = rasptrs[z][y];
270      thispix = &row[x];
271      pfac = 0.0;
272      switch(barrier)
273      {
274         case 6:
275           pfac += *(thispix++)* *(hxp++);
276           pfac += *(thispix++)* *(hxp++);
277 
278         case 5:
279           pfac += *(thispix++)* *(hxp++);
280           pfac += *(thispix++)* *(hxp++);
281 
282         case 4:
283           pfac += *(thispix++)* *(hxp++);
284           pfac += *(thispix++)* *(hxp++);
285 
286         case 3:
287           pfac += *(thispix++)* *(hxp++);
288           pfac += *(thispix++)* *(hxp++);
289 
290         case 2:
291           pfac += *(thispix++)* *(hxp++);
292           pfac += *(thispix++)* *(hxp++);
293 
294        case 1:
295          pfac += *(thispix++)* *(hxp++);
296          pfac += *(thispix++)* *(hxp++);
297          pfac += *(thispix++)* *(hxp++);
298     }
299     return (pfac);
300 }
301 
302 static double sinc_accumi(void *** rasptrs, int z, int y, int x, float *hxp)
303 {
304      int *row;
305      int *thispix;
306      double pfac;
307 
308      row = rasptrs[z][y];
309      thispix = &row[x];
310      pfac = 0.0;
311      switch(barrier)
312      {
313         case 6:
314           pfac += *(thispix++)* *(hxp++);
315           pfac += *(thispix++)* *(hxp++);
316 
317         case 5:
318           pfac += *(thispix++)* *(hxp++);
319           pfac += *(thispix++)* *(hxp++);
320 
321         case 4:
322           pfac += *(thispix++)* *(hxp++);
323           pfac += *(thispix++)* *(hxp++);
324 
325         case 3:
326           pfac += *(thispix++)* *(hxp++);
327           pfac += *(thispix++)* *(hxp++);
328 
329         case 2:
330           pfac += *(thispix++)* *(hxp++);
331           pfac += *(thispix++)* *(hxp++);
332 
333        case 1:
334          pfac += *(thispix++)* *(hxp++);
335          pfac += *(thispix++)* *(hxp++);
336          pfac += *(thispix++)* *(hxp++);
337     }
338     return (pfac);
339 }
340 
341 static double sinc_accumf(void *** rasptrs, int z, int y, int x, float *hxp)
342 {
343      float *row;
344      float *thispix;
345      double pfac;
346 
347      row = rasptrs[z][y];
348      thispix = &row[x];
349      pfac = 0.0;
350      switch(barrier)
351      {
352         case 6:
353           pfac += *(thispix++)* *(hxp++);
354           pfac += *(thispix++)* *(hxp++);
355 
356         case 5:
357           pfac += *(thispix++)* *(hxp++);
358           pfac += *(thispix++)* *(hxp++);
359 
360         case 4:
361           pfac += *(thispix++)* *(hxp++);
362           pfac += *(thispix++)* *(hxp++);
363 
364         case 3:
365           pfac += *(thispix++)* *(hxp++);
366           pfac += *(thispix++)* *(hxp++);
367 
368         case 2:
369           pfac += *(thispix++)* *(hxp++);
370           pfac += *(thispix++)* *(hxp++);
371 
372        case 1:
373          pfac += *(thispix++)* *(hxp++);
374          pfac += *(thispix++)* *(hxp++);
375          pfac += *(thispix++)* *(hxp++);
376     }
377     return (pfac);
378 }
379 
380 
381 void seq_voxel_vtype(Vartype vtype)
382 {
383     if (vtype == char_v) seq_pixel = seq_cpixel; 
384     if (vtype == char_v) sinc_accum = sinc_accumc; 
385     if (vtype == uchar_v) seq_pixel = seq_ucpixel; 
386     if (vtype == uchar_v) sinc_accum = sinc_accumuc; 
387     if (vtype == short_v) seq_pixel = seq_spixel; 
388     if (vtype == short_v) sinc_accum = sinc_accums; 
389     if (vtype == ushort_v) seq_pixel = seq_uspixel; 
390     if (vtype == ushort_v) sinc_accum = sinc_accumus; 
391     if (vtype == int_v) seq_pixel = seq_ipixel; 
392     if (vtype == int_v) sinc_accum = sinc_accumi; 
393     if (vtype == float_v) seq_pixel = seq_fpixel; 
394     if (vtype == float_v) sinc_accum = sinc_accumf; 
395 }
396 
397 void seq_init_interp(int nblx, int nbux, int nbly, int nbuy, int nblz, int nbuz)
398 {
399    blx = nblx+barrier;
400    bux = nbux-barrier;
401    bly = nbly+barrier;
402    buy = nbuy-barrier;
403    blz = nblz+barrier;
404    buz = nbuz-barrier;
405 }
406 
407 float this_pixel(void ***rasptrs, int z, int y, int x)
408 {
409    return(seq_pixel(rasptrs,z, y, x));
410 }
411 
412 float seq_interp(void ***rasptrs, Vec3 pos)
413 {
414    return(interp(rasptrs,pos));
415 }
416 
417 float nearest_pixel(void ***rasptrs, Vec3 pos)
418 {
419     int x, y, z;
420     x = tina_int(pos.el[0]);
421     y = tina_int(pos.el[1]);
422     z = tina_int(pos.el[2]);
423     if (x<blx || x>=bux
424     ||  y<bly || y>=buy
425     ||  z<blz || z>=buz )
426        return(outside);
427     return(seq_pixel(rasptrs,z, y, x));
428 }
429 
430 float tri_linear(void ***rasptrs, Vec3 pos)
431 {
432     float xoffset, yoffset, zoffset;
433     float cenpix,newpix,nextpix;
434     float xfac, yfac, zfac;
435     int x , y, z;
436 
437     x = tina_int(pos.el[0]);
438     y = tina_int(pos.el[1]);
439     z = tina_int(pos.el[2]);
440     if (x<blx || x>=bux 
441     ||  y<bly || y>=buy
442     ||  z<blz || z>=buz )
443        return(outside);
444 
445     xoffset = pos.el[0] - x - 0.5;
446     yoffset = pos.el[1] - y - 0.5;
447     zoffset = pos.el[2] - z - 0.5;
448     cenpix = seq_pixel(rasptrs,z,y,x);
449     if (zoffset<0.0) 
450     {
451       nextpix = seq_pixel(rasptrs,z-1,y,x);
452       zfac = cenpix - nextpix;
453     }
454     else
455     {
456       nextpix = seq_pixel(rasptrs,z+1,y,x);
457       zfac = nextpix - cenpix;
458     }
459 
460     if (yoffset<0.0)
461     {
462       nextpix = seq_pixel(rasptrs,z,y-1,x); 
463       yfac = cenpix -nextpix;
464     }
465     else
466     {
467       nextpix = seq_pixel(rasptrs,z,y+1,x); 
468       yfac = nextpix - cenpix;
469     }
470 
471     if (xoffset<0.0)
472     {
473       nextpix = seq_pixel(rasptrs,z,y,x-1);
474       xfac = cenpix - nextpix;
475     }
476     else
477     {
478       nextpix = seq_pixel(rasptrs,z,y,x+1);
479       xfac = nextpix - cenpix;
480     }
481     newpix = cenpix + xfac*xoffset + yfac*yoffset + zfac*zoffset;
482 
483     return(newpix);
484 }
485 
486 static float han_sinc(float delta)
487 {
488    double sinc ,phase,hanning;
489 
490    phase = PI*delta;
491    if(fabs(phase) > 0.000001) sinc = sin(phase)/phase;
492    else sinc = 1.0;
493    hanning = (0.5 + 0.5*cos(phase/(float)(1.0+barrier)));
494    return(hanning*sinc);
495 }
496 
497 static float sinc_interp(void ***rasptrs, Vec3 pos)
498 {
499     float ypos, zpos;
500     float newpix=0.0;
501     float xoffset, yoffset, zoffset;
502     float hxtot=0.0, hytot=0.0, hztot=0.0;
503     float pfac;
504     static float *hx=NULL;                      /* static data! */
505     static float *hy=NULL, *hz =NULL;           /* static data! */
506     int i,j;
507     int x , y, z;
508 
509     x = tina_int(pos.el[0]);
510     y = tina_int(pos.el[1]);
511     z = tina_int(pos.el[2]);
512     if (x<blx || x>=bux 
513     ||  y<bly || y>=buy
514     ||  z<blz || z>=buz )
515        return(outside);
516 
517     xoffset = pos.el[0] - x - 0.5;
518     yoffset = pos.el[1] - y - 0.5;
519     zoffset = pos.el[2] - z - 0.5;
520     if (hx  == NULL)
521     {
522         hx = fvector_alloc(-6,7);
523         hy = fvector_alloc(-6,7);
524         hz = fvector_alloc(-6,7);
525     }
526     for (i=-barrier;i<=barrier;i++)
527     {
528        hz[i] = han_sinc(zoffset-(float)i);
529        hztot += hz[i];
530        hy[i] = han_sinc(yoffset-(float)i);
531        hytot += hy[i];
532        hx[i] = han_sinc(xoffset-(float)i);
533        hxtot += hx[i];
534     }
535     for (i=-barrier;i<=barrier;i++)
536     {
537        zpos = pos.el[2] + (float)i;
538        for (j=-barrier;j<=barrier;j++)
539        {
540           ypos = pos.el[1] + (float)j;
541           pfac = sinc_accum(rasptrs, tina_int(zpos),tina_int(ypos),
542                  tina_int(pos.el[0]-(float)barrier),&hx[-barrier]);
543           newpix += hz[i]*hy[j]*pfac;
544        }
545     }
546     return(newpix/(hxtot*hytot*hztot)); 
547 }
548 
549 
550 Imrect * seq_slicez(float zslice, Imregion *within, Vec3(*proj_func)( ))
551 {
552     Vec3 post;
553     int i,j;
554     float *row1;
555     Imrect *im;
556     Imregion roi;
557 
558     if (imptrs == NULL)
559     {
560        error("no data to display in seq_slicez() \n",warning);
561        return(NULL);
562     }
563     if (within != NULL)
564          roi = *within;
565     else
566     {
567        roi_fill(&roi,imptrlx,imptrly,imptrux,imptruy);
568     }
569 
570     im = im_alloc(roi.uy-roi.ly,roi.ux-roi.lx,&roi,float_v);
571 
572     seq_init_interp(imptrlx,imptrux,
573                 imptrly,imptruy,
574                 imptrlz,imptruz);
575     for (i=roi.ly;i<roi.uy;i++)
576     {
577        row1 = (float*) IM_ROW(im,i);
578        row1 = &row1[roi.lx];
579        for (j=roi.lx;j<roi.ux;j++)
580        {
581            post = proj_func((double)j+0.5,(double)i+0.5,(double)zslice);
582            *(row1++) = interp(imptrs,post);
583        }
584     }
585     return(im);
586 }
587 
588 Imrect * seq_slicey(float slicey, Imregion *within, Vec3(*proj_func)( ))
589 {
590     Vec3 post;
591     int j,k;
592     float *row1;
593     Imrect *im;
594     Imregion roi;
595 
596     if (imptrs == NULL)
597     {
598        error("no data to display in seq_slicey() \n",warning);
599        return(NULL);
600     }
601     if (within !=NULL)
602        roi = *within;
603     else
604     {
605        roi_fill(&roi,imptrlx,imptrlz,imptrux,imptruz);
606     }
607     im = im_alloc(roi.uy-roi.ly,roi.ux-roi.lx,&roi,float_v);
608     seq_init_interp(imptrlx,imptrux,
609                 imptrly,imptruy,
610                 imptrlz,imptruz);
611 
612     for(k=roi.ly;k<roi.uy;k++)
613     {
614         row1 = (float*) IM_ROW(im,k);
615         row1 = &row1[roi.lx];
616         for (j=roi.lx;j<roi.ux;j++)
617         {
618             post = proj_func((double)j+0.5,(double)slicey,(double)k+0.5);
619 
620             *(row1++) = interp(imptrs,post);
621         }
622     }
623     return(im);
624 }
625 
626 Imrect * seq_slicex(float slicex, Imregion *within, Vec3(*proj_func)( ))
627 {
628     Vec3 post;
629     int i,j;
630     float *row1;
631     Imrect *im;
632     Imregion roi;
633 
634     if (imptrs == NULL)
635     {
636        error("no data to display in seq_slicex() \n",warning);
637        return(NULL);
638     }
639     if (within != NULL)
640        roi = *within;
641     else
642     {
643        roi_fill(&roi,imptrlz,imptrly,imptruz,imptruy);
644     }
645     im = im_alloc(roi.uy-roi.ly,roi.ux-roi.lx,&roi,float_v);
646 
647     seq_init_interp(imptrlx,imptrux,
648                 imptrly,imptruy,
649                 imptrlz,imptruz);
650 
651     for (i=roi.ly;i<roi.uy;i++)
652     {
653        row1 = (float*) IM_ROW(im,i);
654        row1 = &row1[roi.lx];
655        for (j=roi.lx;j<roi.ux;j++)
656        {
657            post = proj_func((double)slicex,(double)i+0.5,(double)j+0.5);
658            *(row1++) = interp(imptrs,post);
659        }
660     }
661     return(im);
662 }
663 
664 static int get_total_frames(Sequence *seq)
665 {
666 
667 int          count;
668   List      *im_ptr;
669 
670   if (seq == NULL)
671     return(0);
672 
673   count = seq->offset;
674 
675   if ((im_ptr = get_seq_start_el(seq)) == NULL)
676     return(count);
677 
678   while(im_ptr->next != NULL)
679     {
680 
681       im_ptr = im_ptr->next;
682       count = count +1;
683     }
684 
685   return (count);
686 }
687 
688 /* get_end_frame finds the number of the end frame of a sequence, given the start image
689    in the sequence and the number of the start image. loops over the image sequence until
690    it reaches the end, counting as it goes. */
691 int get_end_frame(Sequence *seq)
692 {
693   int          count;
694   List      *im_ptr;
695 
696   if (seq == NULL)
697     return(0);
698 
699   count = seq->offset;
700 
701   if ((im_ptr = get_seq_start_el(seq)) == NULL)
702     return(count);
703  
704   while(im_ptr->next != NULL)
705     {
706  
707       im_ptr = im_ptr->next;
708       count = count + seq->stride;
709     }
710  
711   return (count);
712 }
713 
714 
715 void ***seq_slice_init(Sequence *seq)
716 {
717     List *frame;
718     Imrect *im, *imref;
719     Imregion roi;
720     int i,j;
721 
722     if (imptrs!=NULL) parray_free(imptrs,imptrlz,imptrly,imptruz,imptruy);
723     imptrs = NULL;
724     frame = get_seq_start_el(seq);
725     if((frame==NULL)||(imref = frame->to)==NULL) return(NULL);
726     seq_voxel_vtype(imref->vtype);
727     roi = *(imref->region);
728 
729     imptrlx = roi.lx;
730     imptrux = roi.ux;
731     imptrly = roi.ly;
732     imptruy = roi.uy;
733     imptrlz = seq->offset;
734     imptruz = get_total_frames(seq) + 1;
735     imptrs = (void ***)parray_alloc(imptrlz,imptrly,imptruz,imptruy);
736     im = imref;
737     for (i=imptrlz;i<imptruz;i++)
738     {
739        for (j=imptrly;j<imptruy;j++)
740        {
741            imptrs[i][j] = (void *)IM_ROW(im,j);
742        }
743        if( frame !=NULL)
744        {
745           frame = frame->next;
746           if (frame!=NULL) im = frame->to;
747        }
748     }
749 
750         return (imptrs);
751 }
752 
753 
754 void ***imptrs_swap(void ***imptrs_new)
755 {
756         void ***imptrs_temp=NULL;
757 
758         imptrs_temp = imptrs;
759         imptrs = imptrs_new;
760         return imptrs_temp;
761 }
762 
763 
764 static float nearest_pixel_mask(void ***rasptrs, Vec3 pos, int *mask_pix)
765 {
766     int x, y, z;
767     x = tina_int(pos.el[0]);
768     y = tina_int(pos.el[1]);
769     z = tina_int(pos.el[2]);
770     if (x<blx || x>=bux ||  y<bly || y>=buy ||  z<blz || z>=buz )
771     {
772        *mask_pix = 0;
773        return(outside);
774     }
775 
776     *mask_pix = 1;
777     return(seq_pixel(rasptrs,z, y, x));
778 }
779 
780 
781 static float tri_linear_mask(void ***rasptrs, Vec3 pos, int *mask_pix)
782 {
783     float xoffset, yoffset, zoffset;
784     float cenpix,newpix,nextpix;
785     float xfac, yfac, zfac;
786     int x , y, z;
787 
788     x = tina_int(pos.el[0]);
789     y = tina_int(pos.el[1]);
790     z = tina_int(pos.el[2]);
791     if (x<blx || x>=bux ||  y<bly || y>=buy || z<blz || z>=buz )
792     {
793        *mask_pix = 0;
794        return(outside);
795     }
796 
797     xoffset = pos.el[0] - x - 0.5;
798     yoffset = pos.el[1] - y - 0.5;
799     zoffset = pos.el[2] - z - 0.5;
800     cenpix = seq_pixel(rasptrs,z,y,x);
801     if (zoffset<0.0) 
802     {
803       nextpix = seq_pixel(rasptrs,z-1,y,x);
804       zfac = cenpix - nextpix;
805     }
806     else
807     {
808       nextpix = seq_pixel(rasptrs,z+1,y,x);
809       zfac = nextpix - cenpix;
810     }
811 
812     if (yoffset<0.0)
813     {
814       nextpix = seq_pixel(rasptrs,z,y-1,x); 
815       yfac = cenpix -nextpix;
816     }
817     else
818     {
819       nextpix = seq_pixel(rasptrs,z,y+1,x); 
820       yfac = nextpix - cenpix;
821     }
822 
823     if (xoffset<0.0)
824     {
825       nextpix = seq_pixel(rasptrs,z,y,x-1);
826       xfac = cenpix - nextpix;
827     }
828     else
829     {
830       nextpix = seq_pixel(rasptrs,z,y,x+1);
831       xfac = nextpix - cenpix;
832     }
833     newpix = cenpix + xfac*xoffset + yfac*yoffset + zfac*zoffset;
834 
835     *mask_pix = 1;
836     return(newpix);
837 }
838 
839 
840 static float sinc_interp_mask(void ***rasptrs, Vec3 pos, int *mask_pix)
841 {
842     float ypos, zpos;
843     float newpix=0.0;
844     float xoffset, yoffset, zoffset;
845     float hxtot=0.0, hytot=0.0, hztot=0.0;
846     float pfac;
847     static float *hx=NULL;                      /* static data! */
848     static float *hy=NULL, *hz =NULL;           /* static data! */
849     int i,j;
850     int x , y, z;
851 
852     x = tina_int(pos.el[0]);
853     y = tina_int(pos.el[1]);
854     z = tina_int(pos.el[2]);
855     if (x<blx || x>=bux ||  y<bly || y>=buy ||  z<blz || z>=buz )
856     {
857        *mask_pix = 0;
858        return(outside);
859     }
860 
861     xoffset = pos.el[0] - x - 0.5;
862     yoffset = pos.el[1] - y - 0.5;
863     zoffset = pos.el[2] - z - 0.5;
864     if (hx  == NULL)
865     {
866         hx = fvector_alloc(-6,7);
867         hy = fvector_alloc(-6,7);
868         hz = fvector_alloc(-6,7);
869     }
870     for (i=-barrier;i<=barrier;i++)
871     {
872        hz[i] = han_sinc(zoffset-(float)i);
873        hztot += hz[i];
874        hy[i] = han_sinc(yoffset-(float)i);
875        hytot += hy[i];
876        hx[i] = han_sinc(xoffset-(float)i);
877        hxtot += hx[i];
878     }
879     for (i=-barrier;i<=barrier;i++)
880     {
881        zpos = pos.el[2] + (float)i;
882        for (j=-barrier;j<=barrier;j++)
883        {
884           ypos = pos.el[1] + (float)j;
885           pfac = sinc_accum(rasptrs, tina_int(zpos),tina_int(ypos),
886                  tina_int(pos.el[0]-(float)barrier),&hx[-barrier]);
887           newpix += hz[i]*hy[j]*pfac;
888        }
889     }
890 
891     *mask_pix = 1;
892     return(newpix/(hxtot*hytot*hztot)); 
893 }
894 
895 
896 void **seq_slicez_mask(float zslice, Imregion *within, Vec3(*proj_func)( ))
897 {
898     Vec3 post;
899     int i,j, *rowm=NULL;
900     float *row1=NULL;
901     Imrect *im=NULL, *im_mask=NULL;
902     Imregion roi;
903     void **im_vector=NULL;
904 
905     if (imptrs == NULL)
906     {
907        error("no data to display in seq_slicez_mask() \n",warning);
908        return(NULL);
909     }
910     if (within != NULL)
911          roi = *within;
912     else
913     {
914        roi_fill(&roi,imptrlx,imptrly,imptrux,imptruy);
915     }
916 
917     im = im_alloc(roi.uy-roi.ly,roi.ux-roi.lx,&roi,float_v);
918     im_mask = im_alloc(roi.uy-roi.ly,roi.ux-roi.lx,&roi,int_v);
919     im_vector = pvector_alloc(0, 2);
920 
921     seq_init_interp(imptrlx,imptrux,
922                 imptrly,imptruy,
923                 imptrlz,imptruz);
924     for (i=roi.ly;i<roi.uy;i++)
925     {
926        row1 = (float*) IM_ROW(im,i);
927        row1 = &row1[roi.lx];
928        rowm = (int*) IM_ROW(im_mask,i);
929        rowm = &rowm[roi.lx];
930 
931        for (j=roi.lx;j<roi.ux;j++)
932        {
933            post = proj_func((double)j+0.5,(double)i+0.5,(double)zslice);
934            *(row1++) = interp_mask(imptrs, post, rowm++);
935        }
936     }
937 
938     im_vector[0] = im;
939     im_vector[1] = im_mask;
940     return im_vector;
941 }
942 
943 
944 void **seq_slicey_mask(float slicey, Imregion *within, Vec3(*proj_func)( ))
945 {
946     Vec3 post;
947     int j,k, *rowm=NULL;
948     float *row1=NULL;
949     Imrect *im=NULL, *im_mask=NULL;
950     Imregion roi;
951     void **im_vector=NULL;
952 
953     if (imptrs == NULL)
954     {
955        error("no data to display in seq_slicey_mask() \n",warning);
956        return(NULL);
957     }
958     if (within !=NULL)
959        roi = *within;
960     else
961     {
962        roi_fill(&roi,imptrlx,imptrlz,imptrux,imptruz);
963     }
964     im = im_alloc(roi.uy-roi.ly,roi.ux-roi.lx,&roi,float_v);
965     im_mask = im_alloc(roi.uy-roi.ly,roi.ux-roi.lx,&roi,int_v);
966     im_vector = pvector_alloc(0, 2);
967 
968     seq_init_interp(imptrlx,imptrux,
969                 imptrly,imptruy,
970                 imptrlz,imptruz);
971 
972     for(k=roi.ly;k<roi.uy;k++)
973     {
974         row1 = (float*) IM_ROW(im,k);
975         row1 = &row1[roi.lx];
976         rowm = (int*) IM_ROW(im_mask,k);
977         rowm = &rowm[roi.lx];
978 
979         for (j=roi.lx;j<roi.ux;j++)
980         {
981             post = proj_func((double)j+0.5,(double)slicey,(double)k+0.5);
982             *(row1++) = interp_mask(imptrs, post, rowm++);
983         }
984     }
985 
986     im_vector[0] = im;
987     im_vector[1] = im_mask;
988     return im_vector;
989 }
990 
991 
992 void **seq_slicex_mask(float slicex, Imregion *within, Vec3(*proj_func)( ))
993 {
994     Vec3 post;
995     int i, j, *rowm=NULL;
996     float *row1=NULL;
997     Imrect *im=NULL, *im_mask=NULL;
998     Imregion roi;
999     void **im_vector=NULL;
1000 
1001     if (imptrs == NULL)
1002     {
1003        error("no data to display in seq_slicex_mask() \n",warning);
1004        return(NULL);
1005     }
1006     if (within != NULL)
1007        roi = *within;
1008     else
1009     {
1010        roi_fill(&roi,imptrlz,imptrly,imptruz,imptruy);
1011     }
1012     im = im_alloc(roi.uy-roi.ly,roi.ux-roi.lx,&roi,float_v);
1013     im_mask = im_alloc(roi.uy-roi.ly,roi.ux-roi.lx,&roi,int_v);
1014     im_vector = pvector_alloc(0, 2);
1015 
1016     seq_init_interp(imptrlx,imptrux,
1017                 imptrly,imptruy,
1018                 imptrlz,imptruz);
1019 
1020     for (i=roi.ly;i<roi.uy;i++)
1021     {
1022        row1 = (float*) IM_ROW(im,i);
1023        row1 = &row1[roi.lx];
1024        rowm = (int*) IM_ROW(im_mask,i);
1025        rowm = &rowm[roi.lx];
1026 
1027        for (j=roi.lx;j<roi.ux;j++)
1028        {
1029            post = proj_func((double)slicex,(double)i+0.5,(double)j+0.5);
1030            *(row1++) = interp_mask(imptrs, post, rowm++);
1031        }
1032     }
1033 
1034     im_vector[0] = im;
1035     im_vector[1] = im_mask;
1036     return im_vector;
1037 }
1038 
1039 
1040 /* mjs 10.10.03 changing seq_interp_choice to return the old value of the barrier, in
1041    order that the barrier can be set and the reset back to original value */
1042 /* Not sure whether this really needs to be done, if problematic please change back */
1043 
1044 int seq_interp_choice(int choice)
1045 {
1046   int old_val = barrier;
1047 
1048    if (choice == 0)
1049    {
1050       interp = nearest_pixel;
1051       interp_mask = nearest_pixel_mask;
1052       barrier = 0;
1053    }
1054    if (choice == 1)
1055    {
1056       interp = tri_linear;
1057       interp_mask = tri_linear_mask;
1058       barrier = 1;
1059    }
1060    if (choice ==2)
1061    {
1062       interp = sinc_interp;
1063       interp_mask = sinc_interp_mask;
1064       barrier = 2;
1065    }
1066    if (choice ==3)
1067    {
1068       interp = sinc_interp;
1069       interp_mask = sinc_interp_mask;
1070       barrier = 3;
1071    }
1072    if (choice ==4)
1073    {
1074       interp = sinc_interp;
1075       interp_mask = sinc_interp_mask;
1076       barrier = 4;
1077    }
1078    if (choice ==5)
1079    {
1080       interp = sinc_interp;
1081       interp_mask = sinc_interp_mask;
1082       barrier = 5;
1083    }
1084    if (choice ==6)
1085    {
1086       interp = sinc_interp;
1087       interp_mask = sinc_interp_mask;
1088       barrier = 6;
1089    }
1090    return(old_val);
1091 }

~ [ 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.