src/snark/bdhk.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/bdhk.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9  */
  10
  11 #include <cstdio>
  12 #include <cmath>
  13
  14 #include "blkdta.h"
  15 #include "blob.h"
  16 #include "geom.h"
  17 #include "fourie.h"
  18 #include "modefl.h"
  19 #include "consts.h"
  20 #include "uiod.h"
  21 #include "anglst.h"
  22 #include "errfac.h"
  23 #include "errpar.h"
  24 #include "pick.h"
  25 #include "pseudo.h"
  26 #include <cstdio>
  27 #include <iostream>
  28 #include "infile.h"
  29 #include "anglst.h"
  30 #include <iomanip>
  31 #include <fstream>
  32 #include "second.h"
  33 #include "wray.h"
  34 #include "raysel.h"
  35 #include "bdhk.h"
  36 using namespace std;
  37
  38 INTEGER bdhk_class::Init()
  39 {
  40
  41         flag_kount = TRUE;
  42
  43         if (RaySel.useray)
  44                 kount = GeoPar.prjnum * GeoPar.usrays;
  45         else
  46                 kount = GeoPar.prjnum * GeoPar.snrays;
  47         //finishNow=FALSE;
  48         arrayRays = new struct NPNR[kount];
  49         orderRays(arrayRays, kount);
  50
  51         //Lastrecon = new REAL[GeoPar.area]; // the array to keep the recon of the least TV seen so far
  52         //total_time=0;
  53         return 0;
  54 }
  55
  56 INTEGER bdhk_class::Reset()
  57 {
  58         delete[] arrayRays;
  59         //delete [] Lastrecon;
  60         return 0;
  61 }
  62
  63 BOOLEAN bdhk_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
  64 {
  65
  66         INTEGER count = 0;
  67
  68         while (flag_kount)
  69         {
  70
  71                 REAL resX = 0;
  72                 REAL* v = new REAL[GeoPar.area];
  73                 REAL* d = new REAL[GeoPar.area];
  74
  75                 if (iter == 1 && count == 0)
  76                 {
  77                         BOOLEAN eol;
  78                         REAL word;
  79
  80                         //checking for BETA
  81                         word = InFile.getnum(TRUE, &eol);
  82                         if (eol)
  83                         {
  84                                 fprintf(output,
  85                                                 "\n          **** error - you must specify 1 parameter list ***");
  86                                 return TRUE;
  87                         }
  88                         BETA = word;
  89                         if (BETA <= Consts.zero)
  90                         {
  91                                 fprintf(output,
  92                                                 "\n          **** error - BETA = %10.6f must be greater than %12.4e***",
  93                                                 word, Consts.zero);
  94                                 return TRUE;
  95                         }
  96                         else
  97                                 fprintf(output, "\n BETA       = %10.6f", word);
  98
  99                         //checking for kount
 100                         word = InFile.getnum(FALSE, &eol);
 101                         if (eol)
 102                         {
 103                                 fprintf(output, "\n          **** error - parameter list***");
 104                                 return TRUE;
 105                         }
 106
 107                         if (word <= 0)
 108                         {
 109                                 if (RaySel.useray)
 110                                         KOUNT = GeoPar.prjnum * GeoPar.usrays;
 111                                 else
 112                                         KOUNT = GeoPar.prjnum * GeoPar.snrays;
 113                                 fprintf(output,
 114                                                 "\n        KOUNT remained in default, KOUNT = %d                                       ***",
 115                                                 KOUNT);
 116
 117                         }
 118                         else
 119                         {
 120                                 KOUNT = word;
 121                                 fprintf(output, "\n KOUNT        = %10.6f", word);
 122                         }
 123
 124                 }  //end if (iter==1 && count==0)
 125
 126                 resX = calculateResidual(recon, list, weight, arrayRays, kount);
 127
 128                 BETA *= 2;
 129
 130                 REAL* grad = new REAL[GeoPar.area];
 131                 REAL* reconBU = new REAL[GeoPar.area];
 132
 133                 BOOLEAN flag = calculateGradientTV(grad, recon);
 134
 135                 if (flag == FALSE)
 136                 {
 137                         //initialiaze the d vector to all zero
 138                         for (int iter1 = 0; iter1 < GeoPar.area; iter1++)
 139                                 d[iter1] = 0;
 140                 }
 141                 else
 142                 {
 143                         //have d vector equal to g/||g||
 144                         REAL norm_grad = calcNorm(grad);
 145
 146                         if (trace == 1)
 147                                 cout << endl << "norm_grad=" << norm_grad << "    " << endl
 148                                                 << endl;
 149
 150                         for (int iter1 = 0; iter1 < GeoPar.area; iter1++)
 151                                 d[iter1] = grad[iter1] / norm_grad;
 152
 153                 }
 154
 155                 REAL resAv;
 156                 REAL TVx = calcTV(recon);
 157
 158                 if (trace == 1)
 159                         cout << endl << endl << "TVx=" << TVx << "    ";
 160
 161                 REAL TVv;
 162                 BOOLEAN logic = TRUE;
 163                 do
 164                 {
 165                         if (trace == 1)
 166                                 cout << endl << "entered DO-WHILE" << endl;
 167
 168                         resAv = resX + 1;
 169                         BETA /= 2;
 170
 171                         if (trace == 1)
 172                                 cout << endl << "BETA=" << BETA << "    " << endl;
 173
 174                         for (int iter1 = 0; iter1 < GeoPar.area; iter1++)
 175                                 v[iter1] = recon[iter1] - BETA * d[iter1];
 176
 177                         TVv = calcTV(v);
 178
 179                         if (trace == 1)
 180                                 cout << "TVv=" << TVv << "    " << endl;
 181 //*************************************************************************
 182 //CHECKING THE CONDITION:
 183 //~~~~~~~~~~~~~~~~~~~~~~~
 184                         if (TVv <= TVx)
 185                         {
 186
 187                                 //Backing Up recon array as v is going to be passed to MyART
 188                                 for (int iter1 = 0; iter1 < GeoPar.area; iter1++)
 189                                 {
 190                                         reconBU[iter1] = recon[iter1];
 191                                         recon[iter1] = v[iter1];
 192                                 }  //now recon contains the vector v and will be passed to MyART
 193
 194                                 if (trace == 1)
 195                                         cout << endl << "resX=" << resX << "   ";
 196
 197                                 MyART(recon, list, weight, arrayRays, kount);
 198                                 count += kount;
 199
 200                                 resAv = calculateResidual(recon, list, weight, arrayRays,
 201                                                 kount);
 202
 203                                 //if (trace == 1)
 204                                 //cout<<endl<<"resX="<<resX<<"   ";
 205                                 if (trace == 1)
 206                                         cout << "resAv=" << resAv << endl << endl;
 207
 208                                 if (resAv >= resX) //only if the condition in the while loop kept as TRUE, we let recon to go back to its previous state. Otherwise stays as X_k+1
 209                                 {
 210                                         for (int iter1 = 0; iter1 < GeoPar.area; iter1++)
 211                                                 recon[iter1] = reconBU[iter1];
 212                                 }
 213                                 else
 214                                         logic = FALSE;
 215
 216                         } // end if (TVv <= TVx)
 217 //**************************************************************************
 218                 } while (logic);
 219
 220                 if (count >= KOUNT)
 221                         flag_kount = FALSE;
 222
 223 //deallocation of memory:
 224                 delete[] d;
 225                 delete[] v;
 226                 delete[] grad;
 227                 delete[] reconBU;
 228
 229                 if (trace == 1)
 230                         cout << "iter " << iter << " ENDED with BETA=" << BETA << endl
 231                                         << endl;
 232
 233         } //end while (flag_kount)
 234
 235 //returning snark iteration
 236         flag_kount = TRUE;
 237         return FALSE;
 238
 239 } //end of snark iteration
 240