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_CMAP.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 ======================================================================
11 one iteration of the maximum a-posteriori reconstruction algorithm
12 for emission tomography
13 ======================================================================
26 BOOLEAN emap_class::MAP(REAL* recon)
28 //------ local variables
29 INTEGER i, j, first, np, nr, idelta, usnr;
45 idelta = (INTEGER) (MAPEM.delta);
47 //====== step 1: compute q for iteration k+1:
48 for (j = 0; j < area; j++)
53 for (np = 0; np < GeoPar.prjnum; np++)
55 for (first = 0; first < idelta; first++)
57 for (usnr = first; usnr < GeoPar.usrays; usnr += idelta)
59 nr = usnr + GeoPar.fusray;
60 yi = MAPEM.prj[(np) * GeoPar.usrays + usnr];
63 if (findqpart(np, nr, yi, recon, MAPEM.qj))
64 return TRUE; //if findqpart return false, EMAP is not applicable, bywei 1/2005
70 // MAPEM.idy = MAPEM.idx = 0, see comments in emap.cpp
73 fprintf(output,"\n #IDY=%d #IDX=%d\n", MAPEM.idy, MAPEM.idx);fflush(output);
78 for (i = MAPEM.idy * GeoPar.nelem;
79 i < GeoPar.nelem * (GeoPar.nelem - 1) + 1;
80 i += GeoPar.nelem * MAPEM.idymax)
82 for (j = i + MAPEM.idx; j < i + GeoPar.nelem; j += MAPEM.idxmax)
85 if (MAPEM.gamma1 > 0.0)
87 MAPEM.qj[j] *= recon[j] * MAPEM.Sjjinv[j] / MAPEM.gamma1;
91 if (MAPEM.Wjj[j] > Consts.zero)
92 MAPEM.qj[j] *= recon[j] / MAPEM.Wjj[j];
94 MAPEM.qj[j] = recon[j]; //keep recon[j] unchanged when qj[j]is zero, by wei 1/2005
101 for (i = MAPEM.idy * Blob.M; i < Blob.M * (Blob.H - 1) + 1;
102 i += Blob.M * MAPEM.idymax)
104 for (j = i + MAPEM.idx; j < i + Blob.M; j += MAPEM.idxmax)
107 if (j % 2 == Blob.pr)
109 if (MAPEM.Wjj[j] > Consts.zero)
110 MAPEM.qj[j] *= recon[j] / MAPEM.Wjj[j];
112 MAPEM.qj[j] = recon[j]; //recon[j] keep unchanged when Wjj[j] is zero, by wei 1/2005
118 if (MAPEM.gamma1 > 0.0)
119 { // blob case will not enter here.hstau
120 //====== step 2: compute "smoothing component" from "recon" to "Sx"
121 applyK(recon, MAPEM.Kx, GeoPar.nelem);
123 // ---Set edges to 0.0
124 for (i = 0; i < GeoPar.nelem; i++)
126 MAPEM.Kx[i] = 0.0; // first row
127 //MAPEM.Kx[(i-1)*GeoPar.nelem+1] = 0.0;
128 MAPEM.Kx[(i + 1) * GeoPar.nelem - 1] = 0.0; // last column
129 MAPEM.Kx[i * GeoPar.nelem] = 0.0; // first column
130 MAPEM.Kx[(GeoPar.nelem - 1) * GeoPar.nelem + i] = 0.0; // last row
133 applyK(MAPEM.Kx, MAPEM.Sx, GeoPar.nelem);
135 //====== step 3: compute x for iteration k+1 from "Sx" and "qj"
137 for (i = MAPEM.idy * GeoPar.nelem;
138 i < GeoPar.nelem * (GeoPar.nelem - 1) + 1;
139 i += GeoPar.nelem * MAPEM.idymax)
141 for (j = i + MAPEM.idx; j < i + GeoPar.nelem; j += MAPEM.idxmax)
143 p = MAPEM.Wjj[j] / (REAL) 9.0 - recon[j]
144 + MAPEM.Sx[j] * MAPEM.Sjjinv[j] / (REAL) 9.0;
145 q = MAPEM.qj[j] / (REAL) 9.0;
146 if ((p <= 0.0) || (q / (p * p) >= 0.005))
148 MAPEM.qj[j] = (REAL) 0.5
149 * (-p + (REAL) sqrt(p * p + 4 * q));
153 MAPEM.qj[j] = q / p * (1 - (REAL) 0.5 * (q / (p * p)));
157 } // not with the blob case .hstau
159 //====== step 4: replace k-th iteration by (k+1)-th iteration in the
163 for (i = MAPEM.idy * GeoPar.nelem;
164 i < GeoPar.nelem * (GeoPar.nelem - 1) + 1;
165 i += GeoPar.nelem * MAPEM.idymax)
167 for (j = i + MAPEM.idx; j < i + GeoPar.nelem; j += MAPEM.idxmax)
169 if (MAPEM.qj[j] < Consts.zero)
171 recon[j] = Consts.zero;
175 recon[j] = MAPEM.qj[j];
183 for (i = MAPEM.idy * Blob.M; i < Blob.M * (Blob.H - 1) + 1;
184 i += Blob.M * MAPEM.idymax)
186 for (j = i + MAPEM.idx; j < i + Blob.M; j += MAPEM.idxmax)
189 if (j % 2 == Blob.pr)
190 { //only the hex grid points
191 if (MAPEM.qj[j] < Consts.zero)
193 recon[j] = Consts.zero;
197 recon[j] = MAPEM.qj[j];
208 c======================================================================
209 c source & dest are nelem*nelem arrays. source contains a picture
210 c to be convolved with the 3*3 kernel whose central element is 1
211 c and the rest are -1/8. The result goes in dest.
212 c======================================================================
215 void emap_class::applyK(REAL* source, REAL* dest, INTEGER nelem)
219 for (i = 0; i < nelem * nelem; i++)
224 for (i = 0; i < nelem - 1; i++)
226 for (j = 0; j < nelem; j++)
228 dest[i + (j) * nelem] += source[i + (j) * nelem + 1];
233 for (i = 0; i < nelem - 1; i++)
235 for (j = 0; j < nelem; j++)
237 dest[i + (j) * nelem + 1] += source[i + (j) * nelem];
242 for (i = 0; i < nelem; i++)
244 for (j = 0; j < nelem - 1; j++)
246 dest[i + (j) * nelem] += source[i + (j) * nelem + nelem];
250 for (i = 0; i < nelem; i++)
252 for (j = 0; j < nelem - 1; j++)
255 dest[i + (j) * nelem + nelem] += source[i + (j) * nelem];
260 for (i = 0; i < nelem - 1; i++)
262 for (j = 0; j < nelem - 1; j++)
264 dest[i + (j) * nelem] += source[i + (j) * nelem + nelem + 1];
268 for (i = 0; i < nelem - 1; i++)
270 for (j = 0; j < nelem - 1; j++)
272 dest[i + (j) * nelem + nelem + 1] += source[i + (j) * nelem];
277 for (i = 0; i < nelem - 1; i++)
279 for (j = 0; j < nelem - 1; j++)
281 dest[i + (j) * nelem + nelem] += source[i + (j) * nelem + 1];
286 for (i = 0; i < nelem - 1; i++)
288 for (j = 0; j < nelem - 1; j++)
290 dest[i + (j) * nelem + 1] += source[i + (j) * nelem + nelem];
295 sscal(nelem * nelem, -1 / 8., dest, 1);
297 for (i = 0; i < nelem * nelem; i++)
299 dest[i] += source[i];