TPCCLIB
Loading...
Searching...
No Matches
Functions
fbp.c File Reference

Image reconstruction using filtered back-projection. More...

#include "libtpcrec.h"

Go to the source code of this file.

Functions

int fbp (float *sinogram, int rays, int views, int dim, int filter, float cutoff, float *image)
int imgFBP (IMG *scn, IMG *img, int imgDim, float zoom, int filter, float cutoff, float shiftX, float shiftY, float rotation, int verbose)
int fbpMakeFilter (float *sinFFT, float *cosFFT, float cutoff, int filter, int lenFFT, int views, float *fbpFilter, int verbose)
void fbp_fft_bidir_complex_radix2 (float *real, float *imag, int direct, int n, float *sine, float *cosine)
void fbp_back_proj (float *prj, float *imgCorner, int imgDim, int view, int views, int rays, float offsX, float offsY, float *sinB, float *sinBrot)
void fbp_back_proj_round (float *prj, float *imgOrigin, int imgDim, int view, int views, int rays, float offsX, float offsY, float bpZoom, float *sinB, float *sinBrot)

Detailed Description

Image reconstruction using filtered back-projection.

Based on the program fbprec (May 1997) written by Sakari Alenius for Sun UNIX workstations.

Author
Vesa Oikonen

Definition in file fbp.c.

Function Documentation

◆ fbp()

int fbp ( float * sinogram,
int rays,
int views,
int dim,
int filter,
float cutoff,
float * image )

Filtered back-projection (FBP) reconstruction of one 2D data matrix given as an array of floats.

Image cannot be zoomed (zoom=1), rotated, or shifted.

See also
imgFBP, reprojection
Returns
Returns 0 if ok.
Parameters
sinogramPointer to float array containing rays*views sinogram values.
raysNr of rays (bins or columns) in sinogram data.
viewsNr of views (rows) in sinogram data.
dimImage x and y dimensions; must be an even number, preferrably the same as number of rays, but there is no reason for it to be any larger than that.
filterFilter: 0=None, 1=Ramp, 2=Butter, 3=Hann, 4=Hamm, 5=Parzen, 6=Shepp.
cutoffNoise cut-off (for example 0.3).
imagePointer to pre-allocated image data; size must be at least dim*dim.

Definition at line 20 of file fbp.c.

37 {
38 if(sinogram==NULL || rays<2 || views<2 || dim<2 || image==NULL) return(1);
39 if(dim%2) return(2);
40 if(cutoff<0.05 || cutoff>0.99) return(2);
41
42 int i, n;
43
44 /* Set FFT array length */
45 int lenFFT=(2*rays)-1;
46 if(lenFFT<256) lenFFT=256; else if(lenFFT<512) lenFFT=512;
47 else if(lenFFT<1024) lenFFT=1024; else if(lenFFT<2048) lenFFT=2048;
48 else if(lenFFT<4096) lenFFT=4096; else if(lenFFT<8192) lenFFT=8192;
49 else if(lenFFT<16384) lenFFT=16384; else return(3);
50
51 /* Set image pixel values to zero */
52 for(i=0; i<dim*dim; i++) image[i]=0.0;
53
54 /* Allocate memory and set data pointers */
55 float *locData=
56 calloc((2*lenFFT + (lenFFT+1) + 5*lenFFT/4 + 3*views/2), sizeof(float));
57 if(locData==NULL) return(4);
58 float *real, *imag;
59 float *sinFFT, *cosFFT, *fbpFilter, *sinBP;
60 real=locData; imag=real+lenFFT;
61 fbpFilter=imag+lenFFT;
62 sinFFT=fbpFilter+(lenFFT+1);
63 cosFFT=sinFFT+lenFFT/4; // sinFFT and cosFFT overlap on purpose
64 sinBP=sinFFT+5*lenFFT/4;
65
66 /* Pre-compute the sine and cosine tables for FFT */
67 {
68 double df, dg;
69 dg=(double)M_PI*2.0/(double)lenFFT;
70 for(i=0, df=0.0; i<5*lenFFT/4; i++, df+=dg) sinFFT[i]=(float)sin(df);
71 }
72 /* Compute reconstruction filter coefficients */
73 n=fbpMakeFilter(sinFFT, cosFFT, cutoff, filter, lenFFT, views, fbpFilter, 0);
74 if(n) {free(locData); return(100+n);}
75 /* Compute the sine tables for back-projection */
76 recSinTables(views, sinBP, NULL, 0.0);
77
78 /* Process each sinogram row (projection, view), two at the same time */
79 for(int pi=0; pi<views; pi+=2) {
80 /* Copy the data to real and imag arrays */
81 for(i=0; i<rays; i++) real[i]=sinogram[pi*rays+i];
82 for(i=0; i<rays; i++) imag[i]=sinogram[pi*rays+rays+i];
83 for(i=rays; i<lenFFT; i++) real[i]=imag[i]=0.0; // zero padding
84 /* Forward FFT */
85 fbp_fft_bidir_complex_radix2(real, imag, 1, lenFFT, sinFFT, cosFFT);
86 /* Filter */
87 for(i=0; i<lenFFT; i++) real[i]*=fbpFilter[i];
88 for(i=0; i<lenFFT; i++) imag[i]*=fbpFilter[i];
89 /* Inverse FFT */
90 fbp_fft_bidir_complex_radix2(real, imag, -1, lenFFT, sinFFT, cosFFT);
91 /* Back-projection */
92 fbp_back_proj_round(real, image+dim*(dim/2-1)+dim/2, dim, pi, views, rays,
93 0.0, 0.0, 1.0, sinBP, sinBP);
94 fbp_back_proj_round(imag, image+dim*(dim/2-1)+dim/2, dim, pi+1, views, rays,
95 0.0, 0.0, 1.0, sinBP, sinBP);
96 } // next two projections
97
98 free(locData);
99 return(0);
100}
fbpMakeFilter
int fbpMakeFilter(float *sinFFT, float *cosFFT, float cutoff, int filter, int lenFFT, int views, float *fbpFilter, int verbose)
Definition fbp.c:382
fbp_fft_bidir_complex_radix2
void fbp_fft_bidir_complex_radix2(float *real, float *imag, int direct, int n, float *sine, float *cosine)
Definition fbp.c:512
fbp_back_proj_round
void fbp_back_proj_round(float *prj, float *imgOrigin, int imgDim, int view, int views, int rays, float offsX, float offsY, float bpZoom, float *sinB, float *sinBrot)
Definition fbp.c:655
recSinTables
void recSinTables(int views, float *sinB, float *sinBrot, float rotation)
Definition recutil.c:12

