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/blob_footp.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 subroutine to perform the forward (1) or backward (-1) "projection" using blob footprints
23 void Blob_class::footp(REAL *blob_array, REAL *pict, INTEGER option)
25 INTEGER mult; //max number of pixels whose center is covered by a blob
26 INTEGER m, h, m1, h1; // (m,h) coord of the center of a blob in the underlying sq grid. m1 * Blob.H + h1 index in the blob array.
27 INTEGER k, j; //(j,k) coord of the center of a pixel
29 INTEGER jc, kc; // index of the nearest pixel
30 INTEGER jf, kf, jl, kl; //upper and lower bound for j and k given a blob
36 REAL ipixsiz = (REAL) 1.0 / GeoPar.pixsiz;
37 REAL dist, distx, disty; // dist is the distance betw a blob center and a pixel center. (distx, disty) is the vector distance
38 REAL dispx, dispy; // (dispx, dispy) is the displacement from the center of a blob to the center of the nearest pixel.
40 mult = (INTEGER) floor(Blob.support / GeoPar.pixsiz + 0.5) + 1;
42 nelem2 = (INTEGER) (GeoPar.nelem / 2);
44 REAL* footx = new REAL[2 * mult + 1];
45 REAL* footy = new REAL[2 * mult + 1];
47 for (k0 = -mult; k0 <= mult; k0++)
49 footy[k0 + mult] = (REAL) k0 * GeoPar.pixsiz;
52 for (j0 = -mult; j0 <= mult; j0++)
54 footx[j0 + mult] = (REAL) j0 * GeoPar.pixsiz;
57 // initialization for forward operation
60 for (k = 0; k < GeoPar.nelem; k++)
62 for (j = 0; j < GeoPar.nelem; j++)
64 pict[k * GeoPar.nelem + j] = 0.0;
69 //initialization for backward operation
70 else if (option == -1)
72 for (j = 0; j <= Blob.area; j++)
81 // forward or backward operation
83 for (h = lhf; h <= lhl; h++)
85 for (m = lmf; m <= lml; m++)
86 { // m and h must have the same parity!!
90 xc = di2 * m; // coord of the center of the blob (m,h)
93 jx = (INTEGER) floor(ipixsiz * xc + 0.5); // coord of the center of
95 ky = (INTEGER) floor(auy + 0.5); // the nearest pixel
97 jc = jx + nelem2; // bring to the right address
98 kc = (INTEGER) floor(-auy + 0.5) + nelem2;
100 // displacement of the footprint with resp to the nearest pixel center
101 dispx = jx * GeoPar.pixsiz - xc;
103 dispy = ky * GeoPar.pixsiz - yc;
105 jf = MAX0(ljf, jc - mult);
106 jl = MIN0(ljl, jc + mult);
108 kf = MAX0(lkf, kc - mult);
109 kl = MIN0(lkl, kc + mult);
111 if (jl < ljf || jf > ljl || kl < lkf || kf > lkl)
117 ind = h1 * Blob.M + m1;
119 for (k = kf; k <= kl; k++)
122 disty = footy[k0 + mult] - dispy; //the y in the x-y coord and the k in the j-k coord are in the opposite directions
123 if (fabs(disty) > Blob.support)
126 for (j = jf; j <= jl; j++)
129 distx = footx[j0 + mult] + dispx;
130 if (fabs(distx) > Blob.support)
133 dist = sqrt(distx * distx + disty * disty);
134 if (dist < Blob.support)
136 n = (INTEGER) round(dist * sampiblob); // changed "rint" to "round". Lajos, Jan 25, 2005
140 pict[k * GeoPar.nelem + j] += (blob_array[ind]
145 else if (option == -1)
147 blob_array[ind] += (pict[k * GeoPar.nelem + j]
159 delete[] footx; // bug 92 - Lajos - 03/02/2005
160 delete[] footy; // bug 92 - Lajos - 03/02/2005