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/trm5.cpp $
5 $LastChangedRevision: 154 $
6 $Date: 2014-07-28 11:32:36 -0400 (Mon, 28 Jul 2014) $
8 ***********************************************************
10 Kullback-Leibler value based stopping criterion
29 void trm5_class::Init()
31 // read desired epsilon from command
33 epsilon = InFile.getnum(FALSE, &eol);
35 if (epsilon <= Consts_class::zero)
37 cout << "\n **** epsilon must be specified and greater than ZERO";
38 cout << "\n **** program aborted\n";
42 // read reporting command (disabled by default)
43 INTEGER rprt = CHAR2INT('r', 'p', 'r', 't');
44 if (InFile.getwrd(FALSE, &eol, &rprt, 1) == rprt)
47 cout << "\n reporting is enabled";
48 cout << "\n reporting file: RPRTklds";
50 skips = InFile.getnum(FALSE, &eol);
53 cout << "\n report skip factor " << skips;
57 cout << "\n reporting on every iteration";
63 cout << "\n reporting is disabled";
66 cout << "\n epsilon = " << epsilon;
69 BOOLEAN trm5_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
72 // if currently blob, pixelize it
75 reconp = new REAL[GeoPar.area];
76 Blob.blob2pix(recon, reconp);
79 // variables needed during computation of the current epsilon
86 // compute current KL distance
87 REAL currentEpsilon = 0;
88 for (INTEGER i = 0; i < GeoPar.prjnum && continu; i++)
90 for (INTEGER j = GeoPar.fusray; j <= GeoPar.lusray; j++)
92 raysum = Anglst.prdta(i, j);
93 psum = pseudo(reconp, i, j, list, weight, &numb, &snorm, TRUE, TRUE);
99 // report on non-applicability and abort
100 std::cout << "\n****Error: Negative raysum in prjfil detected, KLDS is not applicable! Terminating algorithm execution.";
103 else if (!isnan(currentEpsilon) && currentEpsilon!=std::numeric_limits<double>::infinity() && currentEpsilon!=std::numeric_limits<double>::max())
105 if (raysum<=Consts.zero && psum<=Consts.zero)
109 else if (raysum<=Consts.zero)
111 currentEpsilon += psum;
113 else if (psum<=Consts.zero)
115 // set to infinity and terminate
116 if (std::numeric_limits<double>::has_infinity) currentEpsilon = std::numeric_limits<double>::infinity();
117 else currentEpsilon = std::numeric_limits<double>::max();
124 if ((raysum / psum) > Consts.zero) currentEpsilon += (raysum * log(raysum / psum));
125 else currentEpsilon += (raysum * log(Consts.zero));
127 currentEpsilon += psum - raysum;
138 if (currentEpsilon <= epsilon)
143 cout << "\n current epsilon (KL distance) = " << currentEpsilon;
148 RPRTfile = fopen("RPRTklds", "w");
149 fprintf(RPRTfile,"Kullback-Leibler distance stopping criterion reporting output");
150 fprintf(RPRTfile,"\n iter Kullback-Leibler distance");
151 fprintf(RPRTfile,"\n %5i %20.9f", iter, currentEpsilon);
155 RPRTfile = fopen("RPRTklds", "a");
156 fprintf(RPRTfile,"\n %5i %20.9f", iter, currentEpsilon);
166 if (report && (iter == 1 || skips == 0 || (iter % (INTEGER) skips) == 0))
168 cout << "\n current epsilon (KL distance) = " << currentEpsilon;
173 RPRTfile = fopen("RPRTklds", "w");
174 fprintf(RPRTfile,"Kullback-Leibler distance stopping criterion reporting output");
175 fprintf(RPRTfile,"\n iter Kullback-Leibler distance");
176 fprintf(RPRTfile,"\n %5i %20.9f", iter, currentEpsilon);
180 RPRTfile = fopen("RPRTklds", "a");
181 fprintf(RPRTfile,"\n %5i %20.9f", iter, currentEpsilon);