◆ fbp_back_proj()

void fbp_back_proj ( float * prj,
float * imgCorner,
int imgDim,
int view,
int views,
int rays,
float offsX,
float offsY,
float * sinB,
float * sinBrot )

Backprojection with all projection profiles overlapping the image area.

Parameters
prjpointer to projection data
imgCornerpointer to destination image (upper left corner)
imgDimimage dimension
viewprojection view number
viewsnumber of projection views
raysnumber of rays
offsXx offset
offsYy offset
sinBSine tables
sinBrotSinBrot array

Definition at line 599 of file fbp.c.

620 {
621 float si, co, tleft, sir, cor, rayCenter, dif;
622 float t; /* distance between projection ray and origin */
623 int x, y, halfDim, leftX, rightTopLimit, yBottom, ti;
624
625 float *imgp=imgCorner;
626
627 halfDim=imgDim/2;
628 leftX=-halfDim;
629 rightTopLimit=halfDim-1;
630 yBottom=-halfDim;
631 rayCenter=0.5*(float)rays;
632 si=sinB[view]; /* Sine and cosine for x,y-shift */
633 co=sinB[view+views/2];
634 sir=sinBrot[view]; /* Sine and cosine for backprojection */
635 cor=sinBrot[view+views/2];
636 x=leftX; y=rightTopLimit;
637 tleft= rayCenter - (float)y*sir - offsY*si + ((float)(x+1))*cor + offsX*co;
638 for(; y>=yBottom; y--) {
639 t=tleft;
640 for(x=leftX; x<=rightTopLimit; x++, t+=cor) {
641 ti=(int)t; dif=t-(float)ti;
642 //ti=(int)(t+0.5); *imgp++ += prj[ti];
643 *imgp++ += prj[ti] + (prj[ti+1]-prj[ti])*dif;
644 }
645 tleft+=sir;
646 }
647}

Referenced by imgFBP().

◆ fbp_back_proj_round()

void fbp_back_proj_round ( float * prj,
float * imgOrigin,
int imgDim,
int view,
int views,
int rays,
float offsX,
float offsY,
float bpZoom,
float * sinB,
float * sinBrot )

Backprojection with projection profiles NOT overlapping the image area.

The corner pixels are left unchanged.

Parameters
prjpointer to projection data
imgOriginpointer to image origin (center)
imgDimimage dimension
viewprojection view number
viewsnumber of projection views
raysnumber of rays
offsXx offset in pixels
offsYy offset in pixels
bpZoomZoom
sinBSine tables
sinBrotsinBrot array

