src/snark/trm4.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/trm4.cpp $
   5  $LastChangedRevision: 131 $
   6  $Date: 2014-07-14 19:49:33 -0400 (Mon, 14 Jul 2014) $
   7  $Author: olangthaler $
   8  ***********************************************************
   9
  10  Implementation of the indicator function I(x) used for MLEM-STOP
  11
  12 */
  13
  14 #include <cstdio>
  15 #include <cstdlib>
  16 #include <iostream>
  17 #include "blkdta.h"
  18 #include "uiod.h"
  19 #include "trm4.h"
  20 #include "anglst.h"
  21 #include "wray.h"
  22 #include "math.h"
  23 #include "projfile.h"
  24 #include "blob.h"
  25 #include "infile.h"
  26
  27 using namespace std;
  28
  29 void trm4_class::Init()
  30 {
  31         BOOLEAN eol;
  32         // read reporting command (disabled by default)
  33         INTEGER rprt = CHAR2INT('r','p','r','t');
  34         if (InFile.getwrd(FALSE, &eol, &rprt, 1) == rprt)
  35         {
  36                 report = true;
  37                 cout << "\n         reporting is enabled";
  38                 cout << "\n         reporting file: RPRTmlst";
  39
  40                 skips = InFile.getnum(FALSE, &eol);
  41                 if (skips>0)
  42                 {
  43                         cout << "\n         report skip factor " << skips;
  44                 }
  45                 else
  46                 {
  47                         cout << "\n         reporting on every iteration";
  48                 }
  49         }
  50         else
  51         {
  52                 report = false;
  53                 cout << "\n         reporting is disabled";
  54         }
  55 }
  56
  57 BOOLEAN trm4_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
  58 {
  59         // numerator and denominator of the formula
  60         REAL numerator = 0;
  61         REAL denominator = 0;
  62         // temporary variables needed during calculation
  63         REAL sumaijxj = 0;
  64         // # projections
  65         static INTEGER prjnum;
  66         // # the bigger value of usrays or snrays
  67         static INTEGER nrays;
  68         // array for projection data
  69         static REAL* prjdata;
  70         // temporary variable for the count of near-simultaneous arrivals of photon pairs for the ith LOR
  71         REAL bi = 0;
  72         // variables for geometry data
  73         static INTEGER*** pixelnumbers;
  74         static REAL*** weights;
  75         static INTEGER** numpixels;
  76         REAL snorm = 0;
  77
  78         // gather data which only needs to be gathered during the first iteration
  79         if (iter == 1)
  80         {
  81                 prjnum = GeoPar.prjnum;
  82                 nrays = GeoPar.nrays;
  83
  84                 // max # pixels that may be intersected by a ray
  85                 INTEGER pixelcount = GeoPar.nelem * 2;
  86
  87                 prjdata = new REAL[prjnum * nrays];
  88                 for (INTEGER i = 0; i < prjnum; i++)
  89                 {
  90                         ProjFile.ReadProj(i, &prjdata[i * nrays], nrays);
  91                 }
  92
  93                 // allocate memory and read list of aij for each ray
  94                 numpixels = new INTEGER*[prjnum];
  95                 pixelnumbers = new INTEGER**[prjnum];
  96                 weights = new REAL**[prjnum];
  97                 for (INTEGER i = 0; i < prjnum; i++)
  98                 {
  99                         numpixels[i] = new INTEGER[nrays];
 100                         pixelnumbers[i] = new INTEGER*[nrays];
 101                         weights[i] = new REAL*[nrays];
 102                         for (INTEGER j = 0; j < nrays; j++)
 103                         {
 104                                 pixelnumbers[i][j] = new INTEGER[pixelcount];
 105                                 weights[i][j] = new REAL[pixelcount];
 106
 107                             if(Blob.pix_basis)
 108                             {
 109                                 wray(i, j, pixelnumbers[i][j], weights[i][j], &numpixels[i][j], &snorm);
 110                             }
 111                             else
 112                             {
 113                                 Blob.bwray(i, j, pixelnumbers[i][j], weights[i][j], &numpixels[i][j], &snorm);
 114                             }
 115                         }
 116                 }
 117         }
 118
 119         // calculate numerator and denominator
 120         for (INTEGER i = 0; i < prjnum; i++)
 121         {
 122                 for (INTEGER j = 0; j < nrays; j++)
 123                 {
 124                         sumaijxj = 0;
 125                         for (INTEGER k = 0; k < numpixels[i][j]; k++)
 126                         {
 127                                 sumaijxj += (weights[i][j][k]) * recon[pixelnumbers[i][j][k]];
 128                         }
 129                         denominator += sumaijxj;
 130
 131                         bi = prjdata[i * nrays + j];
 132                         numerator += pow(bi - sumaijxj, 2);
 133                 }
 134         }
 135
 136         // calculate I(x) and check if <=1
 137         if (numerator / denominator <= 1)
 138         {
 139                 // free all allocated memory
 140                 delete[] (prjdata);
 141                 for (INTEGER i = 0; i < prjnum; i++)
 142                 {
 143                         for (INTEGER j = 0; j < nrays; j++)
 144                         {
 145                                 delete[] (pixelnumbers[i][j]);
 146                                 delete[] (weights[i][j]);
 147                         }
 148                         delete[] (numpixels[i]);
 149                         delete[] (pixelnumbers[i]);
 150                         delete[] (weights[i]);
 151                 }
 152                 delete[] (numpixels);
 153                 delete[] (pixelnumbers);
 154                 delete[] (weights);
 155
 156                 // additional output
 157                 if (report)
 158                 {
 159                         cout << "\n         MLEM-STOP indicator function I(x)=" << numerator/denominator;
 160
 161                         FILE* RPRTfile;
 162                         if (iter==1)
 163                         {
 164                                 RPRTfile = fopen("RPRTmlst", "w");
 165                                 fprintf(RPRTfile,"MLEM-STOP indicator function I(x) reporting output");
 166                                 fprintf(RPRTfile,"\n    iter         I(x)");
 167                                 fprintf(RPRTfile,"\n %5i  %20.9f", iter, numerator/denominator);
 168                         }
 169                         else
 170                         {
 171                                 RPRTfile = fopen("RPRTmlst", "a");
 172                                 fprintf(RPRTfile,"\n %5i  %20.9f", iter, numerator/denominator);
 173                         }
 174                         fclose(RPRTfile);
 175                 }
 176
 177                 return TRUE;
 178         }
 179         else
 180         {
 181                 // additional output
 182                 if (report && (iter==1 || skips==0 || (iter%(INTEGER)skips)==0))
 183                 {
 184                         cout << "\n         MLEM-STOP indicator function I(x)=" << numerator/denominator;
 185
 186                         FILE* RPRTfile;
 187                         if (iter==1)
 188                         {
 189                                 RPRTfile = fopen("RPRTmlst", "w");
 190                                 fprintf(RPRTfile,"MLEM-STOP indicator function I(x) reporting output");
 191                                 fprintf(RPRTfile,"\n    iter         I(x)");
 192                                 fprintf(RPRTfile,"\n %5i  %20.9f", iter, numerator/denominator);
 193                         }
 194                         else
 195                         {
 196                                 RPRTfile = fopen("RPRTmlst", "a");
 197                                 fprintf(RPRTfile,"\n %5i  %20.9f", iter, numerator/denominator);
 198                         }
 199                         fclose(RPRTfile);
 200                 }
 201
 202                 return FALSE;
 203         }
 204 }