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) $
8 ***********************************************************
16 #include "GeometricBehaviour.h"
17 #include "SARTConfig.h"
25 GeometricBehaviour* sart::geometricBehaviour = NULL;
26 REAL sart::relaxationParameter = 1.0;
27 REAL sart::zeroThreshold = 0.0;
29 INTEGER sart::Init() {
30 this->geometricBehaviour = GeometricBehaviour::getInstance();
34 INTEGER sart::Reset() {
35 if (this->geometricBehaviour != NULL) {
36 delete this->geometricBehaviour;
37 this->geometricBehaviour = NULL;
40 this->relaxationParameter = 1.0;
41 this->zeroThreshold = 0.0;
46 BOOLEAN sart::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter) {
48 // read configuration on first iteration
50 SARTConfig::readSARTConfiguration();
51 this->relaxationParameter = SARTConfig::relaxiation;
52 this->zeroThreshold = SARTConfig::zeroThreshold;
55 INTEGER numberOfPixels;
63 REAL* pixelRayWeight = new REAL[this->geometricBehaviour->getArea()];
64 REAL* pixelRayWeightSum = new REAL[this->geometricBehaviour->getArea()];
66 REAL* correction = new REAL[this->geometricBehaviour->getArea()];
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;
75 // loop through all rays
76 for (int projectionNr = 0; projectionNr < GeoPar.prjnum; projectionNr++) {
77 for (int rayNr = 0; rayNr < GeoPar.nrays; rayNr++) {
79 raySumMeasured = Anglst.prdta(projectionNr, rayNr);
81 this->geometricBehaviour->getRayTrace(projectionNr, rayNr, list,
82 weight, &numberOfPixels, &snorm);
84 // calculate current raysum
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];
95 if (rayWeightSum < this->zeroThreshold) {
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;
108 for (int i = 0; i < this->geometricBehaviour->getArea(); i++) {
110 if (pixelRayWeightSum[i] < this->zeroThreshold) {
114 REAL oldRecon = recon[i];
115 recon[i] += this->relaxationParameter / pixelRayWeightSum[i]
117 if (recon[i] < 0.0) {
122 delete[] pixelRayWeight;
123 delete[] pixelRayWeightSum;