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_bpix.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 subroutine to perform a task similar to bwray() but pixels instead of the rays
23 void Blob_class::bpix(INTEGER nv, INTEGER nh, INTEGER *list, REAL* weight,
24 INTEGER *numb, REAL *snorm)
29 INTEGER mf, hf, ml, hl;
37 REAL dinv2 = (REAL) 2.0 / Blob.delta;
38 REAL dinv2sq3 = dinv2 * Consts.isqrt3;
40 REAL dist, distx, disty; // dist is the distance betw a blob center and a pixel center. (distx, disty) is the vector distance
41 REAL dispx, dispy; // (dispx, dispy) is the displacement from the center of a blob to the center of the nearest pixel.
43 multp = (INTEGER) floor(Blob.support / Blob.delta + 0.5) + 1;
45 H2 = (INTEGER) (Blob.H / 2);
46 M2 = (INTEGER) (Blob.M / 2);
47 nelem2 = (INTEGER) (GeoPar.nelem / 2);
49 REAL* footpx = new REAL[4 * multp + 1];
50 REAL* footpy = new REAL[4 * multp + 1];
52 for (m0 = -2 * multp; m0 <= 2 * multp; m0++)
54 footpx[m0 + 2 * multp] = (REAL) m0 * Blob.di2;
57 for (h0 = -2 * multp; h0 <= 2 * multp; h0++)
59 footpy[h0 + 2 * multp] = (REAL) h0 * Blob.di2sq3;
62 // coord of the center of the pixel (nh, nv)
63 xc = GeoPar.pixsiz * (nh - nelem2);
64 yc = GeoPar.pixsiz * (nelem2 - nv);
66 mx = (INTEGER) floor(dinv2 * xc + 0.5); // coord of the center of
68 hy = (INTEGER) floor(dinv2sq3 * yc + 0.5); // the nearest blob
73 // displacement of the pixel footprint with resp to the nearest blob center
74 dispx = di2 * mx - xc;
75 dispy = di2sq3 * hy - yc;
77 mf = MAX0(lmf, mc - 2 * multp);
78 ml = MIN0(lml, mc + 2 * multp);
80 hf = MAX0(lhf, hc - 2 * multp);
81 hl = MIN0(lhl, hc + 2 * multp);
86 if (ml < lmf || mf > lml || hl < lhf || hf > lhl)
88 for (h = hf; h <= hl; h++)
91 disty = footpy[h0 + 2 * multp] + dispy; //the y in the x-y coord and the h in the m-h coord are in the same direction
92 if (fabs(disty) > Blob.support)
95 for (m = mf; m <= ml; m++)
98 distx = footpx[m0 + 2 * multp] + dispx;
99 if (fabs(distx) > Blob.support)
102 if ((m + h) % 2 == 0)
106 ind = h1 * Blob.M + m1;
108 dist = sqrt(distx * distx + disty * disty);
109 if (dist < Blob.support)
111 n = (INTEGER) round(dist * sampiblob); // changed "rint" to "round". Lajos, Jan 25, 2005
113 weight[(*numb)] = Blob.vtable[n];
115 (*snorm) += (Blob.vtable[n] * Blob.vtable[n]);
122 delete[] footpx; // bug 92 - Lajos - 03/02/2005
123 delete[] footpy; // bug 92 - Lajos - 03/02/2005