Definition at line 655 of file fbp.c.

678 {
679 if(0) printf("fbp_back_proj_round(..., %d, %d, %d, %d, %g, %g, %g, ...)\n",
680 imgDim, view, views, rays, offsX, offsY, bpZoom);
681
682 float *imgp, si, co, sir, cor, toffs, rayCenter, tPow2, yph, dif;
683 float t; /* distance between projection ray and origo */
684 int x, y, halfDim, right_top_limit, xrightround, ybottomround, ti;
685
686 /* Sine and cosine for x,y-shift */
687 si=sinB[view]*offsY;
688 co=sinB[views/2+view]*offsX;
689 /* Sine and cosine for backprojection */
690 sir=sinBrot[view];
691 cor=sinBrot[views/2+view];
692 halfDim=imgDim/2;
693 y=halfDim-2;
694 rayCenter=0.5*(float)rays;
695 if((float)y>(rayCenter*bpZoom)) {
696 y=(int)(rayCenter*bpZoom);
697 ybottomround=-y;
698 } else {
699 ybottomround=-halfDim+1;
700 }
701 toffs=rayCenter-si+co;
702 tPow2=rayCenter*rayCenter*bpZoom*bpZoom;
703 right_top_limit=halfDim-1;
704 for(; y>=ybottomround; y--) {
705 yph=0.5+(float)y;
706 xrightround=(int)sqrt(tPow2 - yph*yph) + 1;
707 if(xrightround>=halfDim) {
708 xrightround=right_top_limit;
709 x=-halfDim;
710 } else {
711 x=-xrightround;
712 }
713 t=toffs-(float)y*sir+((float)(x+1))*cor;
714 imgp=imgOrigin-y*imgDim+x;
715 for(; x<=xrightround; x++, t+=cor) {
716 ti=(int)t; dif=t-(float)ti;
717 *imgp++ += prj[ti] + (prj[ti+1] - prj[ti])*dif;
718 }
719 }
720}

Referenced by fbp(), and imgFBP().

◆ fbp_fft_bidir_complex_radix2()

void fbp_fft_bidir_complex_radix2 ( float * real,
float * imag,
int direct,
int n,
float * sine,
float * cosine )

In-place bidirectional radix-2 discrete FFT of complex data, applying pre-computed sine and cosine tables.

Adapted from subroutine FOUREA listed in Programs for Digital Signal Processing Edited by Digital Signal Processing Committee IEEE Acoustics Speech and Signal Processing Committee Chapter 1 Section 1.1 Page 1.1-4,5

Author
Sakari Alenius
Parameters
realArray of real part complex data
imagArray of imaginary part complex data
directdirection: 0 and 1=forward FFT, -1=inverse FFT
nArray length
sineSine table
cosineCosine table

Definition at line 512 of file fbp.c.

525 {
526 int i, j, m, halfn, mmax, istep, ind;
527 float c, s, treal, timag;
528
529 halfn=n/2;
530
531
532#if(0)
533 for(i=j=1; i<=n; i++) {
534 if(i<j) {
535 treal=real[j-1]; timag=imag[j-1];
536 real[j-1]=real[i-1]; imag[j-1]=imag[i-1];
537 real[i-1]=treal; imag[i-1]=timag;
538 }
539 m=halfn;
540 while(j>m) {j-=m; m++; m/=2;}
541 j+=m;
542 }
543
544 mmax=1;
545 while(mmax<n) {
546 istep=2*mmax;
547 for(m=1; m<=mmax; m++) {
548 ind=(int)((float)(halfn*(m-1))/(float)mmax);
549 c=cosine[ind]; s=(direct>=0)?sine[ind]:-sine[ind];
550 for(i=m; i<=n; i+=istep) {
551 j=i+mmax;
552 treal=real[j-1]*c - imag[j-1]*s;
553 timag=imag[j-1]*c + real[j-1]*s;
554 real[j-1]=real[i-1]-treal;
555 imag[j-1]=imag[i-1]-timag;
556 real[i-1]+=treal;
557 imag[i-1]+=timag;
558 }
559 }
560 mmax=istep;
561 }
562
563#else
564 for(i=j=0; i<n; i++) {
565 if(i<j) {
566 treal=real[j]; timag=imag[j];
567 real[j]=real[i]; imag[j]=imag[i];
568 real[i]=treal; imag[i]=timag;
569 }
570 m=halfn;
571 while(j>=m) {j-=m; m++; m/=2;}
572 j+=m;
573 }
574
575 mmax=1;
576 while(mmax<n) {
577 istep=2*mmax;
578 for(m=0; m<mmax; m++) {
579 ind=(int)((float)(halfn*m)/(float)mmax);
580 c=cosine[ind]; s=(direct>=0)?sine[ind]:-sine[ind];
581 for(i=m; i<n; i+=istep) {
582 j=i+mmax;
583 treal=real[j]*c - imag[j]*s;
584 timag=imag[j]*c + real[j]*s;
585 real[j]=real[i]-treal;
586 imag[j]=imag[i]-timag;
587 real[i]+=treal;
588 imag[i]+=timag;
589 }
590 }
591 mmax=istep;
592 }
593#endif
594}

