src/snark/sart.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/sart.cpp $
   5  $LastChangedRevision: 95 $
   6  $Date: 2014-07-02 19:43:36 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9  */
  10
  11 #include <cstdio>
  12 #include <iostream>
  13 #include "sart.h"
  14 #include "anglst.h"
  15 #include "pseudo.h"
  16 #include "GeometricBehaviour.h"
  17 #include "SARTConfig.h"
  18
  19 using namespace std;
  20
  21 sart::~sart() {
  22         this->Reset();
  23 }
  24
  25 GeometricBehaviour* sart::geometricBehaviour = NULL;
  26 REAL sart::relaxationParameter = 1.0;
  27 REAL sart::zeroThreshold = 0.0;
  28
  29 INTEGER sart::Init() {
  30         this->geometricBehaviour = GeometricBehaviour::getInstance();
  31         return 0;
  32 }
  33
  34 INTEGER sart::Reset() {
  35         if (this->geometricBehaviour != NULL) {
  36                 delete this->geometricBehaviour;
  37                 this->geometricBehaviour = NULL;
  38         }
  39
  40         this->relaxationParameter = 1.0;
  41         this->zeroThreshold = 0.0;
  42
  43         return 0;
  44 }
  45
  46 BOOLEAN sart::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter) {
  47
  48         // read configuration on first iteration
  49         if (iter == 1) {
  50                 SARTConfig::readSARTConfiguration();
  51                 this->relaxationParameter = SARTConfig::relaxiation;
  52                 this->zeroThreshold = SARTConfig::zeroThreshold;
  53         }
  54
  55         INTEGER numberOfPixels;
  56         REAL snorm = 0;
  57
  58         REAL raySumMeasured;
  59
  60         REAL currentRaySum;
  61
  62         REAL rayWeightSum;
  63         REAL* pixelRayWeight = new REAL[this->geometricBehaviour->getArea()];
  64         REAL* pixelRayWeightSum = new REAL[this->geometricBehaviour->getArea()];
  65
  66         REAL* correction = new REAL[this->geometricBehaviour->getArea()];
  67
  68         // initialize values which are calculated once per iteration
  69         for (int i = 0; i < this->geometricBehaviour->getArea(); i++) {
  70                 pixelRayWeight[i] = 0.0;
  71                 pixelRayWeightSum[i] = 0.0;
  72                 correction[i] = 0.0;
  73         }
  74
  75         // loop through all rays
  76         for (int projectionNr = 0; projectionNr < GeoPar.prjnum; projectionNr++) {
  77                 for (int rayNr = 0; rayNr < GeoPar.nrays; rayNr++) {
  78
  79                         raySumMeasured = Anglst.prdta(projectionNr, rayNr);
  80
  81                         this->geometricBehaviour->getRayTrace(projectionNr, rayNr, list,
  82                                         weight, &numberOfPixels, &snorm);
  83
  84                         // calculate current raysum
  85                         currentRaySum = 0.0;
  86                         rayWeightSum = 0.0;
  87
  88                         for (int i = 0; i < numberOfPixels; i++) {
  89                                 currentRaySum += recon[list[i]] * weight[i];
  90                                 rayWeightSum += weight[i];
  91                                 pixelRayWeight[list[i]] = weight[i];
  92                                 pixelRayWeightSum[list[i]] += weight[i];
  93                         }
  94
  95                         if (rayWeightSum < this->zeroThreshold) {
  96                                 continue;
  97                         }
  98
  99                         for (int i = 0; i < numberOfPixels; i++) {
 100                                 correction[list[i]] += pixelRayWeight[list[i]]
 101                                                 * (raySumMeasured - currentRaySum) / rayWeightSum;
 102                                 pixelRayWeight[list[i]] = 0.0;
 103                         }
 104                 }
 105         }
 106
 107         // update pixel
 108         for (int i = 0; i < this->geometricBehaviour->getArea(); i++) {
 109
 110                 if (pixelRayWeightSum[i] < this->zeroThreshold) {
 111                         continue;
 112                 }
 113
 114                 REAL oldRecon = recon[i];
 115                 recon[i] += this->relaxationParameter / pixelRayWeightSum[i]
 116                                 * correction[i];
 117                 if (recon[i] < 0.0) {
 118                         recon[i] = 0.0;
 119                 }
 120         }
 121
 122         delete[] pixelRayWeight;
 123         delete[] pixelRayWeightSum;
 124         delete[] correction;
 125
 126         return false;
 127 }