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

Linux Cross Reference
Tina4/src/tools/coreg/coreg_slices.c

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

  1 #include <stdio.h>
  2 #include <tina/sys.h>
  3 #include <tina/sysfuncs.h>
  4 #include <tina/math.h>
  5 #include <tina/mathfuncs.h>
  6 #include <tina/vision.h>
  7 extern void ***coreg_limits(int *lz, int *uz, int *ly,
  8                          int *uy, int *lx, int *ux);
  9 static void ***imptrs=NULL;
 10 int imptrlx, imptrux,imptrly, imptruy, imptrlz, imptruz;
 11 static int blx, bux, bly, buy, blz, buz;
 12 static float (*interp)(void ***rasptrs, Vec3 pos) = NULL; 
 13 static float (*coreg_pixel)(void ***rasptrs, int z, int y, int x) = NULL;
 14 static double (*sinc_accum)(void ***rasptrs, int z, int y, int x, float *hxp) = NULL;
 15 static int barrier=0;
 16 static Vartype inter_vtype = short_v;
 17 
 18 
 19 static float coreg_cpixel(void ***rasptrs, int z, int y, int x)
 20 {
 21     char *row;
 22     row = rasptrs[z][y];
 23     return((float)row[x]);
 24 }
 25 
 26 static float coreg_ucpixel(void ***rasptrs, int z, int y, int x)
 27 {
 28     unsigned char *row;
 29     row = rasptrs[z][y];
 30     return((float)row[x]);
 31 }
 32 
 33 static float coreg_spixel(void ***rasptrs, int z, int y, int x)
 34 {
 35     short *row;
 36     row = rasptrs[z][y];
 37     return((float)row[x]);
 38 }
 39 
 40 static float coreg_uspixel(void ***rasptrs, int z, int y, int x)
 41 {
 42     unsigned short *row;
 43     row = rasptrs[z][y];
 44     return((float)row[x]);
 45 }
 46 
 47 static float coreg_ipixel(void ***rasptrs, int z, int y, int x)
 48 {
 49     int *row;
 50     row = rasptrs[z][y];
 51     return((float)row[x]);
 52 }
 53 
 54 static float coreg_fpixel(void ***rasptrs,int z, int y, int x)
 55 {
 56     float *row;
 57     row = rasptrs[z][y];
 58     return((float)row[x]);
 59 }
 60 
 61 double sinc_accumuc(void *** rasptrs, int z, int y, int x, float *hxp)
 62 {
 63      unsigned char *row;
 64      unsigned char *thispix;
 65      double pfac;
 66  
 67      row = rasptrs[z][y];
 68      thispix = &row[x];
 69      pfac = 0.0;
 70      switch(barrier)
 71      {
 72         case 6:
 73           pfac += *(thispix++)* *(hxp++);
 74           pfac += *(thispix++)* *(hxp++);
 75             
 76         case 5:
 77           pfac += *(thispix++)* *(hxp++);
 78           pfac += *(thispix++)* *(hxp++);
 79             
 80         case 4:
 81           pfac += *(thispix++)* *(hxp++);
 82           pfac += *(thispix++)* *(hxp++);
 83             
 84         case 3:
 85           pfac += *(thispix++)* *(hxp++);
 86           pfac += *(thispix++)* *(hxp++);
 87             
 88         case 2:
 89           pfac += *(thispix++)* *(hxp++);
 90           pfac += *(thispix++)* *(hxp++);
 91             
 92        case 1:
 93          pfac += *(thispix++)* *(hxp++);
 94          pfac += *(thispix++)* *(hxp++);
 95          pfac += *(thispix++)* *(hxp++);
 96     }       
 97     return (pfac);
 98 }
 99 