Referenced by fbp(), and imgFBP().

◆ fbpMakeFilter()

int fbpMakeFilter ( float * sinFFT,
float * cosFFT,
float cutoff,
int filter,
int lenFFT,
int views,
float * fbpFilter,
int verbose )

Compute filter coefficients for FBP reconstruction.

The filter is made by multiplying a ramp by a window function, according to the filter type.

Returns
Returns 0 if ok.
Parameters
sinFFTSine table for FFT
cosFFTCosine table for FFT
cutoffNoise cut-off
filterFilter code: 0=None, 1=Ramp, 2=Butter, 3=Hann, 4=Hamm, 5=Parzen, 6=Shepp.
lenFFTlength of FFT array; must be 256, 512, 1024, 2048, ...
viewsnr of scan views (dim y)
fbpFilterfilter array for FBP (output); allocated outside with size lenFFT+1
verboseVerbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 382 of file fbp.c.

399 {
400 if(verbose>0)
401 printf("fbpMakeFilter(..., %g, %d, %d, %d, ...)\n",
402 cutoff, filter, lenFFT, views);
403 if(sinFFT==NULL || cosFFT==NULL || fbpFilter==NULL) return(1);
404
405 int i, istop;
406 float *ramp, f;
407 float invNr, bp_ifft_scale, rampDC, rampSlope;
408
409 /* Check arguments */
410 if(lenFFT<2 || views<2) return(1);
411
412 /*
413 * Make the ramp to be used in building the reconstruction filter.
414 * Multiplied later by a window to get the actual filter.
415 * Ramp filter:
416 * filt[i] = i/N 0 <= i <= N/2
417 * filt[i] = (N-i)/N N/2+1 <= i <= N-1
418 */
419 float *locData=calloc(lenFFT, sizeof(float));
420 if(locData==NULL) return(2);
421 ramp=locData;
422
423 invNr=1.0/(float)lenFFT;
424 bp_ifft_scale=M_PI/(float)(views*lenFFT);
425 if(verbose>1) printf("invNr=%g bp_ifft_scale=%g\n", invNr, bp_ifft_scale);
426 /* Construct the direct method (original VAX) to calculate the ramp function: */
427 ramp[0]=0.0;
428 for(i=1; i<=lenFFT/2; i++)
429 ramp[i]=ramp[lenFFT-i]=(float)i*invNr;
430 rampDC=ramp[0];
431 rampSlope=ramp[lenFFT/2]-rampDC; // height of ramp
432 if(verbose>1) printf("rampDC=%g rampSlope=%g\n", rampDC, rampSlope);
433 /* FFT filter is made by multiplying a ramp by a window function, */
434 /* according to the filter type */
435 istop=(int)(cutoff*(float)lenFFT+0.5);
436 if(istop>lenFFT/2) istop=lenFFT/2;
437
438 switch(filter) {
439 case FBP_FILTER_RAMP: /* Ramp */
440 fbpFilter[0]=ramp[0]*bp_ifft_scale;
441 for(i=1; i<=istop; i++)
442 fbpFilter[i]=fbpFilter[lenFFT-i]=ramp[i]*bp_ifft_scale;
443 for(; i<=lenFFT/2; i++)
444 fbpFilter[i]=fbpFilter[lenFFT-i]=0.0;
445 break;
446 case FBP_FILTER_HANN: /* Hann */
447 fbpFilter[0]=ramp[0]*bp_ifft_scale;
448 for(i=1; i<=istop; i++) {
449 f=0.5+0.5*cosf(M_PI*invNr*(float)i/cutoff);
450 fbpFilter[i]=fbpFilter[lenFFT-i]=ramp[i]*f*bp_ifft_scale;
451 }
452 for(; i<=lenFFT/2; i++) fbpFilter[i]=fbpFilter[lenFFT-i]=0.0;
453 break;
454 case FBP_FILTER_HAMM: /* Hamm */
455 fbpFilter[0]=ramp[0];
456 for(i=1; i<=istop; i++) {
457 f=0.54+0.46*cosf(M_PI*invNr*(float)i/cutoff);
458 fbpFilter[i]=fbpFilter[lenFFT-i]=ramp[i]*f*bp_ifft_scale;
459 }
460 for(; i<=lenFFT/2; i++) fbpFilter[i]=fbpFilter[lenFFT-i]=0.0;
461 break;
462 case FBP_FILTER_PARZEN: /* Parzen */
463 fbpFilter[0]=ramp[0];
464 for(i=1; i<=istop/2; i++) {
465 f=(float)i/(float)istop; f=1.0-6.0*f*f*(1.0-f);
466 fbpFilter[i]=fbpFilter[lenFFT-i]=ramp[i]*f*bp_ifft_scale;
467 }
468 for(; i<=istop; i++) {
469 f=1.0-(float)i/(float)istop; f=2.0*f*f*f;
470 fbpFilter[i]=fbpFilter[lenFFT-i]=ramp[i]*f*bp_ifft_scale;
471 }
472 for(; i<=lenFFT/2; i++) fbpFilter[i]=fbpFilter[lenFFT-i]=0.0;
473 break;
474 case FBP_FILTER_SHEPP: /* Shepp */
475 fbpFilter[0]=ramp[0];
476 for(i=1; i<=istop; i++) {
477 f=0.5*(float)i/(float)istop; f=sinf(f*M_PI)/(f*M_PI);
478 fbpFilter[i]=fbpFilter[lenFFT-i]=ramp[i]*f*bp_ifft_scale;
479 }
480 for(; i<=lenFFT/2; i++) fbpFilter[i]=fbpFilter[lenFFT-i]=0.0;
481 break;
482 case FBP_FILTER_BUTTER: /* Butter */
483 fprintf(stderr, "Warning: Butter filter not implemented, using None.\n");
484 /* __attribute__((fallthrough)); */
485 /* FALLTHROUGH */
486 case FBP_FILTER_NONE: /* None */
487 for(i=0; i<lenFFT; i++) fbpFilter[i]=bp_ifft_scale;
488 rampSlope=0.0; /* Reset the ramp slope value */
489 break;
490 default:
491 fprintf(stderr, "Error: invalid filter requested.\n");
492 free(locData);
493 return(1);
494 } /* switch */
495
496 free(locData);
497 return(0);
498}

