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/bdhk.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
38 INTEGER bdhk_class::Init()
44 kount = GeoPar.prjnum * GeoPar.usrays;
46 kount = GeoPar.prjnum * GeoPar.snrays;
48 arrayRays = new struct NPNR[kount];
49 orderRays(arrayRays, kount);
51 //Lastrecon = new REAL[GeoPar.area]; // the array to keep the recon of the least TV seen so far
56 INTEGER bdhk_class::Reset()
59 //delete [] Lastrecon;
63 BOOLEAN bdhk_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
72 REAL* v = new REAL[GeoPar.area];
73 REAL* d = new REAL[GeoPar.area];
75 if (iter == 1 && count == 0)
81 word = InFile.getnum(TRUE, &eol);
85 "\n **** error - you must specify 1 parameter list ***");
89 if (BETA <= Consts.zero)
92 "\n **** error - BETA = %10.6f must be greater than %12.4e***",
97 fprintf(output, "\n BETA = %10.6f", word);
100 word = InFile.getnum(FALSE, &eol);
103 fprintf(output, "\n **** error - parameter list***");
110 KOUNT = GeoPar.prjnum * GeoPar.usrays;
112 KOUNT = GeoPar.prjnum * GeoPar.snrays;
114 "\n KOUNT remained in default, KOUNT = %d ***",
121 fprintf(output, "\n KOUNT = %10.6f", word);
124 } //end if (iter==1 && count==0)
126 resX = calculateResidual(recon, list, weight, arrayRays, kount);
130 REAL* grad = new REAL[GeoPar.area];
131 REAL* reconBU = new REAL[GeoPar.area];
133 BOOLEAN flag = calculateGradientTV(grad, recon);
137 //initialiaze the d vector to all zero
138 for (int iter1 = 0; iter1 < GeoPar.area; iter1++)
143 //have d vector equal to g/||g||
144 REAL norm_grad = calcNorm(grad);
147 cout << endl << "norm_grad=" << norm_grad << " " << endl
150 for (int iter1 = 0; iter1 < GeoPar.area; iter1++)
151 d[iter1] = grad[iter1] / norm_grad;
156 REAL TVx = calcTV(recon);
159 cout << endl << endl << "TVx=" << TVx << " ";
162 BOOLEAN logic = TRUE;
166 cout << endl << "entered DO-WHILE" << endl;
172 cout << endl << "BETA=" << BETA << " " << endl;
174 for (int iter1 = 0; iter1 < GeoPar.area; iter1++)
175 v[iter1] = recon[iter1] - BETA * d[iter1];
180 cout << "TVv=" << TVv << " " << endl;
181 //*************************************************************************
182 //CHECKING THE CONDITION:
183 //~~~~~~~~~~~~~~~~~~~~~~~
187 //Backing Up recon array as v is going to be passed to MyART
188 for (int iter1 = 0; iter1 < GeoPar.area; iter1++)
190 reconBU[iter1] = recon[iter1];
191 recon[iter1] = v[iter1];
192 } //now recon contains the vector v and will be passed to MyART
195 cout << endl << "resX=" << resX << " ";
197 MyART(recon, list, weight, arrayRays, kount);
200 resAv = calculateResidual(recon, list, weight, arrayRays,
204 //cout<<endl<<"resX="<<resX<<" ";
206 cout << "resAv=" << resAv << endl << endl;
208 if (resAv >= resX) //only if the condition in the while loop kept as TRUE, we let recon to go back to its previous state. Otherwise stays as X_k+1
210 for (int iter1 = 0; iter1 < GeoPar.area; iter1++)
211 recon[iter1] = reconBU[iter1];
216 } // end if (TVv <= TVx)
217 //**************************************************************************
223 //deallocation of memory:
230 cout << "iter " << iter << " ENDED with BETA=" << BETA << endl
233 } //end while (flag_kount)
235 //returning snark iteration
239 } //end of snark iteration