src/snark/trm5.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/trm5.cpp $
   5  $LastChangedRevision: 154 $
   6  $Date: 2014-07-28 11:32:36 -0400 (Mon, 28 Jul 2014) $
   7  $Author: olangthaler $
   8  ***********************************************************
   9
  10  Kullback-Leibler value based stopping criterion
  11
  12  */
  13
  14 #include <cstdio>
  15 #include <cstdlib>
  16 #include <iostream>
  17 #include <math.h>
  18 #include <limits>
  19 #include "trm5.h"
  20 #include "anglst.h"
  21 #include "projfile.h"
  22 #include "pseudo.h"
  23 #include "consts.h"
  24 #include "infile.h"
  25 #include "blob.h"
  26
  27 using namespace std;
  28
  29 void trm5_class::Init()
  30 {
  31         // read desired epsilon from command
  32         BOOLEAN eol;
  33         epsilon = InFile.getnum(FALSE, &eol);
  34
  35         if (epsilon <= Consts_class::zero)
  36         {
  37                 cout << "\n **** epsilon must be specified and greater than ZERO";
  38                 cout << "\n **** program aborted\n";
  39                 exit(-1);
  40         }
  41
  42         // read reporting command (disabled by default)
  43         INTEGER rprt = CHAR2INT('r', 'p', 'r', 't');
  44         if (InFile.getwrd(FALSE, &eol, &rprt, 1) == rprt)
  45         {
  46                 report = true;
  47                 cout << "\n         reporting is enabled";
  48                 cout << "\n         reporting file: RPRTklds";
  49
  50                 skips = InFile.getnum(FALSE, &eol);
  51                 if (skips > 0)
  52                 {
  53                         cout << "\n         report skip factor " << skips;
  54                 }
  55                 else
  56                 {
  57                         cout << "\n         reporting on every iteration";
  58                 }
  59         }
  60         else
  61         {
  62                 report = false;
  63                 cout << "\n         reporting is disabled";
  64         }
  65
  66         cout << "\n         epsilon = " << epsilon;
  67 }
  68
  69 BOOLEAN trm5_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
  70 {
  71         REAL* reconp = recon;
  72         // if currently blob, pixelize it
  73         if (!Blob.pix_basis)
  74         {
  75                 reconp = new REAL[GeoPar.area];
  76                 Blob.blob2pix(recon, reconp);
  77         }
  78
  79         // variables needed during computation of the current epsilon
  80         INTEGER numb;
  81         REAL snorm;
  82         REAL raysum;
  83         REAL psum;
  84         BOOLEAN continu=true;
  85
  86         // compute current KL distance
  87         REAL currentEpsilon = 0;
  88         for (INTEGER i = 0; i < GeoPar.prjnum && continu; i++)
  89         {
  90                 for (INTEGER j = GeoPar.fusray; j <= GeoPar.lusray; j++)
  91                 {
  92                         raysum = Anglst.prdta(i, j);
  93                         psum = pseudo(reconp, i, j, list, weight, &numb, &snorm, TRUE, TRUE);
  94
  95                         if (numb != 0)
  96                         {
  97                                 if (raysum<0)
  98                                 {
  99                                         // report on non-applicability and abort
 100                                         std::cout << "\n****Error: Negative raysum in prjfil detected, KLDS is not applicable! Terminating algorithm execution.";
 101                                         return true;
 102                                 }
 103                                 else if (!isnan(currentEpsilon) && currentEpsilon!=std::numeric_limits<double>::infinity() && currentEpsilon!=std::numeric_limits<double>::max())
 104                                 {
 105                                         if (raysum<=Consts.zero && psum<=Consts.zero)
 106                                         {
 107                                                 // do nothing
 108                                         }
 109                                         else if (raysum<=Consts.zero)
 110                                         {
 111                                                 currentEpsilon += psum;
 112                                         }
 113                                         else if (psum<=Consts.zero)
 114                                         {
 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();
 118
 119                                                 continu=false;
 120                                                 break;
 121                                         }
 122                                         else
 123                                         {
 124                                                 if ((raysum / psum) > Consts.zero) currentEpsilon += (raysum * log(raysum / psum));
 125                                                 else currentEpsilon += (raysum * log(Consts.zero));
 126
 127                                                 currentEpsilon += psum - raysum;
 128                                         }
 129                                 }
 130                         }
 131                 }
 132         }
 133         if (!Blob.pix_basis)
 134         {
 135                 delete[] reconp;
 136         }
 137
 138         if (currentEpsilon <= epsilon)
 139         {
 140                 // additional output
 141                 if (report)
 142                 {
 143                         cout << "\n         current epsilon (KL distance) = " << currentEpsilon;
 144
 145                         FILE* RPRTfile;
 146                         if (iter==1)
 147                         {
 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);
 152                         }
 153                         else
 154                         {
 155                                 RPRTfile = fopen("RPRTklds", "a");
 156                                 fprintf(RPRTfile,"\n %5i  %20.9f", iter, currentEpsilon);
 157                         }
 158                         fclose(RPRTfile);
 159                 }
 160
 161                 return true;
 162         }
 163         else
 164         {
 165                 // additional output
 166                 if (report && (iter == 1 || skips == 0 || (iter % (INTEGER) skips) == 0))
 167                 {
 168                         cout << "\n         current epsilon (KL distance) = " << currentEpsilon;
 169
 170                         FILE* RPRTfile;
 171                         if (iter==1)
 172                         {
 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);
 177                         }
 178                         else
 179                         {
 180                                 RPRTfile = fopen("RPRTklds", "a");
 181                                 fprintf(RPRTfile,"\n %5i  %20.9f", iter, currentEpsilon);
 182                         }
 183                         fclose(RPRTfile);
 184                 }
 185
 186                 return false;
 187         }
 188 }