2 ***********************************************************
3 $SNARK_Header: S N A R K 1 4 - A PICTURE RECONSTRUCTION PROGRAM $
4 $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/lino.cpp $
5 $LastChangedRevision: 80 $
6 $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
8 ***********************************************************
10 ******************************************************************
12 - LINOGRAM RECONSTRUCTION ALGORITHM -
14 THIS PROGRAM IMPLEMENTS AN IMAGE RECONSTRUCTION TECHNIQUE
15 DERIVED BY PAUL EDHOLM AND GABOR T. HERMAN. DATA IS COLLECTED
16 IN A SPECIAL GEOMETRY,AND THE IMAGE IS THEN RECONSTRUCTED
17 (ESSENTIALLY) BY USING 3 ONE-DIMENSIONAL FOURIER TRANSFORMS.
19 VERSION 1.0 : LYLE BERKOWITZ AUGUST, 1986
20 1.1 : DAVID A. ROBERTS MAY, 1987
21 : RESTRUCTURED PROGRAM
22 ADDED GRIDDING ALGORITHM
23 1.2 : DAVID A. ROBERTS JUNE, 1987
24 : USE OF RTFORT TO IMPROVE SPEED
25 1.3 : DAVID A. ROBERTS JULY, 1987
26 : USE OF CHIRP-Z ALGORITHM
28 1.4 : D. ODHNER JAN, 1989
29 ******************************************************************
48 BOOLEAN lino_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
51 static const REAL padwgt = 3.0;
55 INTEGER i, j, k, r1, ii, iii, iflag, i3, kk;
56 INTEGER b, word, filter, n3[3], np, np1, np2;
59 INTEGER glnum, gl1num, imgnum;
60 REAL w, r, theta, sinth, costh, dold, dnew, dist, rinc;
61 REAL sum, diff, c, cutoff, factr1, factr2, factr3;
70 // modified calls to getwrd() to use the 4-parameter version. Lajos, Dec 13, 2004
71 static const INTEGER fname[3] =
72 { CHAR2INT('b', 'a', 'n', 'd'), CHAR2INT('s', 'i', 'n', 'c'), CHAR2INT('c',
77 //commented out since it is never used. Lajos, Dec 16, 2004
78 //static const INTEGER pad = CHAR2INT('p','a','d',' ');
82 //-------------------------------------------------------------------
83 //... START MAIN PROGRAM
84 //------------------------------------------------------------------
90 // READ IN FILTER NAME FROM USER INPUT FILE
92 word = InFile.getwrd(TRUE, &eol, fname, 3);
93 for (filter = 0; filter < 3; filter++)
95 if (word == fname[filter])
99 fprintf(output, "\n **** %4s is not a valid filter name", int2str(word));
100 fprintf(output, "\n **** LINOGRAM execution aborted\n");
103 L0003: cutoff = InFile.getnum(FALSE, &eol);
107 fprintf(output, "\n **** not enough parameters");
108 fprintf(output, "\n **** LINOGRAM execution aborted\n");
113 cutoff = (REAL) 0.0001;
115 fprintf(output, "\n linogram reconstruction algorithm");
116 fprintf(output, "\n using filter %4s", int2str(word));
117 fprintf(output, "\n cutoff = %9.4f", cutoff);
119 // INITIALIZE LOOP INVARIANTS
121 nnew = (GeoPar.nelem - 1) / 2;
125 // Decide whether to pad
128 while ((1 << k) < (GeoPar.nelem + 1))
133 if (((1 << k) + padwgt * 2 * k)
134 <= (GeoPar.nelem + 1 + padwgt * sumfac(GeoPar.nelem + 1)))
138 while (n3[1] < (GeoPar.nelem + 1))
146 n3[1] = GeoPar.nelem + 1;
149 dold = (REAL) 2.0 / ((REAL) (GeoPar.usrays));
150 dist = 1 / (2 * n3[1] * GeoPar.pixsiz);
151 gl1num = 2 * n3[1] + 2;
152 glnum = gl1num * GeoPar.usrays;
153 imgnum = GeoPar.nelem * gl1num;
155 // G IS THE BASE OF PROJECTION DATA
157 g = new REAL[GeoPar.usrays + 2];
159 //******************************************************C
160 //... BEGIN THE RECONSTRUCTION
161 //******************************************************C
163 // FIRST LOOP ON SETS OF PROJECTION ANGLES
164 // IFLAG = 0 FOR THE FIRST SET OF ANGLES (T)
165 // IFLAG = 1 FOR THE SECOND SET OF ANGLES (C)
167 for (iflag = 0; iflag <= 1; iflag++)
169 // GET THE FIRST AND LAST PROJECTION NUMBERS
179 np2 = 2 * GeoPar.usrays;
182 img = new REAL[imgnum];
183 for (ii = 0; ii < imgnum; ii++)
188 gl = new REAL[glnum];
189 for (ii = 0; ii < glnum; ii++)
194 gl1 = new REAL[gl1num];
195 for (ii = 0; ii < gl1num; ii++)
201 for (np = np1; np < np2; np++)
203 for (ii = 0; ii < GeoPar.usrays + 2; ii++)
208 ProjFile.ReadProj(np, g, GeoPar.usrays);
210 Anglst.getang(np, &theta, &sinth, &costh);
212 //-------------------------------------------------------------------C
213 //... 1) READ IN THE PROJECTION DATA BACKWARDS
214 //... 2) MAKE SURE GL1 ARRAY IS IN PROPER FORMAT BEFORE FFT
215 //... (I.E. WITH THE ZERO VALUE AT THE BEGINNING)
216 //-------------------------------------------------------------------C
219 for (i = GeoPar.nelem; i >= 0; i--)
225 for (i = GeoPar.usrays + 1; i <= 2 * n3[1]; i++)
231 for (i = GeoPar.usrays; i >= GeoPar.nelem + 2; i--)
237 //------------------------------------------------------------C
238 //... TAKE FIRST 1-D FOURIER TRANSFORM; USRAY-DIRECTION
239 //... ** N3(2) IS ASSIGNED HALF THE NUMBER OF ELEMENTS **
240 //------------------------------------------------------------C
244 //------------------------------------------------------------C
245 //... MULTIPLY BY FILTER AND WEIGHT (AND SAMPLE SPACING)
246 //------------------------------------------------------------C
248 w = (REAL) fabs(sinth);
250 w = (REAL) fabs(costh);
251 c = cutoff / (2 * GeoPar.pinc * w);
255 for (r1 = 0; r1 <= n3[1]; r1++)
257 if (r <= Consts.zero)
264 factr1 = GeoPar.pinc * w * w;
266 factr2 = qfilt(r, c, 1) / r;
269 factr2 = qfilt(r, c, 1) / r;
270 factr2 = factr2 * factr2;
272 factr3 = qfilt(r, c, filter);
273 gl1[j] = factr1 * factr2 * factr3 * gl1[j];
274 gl1[j + 1] = factr1 * factr2 * factr3 * gl1[j + 1];
277 //-------------------------------------------------------------------C
278 //... STORE DATA IN NEW ARRAY (2-DIMENSIONAL): [NELEM+2,4NELEM+2]
279 //-------------------------------------------------------------------C
281 gl[(r1 + 1) * 2 * GeoPar.usrays - b] = gl1[j];
282 gl[(r1 + 1) * 2 * GeoPar.usrays - b + 1] = gl1[j + 1];
289 //************************************C
290 // END OF LOOP ON PROJECTIONS
291 //************************************C
293 //------------------------------------------------------------------
294 //... PERFORM GRIDDING (INCORPORATES A 1-D FFT IN ANGULAR VARIABLE)
295 //... STORE OUTPUT OF GRIDDING IN TWO DIMENSIONAL ARRAY:
296 //... [NELEM+1,2*NELEM]
297 //-------------------------------------------------------------------
299 otemp = new REAL[2 * GeoPar.nelem];
302 for (n = 0; n < n3[1] + 1; n++)
307 // SPECIAL CASE WHEN n = 1
308 // DUE TO NATURE OF FILTERING, DC VALUE WILL ALWAYS BE ZERO.
310 for (kk = 0; kk < 2 * GeoPar.nelem; kk++)
318 // **** CHANGED DNEW FOR RTFORT! ****
319 // (DUE TO PADDING WITH AN EXTRA ZERO)
321 dnew = ((REAL) (n)) / (2 * n3[1]);
323 // **** USE CHIRP-Z ****
325 temp = &(gl[(n) * 2 * GeoPar.usrays]);
326 czt(temp, otemp, dold, dnew, GeoPar.nelem, nnew);
329 //-------------------------------------------------------------------
330 //... WRITE 'OTEMP' INTO THE 2D ARRAY 'IMG' ON A COLUMN-BY-COLUMN BASIS
331 //-------------------------------------------------------------------
334 for (k = 0; k < GeoPar.nelem; k++)
336 img[k * gl1num + iii] = otemp[ii];
337 img[k * gl1num + iii + 1] = otemp[ii + 1];
342 //-----------------------------------------------------------------*
343 //... PERFORM AN INVERSE FFT ALONG ROWS IN THE ARRAY 'IMG'.
344 //... THE FOURIER SAMPLES ARE ARRANGED IN EACH ROW STARTING FROM
346 //-----------------------------------------------------------------*
348 // LOOP ON ROWS OF 'IMG' (SHOULD BE 'NELEM' ROWS)
350 for (n = 0; n < GeoPar.nelem; n++)
352 rtfort(&(img[n * gl1num]), n3, -1);
354 //----------------------------------------------C
355 // MULTIPLY BY SAMPLE SPACING AFTER FFT:
356 //----------------------------------------------C
358 for (k = 0; k < 2 * n3[1]; k++)
360 img[(n) * gl1num + k] /= (GeoPar.pixsiz);
364 //-----------------------------------------------------------------*
365 //... WE NOW HAVE THE VALUES ALONG ROWS/COLUMNS IN THE IMAGE
366 //... IN THE ROWS OF 'IMG'. ALL WE HAVE TO DO NOW IS OUTPUT
367 //... THESE VALUES IN THE PROPER ORDER INTO THE ARRAY 'RECON'.
368 //-----------------------------------------------------------------*
370 // LOOP ON ROWS OF 'IMG' (SHOULD BE 'NELEM' ROWS)
373 for (n = GeoPar.nelem - 1; n >= 0; n--)
376 // LOOP ON LAST 'N' COLUMNS IN EACH ROW [N = (NELEM-1)/2]
379 for (k = 2 * n3[1] - nnew; k < 2 * n3[1]; k++)
383 recon[(n) * GeoPar.nelem + ii] += img[(iii) * gl1num + k];
387 recon[(GeoPar.nelem - ii - 1) * GeoPar.nelem + n] +=
388 img[(iii) * gl1num + k];
393 // LOOP ON FIRST 'N+1' COLUMNS IN EACH ROW
394 // (USING THE SAME INDEX, 'ii', TO INDEX RECON ARRAY)
396 for (k = 0; k < nnew + 1; k++)
400 recon[(n) * GeoPar.nelem + ii] += img[(iii) * gl1num + k];
404 recon[(GeoPar.nelem - ii - 1) * GeoPar.nelem + n] +=
405 img[(iii) * gl1num + k];
418 //*************************************************C
419 // END OF LOOP ON SETS OF PROJECTION ANGLES
420 //*************************************************C
425 for (n = 0; n < GeoPar.area; n++)
429 sum /= (REAL) (GeoPar.nelem * GeoPar.nelem);
431 diff = (REAL) fabs(GeoPar.aveden - sum);
432 for (n = 0; n < GeoPar.area; n++)