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/trm6.cpp $
5 $LastChangedRevision: 145 $
6 $Date: 2014-07-17 19:49:25 -0400 (Thu, 17 Jul 2014) $
8 ***********************************************************
10 Residual based stopping criterion
28 void trm6_class::Init()
30 // read desired epsilon from command
32 epsilon = InFile.getnum(FALSE, &eol);
34 if (epsilon < Consts_class::zero)
36 cout << "\n **** epsilon must be specified and >= ZERO";
37 cout << "\n **** program aborted\n";
41 // read reporting command (disabled by default)
42 INTEGER rprt = CHAR2INT('r', 'p', 'r', 't');
43 if (InFile.getwrd(FALSE, &eol, &rprt, 1) == rprt)
46 cout << "\n reporting is enabled";
47 cout << "\n reporting file: RPRTresi";
49 skips = InFile.getnum(FALSE, &eol);
52 cout << "\n report skip factor " << skips;
56 cout << "\n reporting on every iteration";
62 cout << "\n reporting is disabled";
65 cout << "\n epsilon = " << epsilon;
68 BOOLEAN trm6_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
71 // if currently blob, pixelize it
74 reconp = new REAL[GeoPar.area];
75 Blob.blob2pix(recon, reconp);
78 // variables needed during computation of the current residual
89 // compute current residual
90 for (int i = 0; i < GeoPar.prjnum; i++)
92 REAL rinc = GeoPar.pinc;
93 Anglst.getang(i, &theta, &sinth, &costh);
97 rinc = GeoPar.pinc * MAX0(fabs(sinth), fabs(costh));
100 for (int j = GeoPar.fusray; j <= GeoPar.lusray; j++)
102 raysum = Anglst.prdta(i, j);
103 psum = pseudo(reconp, i, j, list, weight, &numb, &snorm, linef,
108 if (linef && !GeoPar.line)
110 if (!linef && GeoPar.line)
112 rdis += pow((psum - raysum), 2);
116 REAL currentEpsilon = sqrt(rdis);
122 if (currentEpsilon <= epsilon)
127 cout << "\n current epsilon (residual) = " << currentEpsilon;
132 RPRTfile = fopen("RPRTresi", "w");
133 fprintf(RPRTfile,"Residual based stopping criterion reporting output");
134 fprintf(RPRTfile,"\n iter Residual");
135 fprintf(RPRTfile,"\n %5i %20.9f", iter, currentEpsilon);
139 RPRTfile = fopen("RPRTresi", "a");
140 fprintf(RPRTfile,"\n %5i %20.9f", iter, currentEpsilon);
150 if (report && (iter == 1 || skips == 0 || (iter % (INTEGER) skips) == 0))
152 cout << "\n current epsilon (residual) = " << currentEpsilon;
157 RPRTfile = fopen("RPRTresi", "w");
158 fprintf(RPRTfile,"Residual based stopping criterion reporting output");
159 fprintf(RPRTfile,"\n iter Residual");
160 fprintf(RPRTfile,"\n %5i %20.9f", iter, currentEpsilon);
164 RPRTfile = fopen("RPRTresi", "a");
165 fprintf(RPRTfile,"\n %5i %20.9f", iter, currentEpsilon);