100 double sinc_accumc(void *** rasptrs, int z, int y, int x, float *hxp)
101 {
102      char *row;
103      char *thispix;
104      double pfac;
105 
106      row = rasptrs[z][y];
107      thispix = &row[x];
108      pfac = 0.0;
109      switch(barrier)
110      {
111         case 6:
112           pfac += *(thispix++)* *(hxp++);
113           pfac += *(thispix++)* *(hxp++);
114 
115         case 5:
116           pfac += *(thispix++)* *(hxp++);
117           pfac += *(thispix++)* *(hxp++);
118 
119         case 4:
120           pfac += *(thispix++)* *(hxp++);
121           pfac += *(thispix++)* *(hxp++);
122 
123         case 3:
124           pfac += *(thispix++)* *(hxp++);
125           pfac += *(thispix++)* *(hxp++);
126 
127         case 2:
128           pfac += *(thispix++)* *(hxp++);
129           pfac += *(thispix++)* *(hxp++);
130 
131        case 1:
132          pfac += *(thispix++)* *(hxp++);
133          pfac += *(thispix++)* *(hxp++);
134          pfac += *(thispix++)* *(hxp++);
135     }
136     return (pfac);
137 }
138 
139 double sinc_accums(void *** rasptrs, int z, int y, int x, float *hxp)
140 {
141      short *row;
142      short *thispix;
143      double pfac;
144 
145      row = rasptrs[z][y];
146      thispix = &row[x];
147      pfac = 0.0;
148      switch(barrier)
149      {
150         case 6:
151           pfac += *(thispix++)* *(hxp++);
152           pfac += *(thispix++)* *(hxp++);
153 
154         case 5:
155           pfac += *(thispix++)* *(hxp++);
156           pfac += *(thispix++)* *(hxp++);
157 
158         case 4:
159           pfac += *(thispix++)* *(hxp++);
160           pfac += *(thispix++)* *(hxp++);
161 
162         case 3:
163           pfac += *(thispix++)* *(hxp++);
164           pfac += *(thispix++)* *(hxp++);
165 
166         case 2:
167           pfac += *(thispix++)* *(hxp++);
168           pfac += *(thispix++)* *(hxp++);
169 
170        case 1:
171          pfac += *(thispix++)* *(hxp++);
172          pfac += *(thispix++)* *(hxp++);
173          pfac += *(thispix++)* *(hxp++);
174     }
175     return (pfac);
176 }
177 
178 double sinc_accumus(void *** rasptrs, int z, int y, int x, float *hxp)
179 {
180      unsigned short *row;
181      unsigned short *thispix;
182      double pfac;
183 
184      row = rasptrs[z][y];
185      thispix = &row[x];
186      pfac = 0.0;
187      switch(barrier)
188      {
189         case 6:
190           pfac += *(thispix++)* *(hxp++);
191           pfac += *(thispix++)* *(hxp++);
192 
193         case 5:
194           pfac += *(thispix++)* *(hxp++);
195           pfac += *(thispix++)* *(hxp++);
196 
197         case 4:
198           pfac += *(thispix++)* *(hxp++);
199           pfac += *(thispix++)* *(hxp++);
200 
201         case 3:
202           pfac += *(thispix++)* *(hxp++);
203           pfac += *(thispix++)* *(hxp++);
204 
205         case 2:
206           pfac += *(thispix++)* *(hxp++);
207           pfac += *(thispix++)* *(hxp++);
208 
209        case 1:
210          pfac += *(thispix++)* *(hxp++);
211          pfac += *(thispix++)* *(hxp++);
212          pfac += *(thispix++)* *(hxp++);
213     }
214     return (pfac);
215 }
216 
217 double sinc_accumi(void *** rasptrs, int z, int y, int x, float *hxp)
218 {
219      int *row;
220      int *thispix;
221      double pfac;
222 
223      row = rasptrs[z][y];
224      thispix = &row[x];
225      pfac = 0.0;
226      switch(barrier)
227      {
228         case 6:
229           pfac += *(thispix++)* *(hxp++);
230           pfac += *(thispix++)* *(hxp++);
231 
232         case 5:
233           pfac += *(thispix++)* *(hxp++);
234           pfac += *(thispix++)* *(hxp++);
235 
236         case 4:
237           pfac += *(thispix++)* *(hxp++);
238           pfac += *(thispix++)* *(hxp++);
239 
240         case 3:
241           pfac += *(thispix++)* *(hxp++);
242           pfac += *(thispix++)* *(hxp++);
243 
244         case 2:
245           pfac += *(thispix++)* *(hxp++);
246           pfac += *(thispix++)* *(hxp++);
247 
248        case 1:
249          pfac += *(thispix++)* *(hxp++);
250          pfac += *(thispix++)* *(hxp++);
251          pfac += *(thispix++)* *(hxp++);
252     }
253     return (pfac);
254 }
255 
256 double sinc_accumf(void *** rasptrs, int z, int y, int x, float *hxp)
257 {
258      float *row;
259      float *thispix;
260      double pfac;
261 
262      row = rasptrs[z][y];
263      thispix = &row[x];
264      pfac = 0.0;
265      switch(barrier)
266      {
267         case 6:
268           pfac += *(thispix++)* *(hxp++);
269           pfac += *(thispix++)* *(hxp++);
270 
271         case 5:
272           pfac += *(thispix++)* *(hxp++);
273           pfac += *(thispix++)* *(hxp++);
274 
275         case 4:
276           pfac += *(thispix++)* *(hxp++);
277           pfac += *(thispix++)* *(hxp++);
278 
279         case 3:
280           pfac += *(thispix++)* *(hxp++);
281           pfac += *(thispix++)* *(hxp++);
282 
283         case 2:
284           pfac += *(thispix++)* *(hxp++);
285           pfac += *(thispix++)* *(hxp++);
286 
287        case 1:
288          pfac += *(thispix++)* *(hxp++);
289          pfac += *(thispix++)* *(hxp++);
290          pfac += *(thispix++)* *(hxp++);
291     }
292     return (pfac);
293 }
294 
295 
296 void set_voxel_vtype(Vartype vtype)
297 {
298     if (vtype == char_v) coreg_pixel = coreg_cpixel; 
299     if (vtype == char_v) sinc_accum = sinc_accumc; 
300     if (vtype == uchar_v) coreg_pixel = coreg_ucpixel; 
301     if (vtype == uchar_v) sinc_accum = sinc_accumuc; 
302     if (vtype == short_v) coreg_pixel = coreg_spixel; 
303     if (vtype == short_v) sinc_accum = sinc_accums; 
304     if (vtype == ushort_v) coreg_pixel = coreg_uspixel; 
305     if (vtype == ushort_v) sinc_accum = sinc_accumus; 
306     if (vtype == int_v) coreg_pixel = coreg_ipixel; 
307     if (vtype == int_v) sinc_accum = sinc_accumi; 
308     if (vtype == float_v) coreg_pixel = coreg_fpixel; 
309     if (vtype == float_v) sinc_accum = sinc_accumf; 
310 }
311 
312 void init_interp(int nblx, int nbux, int nbly, int nbuy, int nblz, int nbuz)
313 {
314    blx = nblx;
315    bux = nbux;
316    bly = nbly;
317    buy = nbuy;
318    blz = nblz;
319    buz = nbuz;
320 }
321 
322 float this_pixel(void ***rasptrs, int z, int y, int x)
323 {
324    return(coreg_pixel(rasptrs,z, y, x));
325 }
326 
327 float nearest_pixel(void ***rasptrs, Vec3 pos)
328 {
329     int x, y, z;
330     x = tina_int(pos.el[0]);
331     y = tina_int(pos.el[1]);
332     z = tina_int(pos.el[2]);
333     if (x<blx || x>=bux 
334     ||  y<bly || y>=buy
335     ||  z<blz || z>=buz ) 
336        return(0.0);
337     return(coreg_pixel(rasptrs,z, y, x));
338 }
339 
340 float bi_linear(void ***rasptrs, Vec3 pos)
341 {
342     float xoffset, yoffset, zoffset;
343     float cenpix,newpix,nextpix;
344     float xfac, yfac, zfac;
345     int x , y, z;
346 
347     x = tina_int(pos.el[0]);
348     y = tina_int(pos.el[1]);
349     z = tina_int(pos.el[2]);
350     if (x<blx || x>=bux 
351     ||  y<bly || y>=buy
352     ||  z<blz || z>=buz )        
353        return(0.0);
354 
355     xoffset = pos.el[0] - x - 0.5;
356     yoffset = pos.el[1] - y - 0.5;
357     zoffset = pos.el[2] - z - 0.5;
358     cenpix = coreg_pixel(rasptrs,z,y,x);
359 /*    if (cenpix == 0.0) return(0.0);
360 */
361     if (zoffset<0.0) 
362     {
363       nextpix = coreg_pixel(rasptrs,z-1,y,x);
364 /*      if (nextpix ==0.0) return(0.0);
365 */
366       zfac = cenpix - nextpix;
367     }
368     else
369     {
370       nextpix = coreg_pixel(rasptrs,z+1,y,x);
371 /*      if (nextpix ==0.0) return(0.0);
372 */
373       zfac = nextpix - cenpix;
374     }
375 
376     if (yoffset<0.0)
377     {
378       nextpix = coreg_pixel(rasptrs,z,y-1,x); 
379 /*      if (nextpix ==0.0) return(0.0);
380 */
381       yfac = cenpix -nextpix;
382     }
383     else
384     {
385       nextpix = coreg_pixel(rasptrs,z,y+1,x); 
386 /*      if (nextpix ==0.0) return(0.0);
387 */
388       yfac = nextpix - cenpix;
389     }
390 
391     if (xoffset<0.0)
392     {
393       nextpix = coreg_pixel(rasptrs,z,y,x-1);
394 /*      if (nextpix ==0.0) return(0.0);
395 */
396       xfac = cenpix - nextpix;
397     }
398     else
399     {
400       nextpix = coreg_pixel(rasptrs,z,y,x+1);
401 /*      if (nextpix ==0.0) return(0.0);
402 */
403       xfac = nextpix - cenpix;
404     }
405     newpix = cenpix + xfac*xoffset + yfac*yoffset + zfac*zoffset;
406 
407     return(newpix);
408 }
409 
410 float han_sinc(float delta)
411 {
412    double sinc ,phase,hanning;
413 
414    phase = PI*delta;
415    if(fabs(phase) > 0.000001) sinc = sin(phase)/phase;
416    else sinc = 1.0;
417    hanning = (0.5 + 0.5*cos(phase/(float)(1.0+barrier)));
418    return(hanning*sinc);
419 }
420 
421 float sinc_interp(void ***rasptrs, Vec3 pos)
422 {
423     float xpos, ypos, zpos;
424     float newpix=0.0;
425     float xoffset, yoffset, zoffset;
426     float hxtot=0.0, hytot=0.0, hztot=0.0;
427     float pfac;
428     static float *hx=NULL;
429     static float *hy=NULL, *hz =NULL;
430     float *hxp;
431     int i,j,k;
432     int x , y, z;
433 
434     x = tina_int(pos.el[0]);
435     y = tina_int(pos.el[1]);
436     z = tina_int(pos.el[2]);
437     if (x<blx || x>=bux 
438     ||  y<bly || y>=buy
439     ||  z<blz || z>=buz )        
440        return(0.0);
441    
442     xoffset = pos.el[0] - x - 0.5;
443     yoffset = pos.el[1] - y - 0.5;
444     zoffset = pos.el[2] - z - 0.5;
445     if (hx  == NULL)
446     {
447         hx = fvector_alloc(-6,7);
448         hy = fvector_alloc(-6,7);
449         hz = fvector_alloc(-6,7);
450     }
451     for (i=-barrier;i<=barrier;i++)
452     {
453        hz[i] = han_sinc(zoffset-(float)i);
454        hztot += hz[i];
455        hy[i] = han_sinc(yoffset-(float)i);
456        hytot += hy[i];
457        hx[i] = han_sinc(xoffset-(float)i);
458        hxtot += hx[i];
459     }
460     for (i=-barrier;i<=barrier;i++)
461     {
462        zpos = pos.el[2] + (float)i;
463        for (j=-barrier;j<=barrier;j++)
464        {
465           ypos = pos.el[1] + (float)j;
466           pfac = sinc_accum(rasptrs, tina_int(zpos),tina_int(ypos),
467                  tina_int(pos.el[0]-(float)barrier),&hx[-barrier]);
468           newpix += hz[i]*hy[j]*pfac;
469        }
470     }
471     return(newpix/(hxtot*hytot*hztot)); 
472 }
473 
474 void set_interp_choice(int choice)
475 {
476    if (choice == 0)
477    {
478       interp = nearest_pixel;
479       barrier = 0;
480    }
481    if (choice == 1)
482    {
483       interp = bi_linear;
484       barrier = 1;
485    }
486    if (choice ==2)
487    {
488       interp = sinc_interp;
489       barrier = 2;
490    }
491    if (choice ==3)
492    {
493       interp = sinc_interp;
494       barrier = 3;
495    }
496    if (choice ==4)
497    {
498       interp = sinc_interp;
499       barrier = 4;
500    }
501    if (choice ==5)
502    {
503       interp = sinc_interp;
504       barrier = 5;
505    }
506    if (choice ==6)
507    {
508       interp = sinc_interp;
509       barrier = 6;
510    }
511 }
512 
513 Imrect * coreg_slicez(float zslice, Imregion *within)
514 {
515     Vec3 post;
516     Vec3 coreg_bproj(double x, double y, double z);
517     int i,j,k;
518     float *row1;
519     short *row;
520     Imrect *im;
521     Imregion roi;
522 
523     imptrs = coreg_limits(&imptrlz,&imptruz,&imptrly,
524                           &imptruy,&imptrlx,&imptrux);
525     if (imptrs == NULL)
526     {
527        error("no data to display in coreg_slicez() \n",warning);
528        return(NULL);
529     }
530     if (within != NULL)
531          roi = *within;
532     else
533     {
534        roi_fill(&roi,imptrlx,imptrly,imptrux,imptruy);
535     }
536 
537     im = im_alloc(roi.uy-roi.ly,roi.ux-roi.lx,&roi,float_v);
538 
539     init_interp(imptrlx+barrier,imptrux-barrier,
540                 imptrly+barrier,imptruy-barrier,
541                 imptrlz+barrier,imptruz-barrier);
542     for (i=roi.ly;i<roi.uy;i++)
543     {
544        row1 = (float*) IM_ROW(im,i);
545        row1 = &row1[roi.lx];
546        for (j=roi.lx;j<roi.ux;j++)
547        {
548            post = coreg_bproj((double)j+0.5,(double)i+0.5,(double)zslice);
549            *(row1++) = interp(imptrs,post); 
550        }
551     }
552     return(im);
553 }
554 
555 Imrect * coreg_slicey(float slicey, Imregion *within)
556 {
557     Vec3 coreg_bproj(double x, double y, double z);    
558     Vec3 post;
559     int i,j,k;
560     float *row1;
561     short *row;
562     Imrect *im;
563     Imregion roi;
564 
565     imptrs = coreg_limits(&imptrlz,&imptruz,&imptrly,
566                           &imptruy,&imptrlx,&imptrux);
567     if (imptrs == NULL)
568     {
569        error("no data to display in coreg_slicey() \n",warning);
570        return(NULL);
571     }
572     if (within !=NULL)
573        roi = *within;
574     else
575     {
576        roi_fill(&roi,imptrlx,imptrlz,imptrux,imptruz);
577     }
578     im = im_alloc(roi.uy-roi.ly,roi.ux-roi.lx,&roi,float_v);
579     init_interp(imptrlx+barrier,imptrux-barrier,
580                 imptrly+barrier,imptruy-barrier,
581                 imptrlz+barrier,imptruz-barrier);
582 
583     for(k=roi.ly;k<roi.uy;k++)
584     {
585         row1 = (float*) IM_ROW(im,k);
586         row1 = &row1[roi.lx];
587         for (j=roi.lx;j<roi.ux;j++)
588         {
589             post = coreg_bproj((double)j+0.5,(double)slicey,(double)k+0.5);
590 
591             *(row1++) = interp(imptrs,post); 
592         }
593     }
594     return(im);
595 }
596 
597 Imrect * coreg_slicex(float slicex, Imregion *within)
598 {
599     Vec3 coreg_bproj(double x, double y, double z);   
600     Vec3 post;
601     int i,j,k;
602     float *row1;
603     short *row;
604     Imrect *im;
605     Imregion roi;
606 
607 
608     imptrs = coreg_limits(&imptrlz,&imptruz,&imptrly,
609                           &imptruy,&imptrlx,&imptrux);
610     if (imptrs == NULL)
611     {
612        error("no data to display in coreg_slicex() \n",warning);
613        return(NULL);
614     }
615     if (within != NULL)
616        roi = *within;
617     else
618     {
619        roi_fill(&roi,imptrlz,imptrly,imptruz,imptruy);
620     }
621     im = im_alloc(roi.uy-roi.ly,roi.ux-roi.lx,&roi,float_v);
622 
623     init_interp(imptrlx+barrier,imptrux-barrier,
624                 imptrly+barrier,imptruy-barrier,
625                 imptrlz+barrier,imptruz-barrier);
626 
627     for (i=roi.ly;i<roi.uy;i++)
628     {
629        row1 = (float*) IM_ROW(im,i);
630        row1 = &row1[roi.lx];
631        for (j=roi.lx;j<roi.ux;j++)
632        {
633            post = coreg_bproj((double)slicex,(double)i+0.5,(double)j+0.5);
634            *(row1++) = interp(imptrs,post);
635        }
636     }
637     return(im);
638 }
639 
640 

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