src/snark/trm6.cpp

   1 /*
   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) $
   7  $Author: olangthaler $
   8  ***********************************************************
   9
  10  Residual based stopping criterion
  11
  12  */
  13
  14 #include <cstdio>
  15 #include <cstdlib>
  16 #include <iostream>
  17 #include <math.h>
  18 #include "trm6.h"
  19 #include "anglst.h"
  20 #include "projfile.h"
  21 #include "pseudo.h"
  22 #include "consts.h"
  23 #include "infile.h"
  24 #include "blob.h"
  25
  26 using namespace std;
  27
  28 void trm6_class::Init()
  29 {
  30         // read desired epsilon from command
  31         BOOLEAN eol;
  32         epsilon = InFile.getnum(FALSE, &eol);
  33
  34         if (epsilon < Consts_class::zero)
  35         {
  36                 cout << "\n **** epsilon must be specified and >= ZERO";
  37                 cout << "\n **** program aborted\n";
  38                 exit(-1);
  39         }
  40
  41         // read reporting command (disabled by default)
  42         INTEGER rprt = CHAR2INT('r', 'p', 'r', 't');
  43         if (InFile.getwrd(FALSE, &eol, &rprt, 1) == rprt)
  44         {
  45                 report = true;
  46                 cout << "\n         reporting is enabled";
  47                 cout << "\n         reporting file: RPRTresi";
  48
  49                 skips = InFile.getnum(FALSE, &eol);
  50                 if (skips > 0)
  51                 {
  52                         cout << "\n         report skip factor " << skips;
  53                 }
  54                 else
  55                 {
  56                         cout << "\n         reporting on every iteration";
  57                 }
  58         }
  59         else
  60         {
  61                 report = false;
  62                 cout << "\n         reporting is disabled";
  63         }
  64
  65         cout << "\n         epsilon = " << epsilon;
  66 }
  67
  68 BOOLEAN trm6_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
  69 {
  70         REAL* reconp = recon;
  71         // if currently blob, pixelize it
  72         if (!Blob.pix_basis)
  73         {
  74                 reconp = new REAL[GeoPar.area];
  75                 Blob.blob2pix(recon, reconp);
  76         }
  77
  78         // variables needed during computation of the current residual
  79         INTEGER numb = 0;
  80         REAL snorm = 0;
  81         REAL raysum = 0;
  82         REAL psum = 0;
  83         REAL rdis = 0;
  84         REAL theta = 0;
  85         REAL sinth = 0;
  86         REAL costh = 0;
  87         BOOLEAN linef = true;
  88
  89         // compute current residual
  90         for (int i = 0; i < GeoPar.prjnum; i++)
  91         {
  92                 REAL rinc = GeoPar.pinc;
  93                 Anglst.getang(i, &theta, &sinth, &costh);
  94
  95                 if (GeoPar.vri)
  96                 {
  97                         rinc = GeoPar.pinc * MAX0(fabs(sinth), fabs(costh));
  98                 }
  99
 100                 for (int j = GeoPar.fusray; j <= GeoPar.lusray; j++)
 101                 {
 102                         raysum = Anglst.prdta(i, j);
 103                         psum = pseudo(reconp, i, j, list, weight, &numb, &snorm, linef,
 104                                         TRUE);
 105
 106                         if (numb != 0)
 107                         {
 108                                 if (linef && !GeoPar.line)
 109                                         raysum /= rinc;
 110                                 if (!linef && GeoPar.line)
 111                                         raysum *= rinc;
 112                                 rdis += pow((psum - raysum), 2);
 113                         }
 114                 }
 115         }
 116         REAL currentEpsilon = sqrt(rdis);
 117         if (!Blob.pix_basis)
 118         {
 119                 delete[] reconp;
 120         }
 121
 122         if (currentEpsilon <= epsilon)
 123         {
 124                 // additional output
 125                 if (report)
 126                 {
 127                         cout << "\n         current epsilon (residual) = " << currentEpsilon;
 128
 129                         FILE* RPRTfile;
 130                         if (iter==1)
 131                         {
 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);
 136                         }
 137                         else
 138                         {
 139                                 RPRTfile = fopen("RPRTresi", "a");
 140                                 fprintf(RPRTfile,"\n %5i  %20.9f", iter, currentEpsilon);
 141                         }
 142                         fclose(RPRTfile);
 143                 }
 144
 145                 return true;
 146         }
 147         else
 148         {
 149                 // additional output
 150                 if (report && (iter == 1 || skips == 0 || (iter % (INTEGER) skips) == 0))
 151                 {
 152                         cout << "\n         current epsilon (residual) = " << currentEpsilon;
 153
 154                         FILE* RPRTfile;
 155                         if (iter==1)
 156                         {
 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);
 161                         }
 162                         else
 163                         {
 164                                 RPRTfile = fopen("RPRTresi", "a");
 165                                 fprintf(RPRTfile,"\n %5i  %20.9f", iter, currentEpsilon);
 166                         }
 167                         fclose(RPRTfile);
 168                 }
 169
 170                 return false;
 171         }
 172 }