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/emap_CInitMAP.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 ===========================================================================
11 initialize working space for the MAP algorithm
12 ===========================================================================
26 REAL getSjj(INTEGER j, INTEGER nelem);
28 void emap_class::InitMAP(REAL* recon)
32 //------ local variables
33 INTEGER i, first, np, nr, idelta, numb, list[4096], usnr;
34 REAL snorm, weight[4096];
47 idelta = (INTEGER) (MAPEM.delta);
49 //------ set image to average density and initialize the work area for
50 // the total length of all projections through each pixel
53 for (i = 0; i < area; i++)
57 MAPEM.Sjjinv[i] = 1 / getSjj(i, GeoPar.nelem); //blob case will not use Sjj Hstau 6/25/2003
61 { //not needed in the blob case Hstau 6/3/2003
62 for (i = 0; i < area; i++)
69 // compute the total length of all projections through each pixel
71 for (np = 0; np < GeoPar.prjnum; np++)
73 for (first = 0; first < idelta; first++)
75 for (nr = first; nr < GeoPar.nrays; nr += idelta)
79 wray(np, nr, list, weight, &numb, &snorm);
83 Blob.bwray(np, nr, list, weight, &numb, &snorm);
86 for (i = 0; i < numb; i++)
89 if (MAPEM.gamma1 > 0.0)
90 { //when blobs, gamma1 = 0.0
91 MAPEM.Wjj[list[i]] += weight[i] / MAPEM.gamma1
92 * MAPEM.Sjjinv[list[i]];
96 MAPEM.Wjj[list[i]] += weight[i];
103 //------ get projection data
106 for (np = 0; np < GeoPar.prjnum; np++)
108 for (usnr = 0; usnr < GeoPar.usrays; usnr++)
111 //nr = usnr + GeoPar.fusray - 1;
112 nr = usnr + GeoPar.fusray; // mk
114 // T.K. Narayan's modification - 24.Apr.1996
115 // setting negative projection data to 0.0
117 if (Anglst.prdta(np, nr) >= 0.0)
119 //MAPEM.prj[(np-1)*GeoPar.usrays+usnr] = prdta(np,nr);
120 MAPEM.prj[(np) * GeoPar.usrays + usnr] = Anglst.prdta(np, nr);
124 //MAPEM.prj[(np-1)*GeoPar.usrays+usnr] = 0.0;
125 MAPEM.prj[(np) * GeoPar.usrays + usnr] = 0.0;
130 "\n *** WARNING - There exist negative values in the projection data ***"); //bug 127, hstau 7/8/2005
132 "\n *** Negative values set to 0.0 ***");
141 REAL getSjj(INTEGER j, INTEGER nelem)
154 } //add exception when nelem=3,by wei, 1/2005
156 jx = j - ((j - 1) / nelem) * nelem;
157 jy = (j - 1) / nelem + 1;
158 if ((jx > 2) && (jx < nelem - 1))
160 if ((jy > 2) && (jy < nelem - 1))
166 if ((jy == 2) || (jy == nelem - 1))
178 if ((jx == 2) || (jx == nelem - 1))
180 if ((jy > 2) && (jy < nelem - 1))
186 if ((jy == 2) || (jy == nelem - 1))
198 if ((jy > 2) && (jy < nelem - 1))
204 if ((jy == 2) || (jy == nelem - 1))