Referenced by fbp(), and imgFBP().

◆ imgFBP()

int imgFBP ( IMG * scn,
IMG * img,
int imgDim,
float zoom,
int filter,
float cutoff,
float shiftX,
float shiftY,
float rotation,
int verbose )

Filtered back-projection (FBP) reconstruction using data in IMG struct.

See also
fbp, imgReprojection
Returns
Returns 0 if ok.
Parameters
scnSinogram (input) data.
imgImage (output) data; allocated here.
imgDimImage dimension (size, usually 128 or 256); must be an even number.
zoomZoom factor (for example 2.45)
filterFilter code: 0=None, 1=Ramp, 2=Butter, 3=Hann, 4=Hamm, 5=Parzen, 6=Shepp.
cutoffNoise cut-off (for example 0.3).
shiftXPossible shifting in x dimension (mm).
shiftYPossible shifting in y dimension (mm).
rotationPossible image rotation (in degrees).
verboseVerbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 108 of file fbp.c.

130 {
131 if(verbose>0)
132 printf("imgFBP(scn, img, %d, %g, %d, %g, %g, %g, %g)\n",
133 imgDim, zoom, filter, cutoff, shiftX, shiftY, rotation);
134
135 int i, j, ret, plane, frame, fullBP;
136 int lenFFT, halfDim;
137 float *sinB, *sinBrot, *sinFFT, *cosFFT;
138 float *oReal, *oImag, *real, *imag, *scnData, *imgData, *imgOrigin;
139 float *fbpFilter, bpZoom, bpZoomInv, *fptr, *gptr, *row1, *row2, f;
140 float offsX, offsY, pixSize; /* note that pixSize is in cm */
141
142
143 /* Check the arguments */
144 if(scn->status!=IMG_STATUS_OCCUPIED) return(1);
145 if(imgDim<2 || imgDim>4096 || imgDim%2) return(1);
146 if(zoom<0.01 || zoom>1000.) return(1);
147 if(filter<0 || filter>6) return(1);
148 if(cutoff<=0.0 || cutoff>=1.0) return(1);
149 if(scn->dimx<=1 || scn->dimx>16384) return(1);
150
151 /*
152 * Allocate output image
153 */
154 if(verbose>1) {printf("allocating memory for the image\n"); fflush(stdout);}
155 imgEmpty(img);
156 ret=imgAllocate(img, scn->dimz, imgDim, imgDim, scn->dimt);
157 if(ret) return(3);
158
159 /* Set image "header" information */
160 if(verbose>1) {printf("setting image header\n"); fflush(stdout);}
161 img->type=IMG_TYPE_IMAGE;
162 img->unit=CUNIT_COUNTS; /* This will convert to ECAT counts/sec */
163 img->scanStart=scn->scanStart;
164 img->axialFOV=scn->axialFOV;
165 img->transaxialFOV=scn->transaxialFOV;
166 img->sizez=scn->sizez;
167 strcpy(img->studyNr, scn->studyNr);
168 img->sampleDistance=scn->sampleDistance;
169 if(scn->sampleDistance<=0.0) {
170 if(scn->dimx==281) {
171 img->sizez=4.25;
172 scn->sampleDistance=1.95730;
173 scn->axialFOV=150.; scn->transaxialFOV=550.;
174 } else {
175 img->sizez=6.75;
176 scn->sampleDistance=3.12932;
177 scn->axialFOV=108.; scn->transaxialFOV=600.829;
178 }
179 }
180 pixSize=scn->sampleDistance*(float)scn->dimx/((float)imgDim*zoom);
181 img->sizex=img->sizey=pixSize;
182 for(plane=0; plane<scn->dimz; plane++)
183 img->planeNumber[plane]=scn->planeNumber[plane];
184 strcpy(img->radiopharmaceutical, scn->radiopharmaceutical);
185 img->isotopeHalflife=scn->isotopeHalflife;
186 img->decayCorrection=IMG_DC_NONCORRECTED;
187 for(frame=0; frame<scn->dimt; frame++) {
188 img->start[frame]=scn->start[frame]; img->end[frame]=scn->end[frame];
189 img->mid[frame]=0.5*(img->start[frame]+img->end[frame]);
190 }
191 img->isWeight=0;
192
193 /*
194 * Preparations for reconstruction
195 */
196 /* Allocate memory for scan and image data arrays */
197 if(verbose>1) {printf("allocating memory for matrices\n"); fflush(stdout);}
198 i=scn->dimx*scn->dimy + img->dimx*img->dimy;
199 float *matData=calloc(i, sizeof(float));
200 if(matData==NULL) {
201 imgEmpty(img);
202 return(4);
203 }
204 scnData=matData; imgData=matData+(scn->dimx*scn->dimy);
205
206 /* Set FFT array length */
207 lenFFT=(2*scn->dimx)-1;
208 if(lenFFT<256) lenFFT=256; else if(lenFFT<512) lenFFT=512;
209 else if(lenFFT<1024) lenFFT=1024; else if(lenFFT<2048) lenFFT=2048;
210 else if(lenFFT<4096) lenFFT=4096; else if(lenFFT<8192) lenFFT=8192;
211 else if(lenFFT<16384) lenFFT=16384; else return(1);
212 if(verbose>2)
213 printf("rays=%d views=%d lenFFT=%d\n", scn->dimx, scn->dimy, lenFFT);
214
215 /* Allocate memory for pre-computed tables */
216 i=3*scn->dimy/2 + 3*scn->dimy/2 + 5*lenFFT/4;
217 if(verbose>1) printf("allocating memory for pre-computed tables\n");
218 float *tabData=calloc(i, sizeof(float));
219 if(tabData==NULL) {
220 imgEmpty(img); free(matData);
221 return(4);
222 }
223 sinB=tabData; sinBrot=tabData+3*scn->dimy/2;
224 sinFFT=sinBrot+3*scn->dimy/2;
225
226 /* Pre-compute the sine tables for back-projection (note the rotation!) */
227 if(verbose>1) printf("computing sine tables for back-projection\n");
228 recSinTables(scn->dimy, sinB, sinBrot, rotation);
229
230 /* Pre-compute the sine and cosine tables for FFT */
231 if(verbose>1) printf("computing sine and cosine tables for FFT\n");
232 {
233 double df=0.0, dg;
234 dg=(double)M_PI*2.0/(double)lenFFT;
235 for(i=0; i<5*lenFFT/4; i++, df+=dg) sinFFT[i]=(float)sin(df);
236 cosFFT=sinFFT+lenFFT/4;
237 }
238
239 /* Allocate memory for reconstruction filter coefficients */
240 if(verbose>1) printf("allocating memory for reconstruction filters\n");
241 float *filData=calloc((5*lenFFT), sizeof(float));
242 if(filData==NULL) {
243 imgEmpty(img); free(matData); free(tabData);
244 return(4);
245 }
246 fbpFilter=filData; oReal=filData+lenFFT; oImag=oReal+2*lenFFT;
247 real=oReal+lenFFT/2; imag=oImag+lenFFT/2;
248
249 /* Compute reconstruction filter coefficients */
250 ret=fbpMakeFilter(sinFFT, cosFFT, cutoff, filter, lenFFT, scn->dimy,
251 fbpFilter, verbose-1);
252 if(ret) {
253 imgEmpty(img); free(matData); free(tabData); free(filData);
254 return(5);
255 }
256 if(verbose>3) {
257 printf("\nfbpFilter:\n");
258 for(i=0; i<lenFFT; i++) {
259 printf("%g ", fbpFilter[i]);
260 if((i+1)%7==0 || i==lenFFT-1) printf("\n");}
261 printf("\n");
262 }
263
264 /* Set the backprojection zoom and inverse (globals) */
265 bpZoom=zoom*(float)imgDim/(float)scn->dimx;
266 bpZoomInv=1.0/bpZoom;
267 if(verbose>2) printf("bpZoom=%g bpZoomInv=%g\n", bpZoom, bpZoomInv);
268
269 /* Check if image corners need to be processed */
270 if(zoom>M_SQRT2) fullBP=1; /* profile overlaps the image; 1.41421356 */
271 else fullBP=0; /* do not process image corners */
272 if(verbose>1) printf("fullBP := %d\n", fullBP);
273
274 /* Initialize variables used by back-projection */
275 if(verbose>1) printf("initialize variables for back-projection\n");
276 halfDim=imgDim/2;
277 offsX=shiftX/pixSize; offsY=shiftY/pixSize;
278 for(i=0; i<3*(scn->dimy)/2; i++) {
279 sinB[i]*=bpZoomInv;
280 sinBrot[i]*=bpZoomInv;
281 }
282 imgOrigin=imgData+imgDim*(halfDim-1)+halfDim;
283 if(verbose>2) {
284 printf("halfDim=%d offsX=%g offsY=%g\n", halfDim, offsX, offsY);
285 }
286
287 /*
288 * Reconstruct one matrix at a time
289 */
290 if(verbose>1) printf("reconstruct one matrix at a time...\n");
291 for(plane=0; plane<scn->dimz; plane++)
292 for(frame=0; frame<scn->dimt; frame++) {
293
294 if(verbose>3) {
295 printf("reconstructing plane %d frame %d\n",
296 scn->planeNumber[plane], frame+1);
297 fflush(stdout);
298 }
299
300 /* Copy scan data into the array */
301 for(i=0, fptr=scnData; i<scn->dimy; i++)
302 for(j=0; j<scn->dimx; j++)
303 *fptr++=scn->m[plane][i][j][frame];
304 if(verbose>8)
305 printf("scn->m[%d][158][116][%d]=%e\n",
306 plane, frame, scn->m[plane][64][56][frame]);
307
308 /* Initiate image buffer */
309 for(i=0, fptr=imgData; i<imgDim*imgDim; i++) *fptr++=0.0;
310
311 /* Process each sinogram row (projection), two at the same time */
312 for(int view=0; view<scn->dimy; view+=2) {
313
314 /* Fill FFT buffers with two rows */
315 for(j=0; j<2*lenFFT; j++) oReal[j]=0.0; /* zero padding */
316 for(j=0; j<2*lenFFT; j++) oImag[j]=0.0; /* zero padding */
317 row1=scnData+view*scn->dimx;
318 row2=row1+scn->dimx;
319 fptr=real; gptr=imag;
320 for(j=0; j<scn->dimx; j++) {
321 *fptr++ = *row1++;
322 *gptr++ = *row2++;
323 }
324
325 /* Forward FFT */
326 fptr=real; gptr=imag;
327 fbp_fft_bidir_complex_radix2(fptr, gptr, 1, lenFFT, sinFFT, cosFFT);
328
329 /* Due to symmetry properties, filtering the combined data together */
330 for(j=0, fptr=fbpFilter; j<lenFFT; j++) {
331 real[j] *= *fptr;
332 imag[j] *= *fptr++;
333 }
334
335 /* Inverse FFT */
336 fptr=real; gptr=imag;
337 fbp_fft_bidir_complex_radix2(fptr, gptr, -1, lenFFT, sinFFT, cosFFT);
338 fptr=real; gptr=imag;
339 if(fullBP) {
340 fbp_back_proj(fptr, imgData, imgDim, view, scn->dimy, scn->dimx,
341 offsX, offsY, sinB, sinBrot);
342 fbp_back_proj(gptr, imgData, imgDim, view+1, scn->dimy, scn->dimx,
343 offsX, offsY, sinB, sinBrot);
344 } else {
345 fbp_back_proj_round(fptr, imgOrigin, imgDim, view,
346 scn->dimy, scn->dimx, offsX, offsY, bpZoom, sinB, sinBrot);
347 fbp_back_proj_round(gptr, imgOrigin, imgDim, view+1,
348 scn->dimy, scn->dimx, offsX, offsY, bpZoom, sinB, sinBrot);
349 }
350
351 } /* next projection (view) */
352
353 /* Copy the image array to matrix data */
354 /* At the same time, correct for the frame length */
355 f=img->end[frame]-img->start[frame];
356 if(f<1.0) f=1.0;
357 f=1.0/f;
358 for(i=0, fptr=imgData; i<img->dimy; i++)
359 for(j=0; j<img->dimx; j++)
360 img->m[plane][i][j][frame] = (*fptr++)*f;
361 if(verbose>8)
362 printf("img->m[%d][64][56][%d]=%e\n",
363 plane, frame, img->m[plane][64][56][frame]);
364
365 } /* next matrix */
366
367 /* Free the memory */
368 free(matData); free(tabData); free(filData);
369
370 if(verbose>1) printf("imgFBP() done.\n");
371 return(0);
372}
fbp_back_proj
void fbp_back_proj(float *prj, float *imgCorner, int imgDim, int view, int views, int rays, float offsX, float offsY, float *sinB, float *sinBrot)
Definition fbp.c:599
imgAllocate
int imgAllocate(IMG *image, int planes, int rows, int columns, int frames)
Definition img.c:194
imgEmpty
void imgEmpty(IMG *image)
Definition img.c:121
IMG_STATUS_OCCUPIED
#define IMG_STATUS_OCCUPIED
Definition libtpcimgio.h:1479
IMG_DC_NONCORRECTED
#define IMG_DC_NONCORRECTED
Definition libtpcimgio.h:1533
IMG_TYPE_IMAGE
#define IMG_TYPE_IMAGE
Definition libtpcimgio.h:1493
IMG::sizex
float sizex
Definition libtpcimgio.h:1714
IMG::dimx
unsigned short int dimx
Definition libtpcimgio.h:1780
IMG::type
char type
Definition libtpcimgio.h:1704
IMG::sampleDistance
float sampleDistance
Definition libtpcimgio.h:1712
IMG::m
float **** m
Definition libtpcimgio.h:1799
IMG::decayCorrection
char decayCorrection
Definition libtpcimgio.h:1688
IMG::transaxialFOV
float transaxialFOV
Definition libtpcimgio.h:1710
IMG::unit
char unit
Definition libtpcimgio.h:1673
IMG::status
char status
Definition libtpcimgio.h:1665
IMG::scanStart
time_t scanStart
Definition libtpcimgio.h:1692
IMG::dimt
unsigned short int dimt
Definition libtpcimgio.h:1778
IMG::planeNumber
int * planeNumber
Definition libtpcimgio.h:1809
IMG::sizey
float sizey
Definition libtpcimgio.h:1716
IMG::start
float * start
Definition libtpcimgio.h:1815
IMG::dimz
unsigned short int dimz
Definition libtpcimgio.h:1784
IMG::dimy
unsigned short int dimy
Definition libtpcimgio.h:1782
IMG::end
float * end
Definition libtpcimgio.h:1817
IMG::radiopharmaceutical
char radiopharmaceutical[32]
Definition libtpcimgio.h:1683
IMG::isotopeHalflife
float isotopeHalflife
Definition libtpcimgio.h:1685
IMG::studyNr
char studyNr[MAX_STUDYNR_LEN+1]
Definition libtpcimgio.h:1677
IMG::axialFOV
float axialFOV
Definition libtpcimgio.h:1708
IMG::mid
float * mid
Definition libtpcimgio.h:1819
IMG::isWeight
char isWeight
Definition libtpcimgio.h:1825
IMG::sizez
float sizez
Definition libtpcimgio.h:1718
Generated on for TPCCLIB by doxygen 1.15.0