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/trm4.cpp $
5 $LastChangedRevision: 131 $
6 $Date: 2014-07-14 19:49:33 -0400 (Mon, 14 Jul 2014) $
8 ***********************************************************
10 Implementation of the indicator function I(x) used for MLEM-STOP
29 void trm4_class::Init()
32 // read reporting command (disabled by default)
33 INTEGER rprt = CHAR2INT('r','p','r','t');
34 if (InFile.getwrd(FALSE, &eol, &rprt, 1) == rprt)
37 cout << "\n reporting is enabled";
38 cout << "\n reporting file: RPRTmlst";
40 skips = InFile.getnum(FALSE, &eol);
43 cout << "\n report skip factor " << skips;
47 cout << "\n reporting on every iteration";
53 cout << "\n reporting is disabled";
57 BOOLEAN trm4_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
59 // numerator and denominator of the formula
62 // temporary variables needed during calculation
65 static INTEGER prjnum;
66 // # the bigger value of usrays or snrays
68 // array for projection data
70 // temporary variable for the count of near-simultaneous arrivals of photon pairs for the ith LOR
72 // variables for geometry data
73 static INTEGER*** pixelnumbers;
74 static REAL*** weights;
75 static INTEGER** numpixels;
78 // gather data which only needs to be gathered during the first iteration
81 prjnum = GeoPar.prjnum;
84 // max # pixels that may be intersected by a ray
85 INTEGER pixelcount = GeoPar.nelem * 2;
87 prjdata = new REAL[prjnum * nrays];
88 for (INTEGER i = 0; i < prjnum; i++)
90 ProjFile.ReadProj(i, &prjdata[i * nrays], nrays);
93 // allocate memory and read list of aij for each ray
94 numpixels = new INTEGER*[prjnum];
95 pixelnumbers = new INTEGER**[prjnum];
96 weights = new REAL**[prjnum];
97 for (INTEGER i = 0; i < prjnum; i++)
99 numpixels[i] = new INTEGER[nrays];
100 pixelnumbers[i] = new INTEGER*[nrays];
101 weights[i] = new REAL*[nrays];
102 for (INTEGER j = 0; j < nrays; j++)
104 pixelnumbers[i][j] = new INTEGER[pixelcount];
105 weights[i][j] = new REAL[pixelcount];
109 wray(i, j, pixelnumbers[i][j], weights[i][j], &numpixels[i][j], &snorm);
113 Blob.bwray(i, j, pixelnumbers[i][j], weights[i][j], &numpixels[i][j], &snorm);
119 // calculate numerator and denominator
120 for (INTEGER i = 0; i < prjnum; i++)
122 for (INTEGER j = 0; j < nrays; j++)
125 for (INTEGER k = 0; k < numpixels[i][j]; k++)
127 sumaijxj += (weights[i][j][k]) * recon[pixelnumbers[i][j][k]];
129 denominator += sumaijxj;
131 bi = prjdata[i * nrays + j];
132 numerator += pow(bi - sumaijxj, 2);
136 // calculate I(x) and check if <=1
137 if (numerator / denominator <= 1)
139 // free all allocated memory
141 for (INTEGER i = 0; i < prjnum; i++)
143 for (INTEGER j = 0; j < nrays; j++)
145 delete[] (pixelnumbers[i][j]);
146 delete[] (weights[i][j]);
148 delete[] (numpixels[i]);
149 delete[] (pixelnumbers[i]);
150 delete[] (weights[i]);
152 delete[] (numpixels);
153 delete[] (pixelnumbers);
159 cout << "\n MLEM-STOP indicator function I(x)=" << numerator/denominator;
164 RPRTfile = fopen("RPRTmlst", "w");
165 fprintf(RPRTfile,"MLEM-STOP indicator function I(x) reporting output");
166 fprintf(RPRTfile,"\n iter I(x)");
167 fprintf(RPRTfile,"\n %5i %20.9f", iter, numerator/denominator);
171 RPRTfile = fopen("RPRTmlst", "a");
172 fprintf(RPRTfile,"\n %5i %20.9f", iter, numerator/denominator);
182 if (report && (iter==1 || skips==0 || (iter%(INTEGER)skips)==0))
184 cout << "\n MLEM-STOP indicator function I(x)=" << numerator/denominator;
189 RPRTfile = fopen("RPRTmlst", "w");
190 fprintf(RPRTfile,"MLEM-STOP indicator function I(x) reporting output");
191 fprintf(RPRTfile,"\n iter I(x)");
192 fprintf(RPRTfile,"\n %5i %20.9f", iter, numerator/denominator);
196 RPRTfile = fopen("RPRTmlst", "a");
197 fprintf(RPRTfile,"\n %5i %20.9f", iter, numerator/denominator);