src/snark/noise.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/noise.cpp $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9  */
  10
  11 #include <cmath>
  12 #include <cstdio>
  13 #include <cstdlib>
  14 #include <DIGGauss.h>
  15 #include <DIGPoisson.h>
  16 #include "consts.h"
  17 #include "geom.h"
  18 #include "noise.h"
  19 #include "spctrm.h"
  20 #include "uiod.h"
  21
  22 // bug221 - moved noise code to this class - swr - 03/03/07
  23
  24 Noise_class NoisePar = Noise_class();
  25
  26 Noise_class::Noise_class()
  27 {
  28         quanin = 0;
  29         quanmn = 1.0;
  30         quancm = 1.0;
  31         sctnfl = false;
  32         sctnpk = 0.0;
  33         sctnwd = 0.0;
  34         ultnfl = false;
  35         ultnmn = 0.0;
  36         ultnsd = 0.0;
  37         addnfl = false;
  38         addnmn = 0.0;
  39         addnsd = 0.0;
  40         oldnp = -1;
  41         sctcnt = NULL;
  42         sctwts = NULL;
  43         calray = NULL;
  44 }
  45
  46 void Noise_class::resetNoise()
  47 {
  48         // jklukowska, bug 263, in case the noise was set
  49         // previously, we need to reset it to initial values
  50         quanin = 0;
  51         quanmn = 1.0;
  52         quancm = 1.0;
  53         sctnfl = false;
  54         sctnpk = 0.0;
  55         sctnwd = 0.0;
  56         ultnfl = false;
  57         ultnmn = 0.0;
  58         ultnsd = 0.0;
  59         addnfl = false;
  60         addnmn = 0.0;
  61         addnsd = 0.0;
  62         oldnp = -1;
  63         sctcnt = NULL;
  64         sctwts = NULL;
  65         calray = NULL;
  66 }
  67
  68 Noise_class::~Noise_class()
  69 {
  70         if (sctcnt != NULL)
  71                 delete[] sctcnt;
  72         if (sctwts != NULL)
  73                 delete[] sctwts;
  74         if (calray != NULL)
  75                 delete[] calray;
  76 }
  77
  78 void Noise_class::setAdditiveNoise(REAL mean, REAL std)
  79 {
  80         addnfl = true;
  81         addnmn = mean;
  82         addnsd = std;
  83 }
  84
  85 REAL Noise_class::additiveNoise(REAL value)
  86 {
  87         return addnfl ? value + Gauss(addnmn, addnsd) : value;
  88 }
  89
  90 void Noise_class::setMultiplicativeNoise(REAL mean, REAL std)
  91 {
  92         ultnfl = true;
  93         ultnmn = mean;
  94         ultnsd = std;
  95 }
  96
  97 REAL Noise_class::multiplicativeNoise(REAL value)
  98 {
  99         return ultnfl ? value * Gauss(ultnmn, ultnsd) : value;
 100 }
 101
 102 void Noise_class::setScatter(REAL peak, REAL width)
 103 {
 104         sctnfl = true;
 105         sctnpk = peak;
 106         sctnwd = width;
 107         // PRE-COMPUTE THE WEIGHTS FOR SCATTER CALCULATION
 108         // FIRST CONPUTE THE EXTENT OF SCATTER IN TERMS OF NUMBER OF RAYS
 109         // WE NEED SPACE = USRAYS LONG FOR CONVOLUTING
 110
 111         sctcnt = new REAL[GeoPar.usrays];
 112         nsctr = (int) (width / GeoPar.pinc + 1.0);
 113
 114         sctwts = new REAL[nsctr];
 115         for (INTEGER n = 0; n < nsctr; n++)
 116         {
 117                 sctwts[n] = (REAL) (peak * (1.0 - n * GeoPar.pinc / width));
 118         }
 119         sctwts[0] += 1.0;
 120 }
 121
 122 void Noise_class::scatterNoise(REAL* projection)
 123 {
 124         if (sctnfl)
 125         {
 126
 127                 // CONVOLUTE WITH THE SCATTER FUNCTION
 128
 129                 for (INTEGER nr = 0; nr < GeoPar.usrays; nr++)
 130                 {
 131                         sctcnt[nr] = projection[nr];
 132                 }
 133
 134                 for (INTEGER nr = 0; nr < GeoPar.usrays; nr++)
 135                 {
 136                         REAL sum = sctcnt[nr] * sctwts[0];
 137                         REAL vsum = sctwts[0];
 138
 139                         for (INTEGER nrsct = 1; nrsct < nsctr; nrsct++)
 140                         {
 141                                 INTEGER nshift = nrsct;
 142
 143                                 // CONVOLVE CAREFULLY
 144
 145                                 if ((nr - nshift) >= 0)
 146                                 {
 147                                         vsum += sctwts[nrsct];
 148                                         sum += sctcnt[nr - nshift] * sctwts[nrsct];
 149                                 }
 150
 151                                 if ((nr + nshift) < GeoPar.usrays)
 152                                 {
 153                                         vsum += sctwts[nrsct];
 154                                         sum += sctcnt[nr + nshift] * sctwts[nrsct];
 155                                 }
 156
 157                         }
 158                         projection[nr] = sum / vsum;
 159                 }
 160         }
 161
 162 }
 163
 164 void Noise_class::setQuantumNoise(INTEGER type, REAL mean, REAL calibration)
 165 {
 166         quanin = type;
 167         quanmn = mean;
 168         quancm = calibration;
 169         setQuantumNoise();
 170 }
 171
 172 void Noise_class::setQuantumNoise()
 173 {
 174         // GET AVERAGE BACKGROUND ABSORPTION
 175         back = 0.0;
 176         if (quanin == PET)
 177         {
 178                 for (INTEGER l = 0; l < Spctrm.nergy; l++)
 179                 {
 180                         back += Spctrm.engwt[l] * Spctrm.backgr[l];
 181                 }
 182         }
 183         else
 184         {
 185                 for (INTEGER l = 0; l < Spctrm.nergy; l++)
 186                 {
 187                         back += (REAL) (Spctrm.engwt[l] * exp(-Spctrm.backgr[l]));
 188                 }
 189         }
 190         c = back;
 191         // COMPUTE THE MEAN NUMBER OF PHOTONS EXITING FOR ACTUAL AND
 192         // CALIBRATION MEASUREMENTS
 193         aback = quanmn * back;
 194         if (quanin == PET)
 195                 aback = back;
 196         cback = quancm * aback;
 197
 198         // IF CALIBRATION TYPE 2 GET THE USRAYS CALIBRATION VALUES AND STORE
 199         // IN WORK AREA
 200         if (quanin == CONSTANT_RAY)
 201         {
 202
 203                 calray = new REAL[GeoPar.nrays];
 204
 205                 for (INTEGER nr = 0; nr < GeoPar.nrays; nr++)
 206                 {
 207                         calray[nr] = Gauss(cback, (REAL) sqrt(cback))
 208                                         / Gauss(cback, (REAL) sqrt(cback));
 209 #ifdef FFCOMPARE
 210
 211                         fprintf(output,"\n calray= %36.30f", calray[nr]);
 212 #endif
 213                 }
 214         }
 215         oldnp = -1;
 216
 217 }
 218
 219 REAL Noise_class::applyQuantumNoise(INTEGER np, INTEGER nr, REAL sum)
 220 {
 221         if (quanin == NONOISE)
 222                 return sum;
 223         if (quanin == PET)
 224         {
 225                 if (quanmn >= 0)
 226                         sum = (REAL) Poisson(sum);
 227         }
 228         else
 229         {
 230                 //    if (quanin > NONOISE) sum = (REAL) Poisson(sum);
 231                 if (sum < (REAL) 100.0)
 232                         sum = (REAL) Poisson(sum);
 233                 if (sum >= (REAL) 100.0)
 234                         sum = Gauss(sum, (REAL) sqrt(sum));
 235         }
 236
 237         // CHECK IF THE OBJECT IS OPAQUE FOR THIS RAY; TOO MUCH ABSORPTION
 238
 239         if ((sum <= Consts.zero * aback) && (quanin != PET))
 240         {
 241                 fprintf(output,
 242                                 "\n **** in projection %5i ray %5i no photons get through", np,
 243                                 nr);
 244                 fprintf(output,
 245                                 "\n **** in order to do tomography, more photons are needed");
 246                 fprintf(output, "\n **** program aborted\n");
 247
 248                 exit(-1);
 249         }
 250         return sum;
 251 }
 252
 253 REAL Noise_class::raysum(REAL raysum)
 254 {
 255         return quanin == PET ? raysum : exp(-raysum);
 256 }
 257
 258 REAL Noise_class::computeAttenuation(INTEGER np, INTEGER nr, REAL ind)
 259 {
 260         if (np != oldnp && quanin == CONSTANT_PROJECTION)
 261         {
 262                 c = Gauss(cback, (REAL) sqrt(cback)) / Gauss(cback, (REAL) sqrt(cback));
 263                 oldnp = np;
 264         }
 265         if (quanin == PET)
 266         {
 267                 if (quanmn < 0)
 268                 {
 269                         ind += aback;
 270                 }
 271                 else
 272                 {
 273                         ind += Poisson(aback);
 274                 }
 275         }
 276         else
 277         {
 278                 REAL a = ind;
 279                 if (quanin > NONOISE)
 280                         a /= Gauss(aback, (REAL) sqrt(aback));
 281                 if (quanin == CONSTANT_RAY)
 282                         c = calray[nr];
 283                 if (quanin == VARIABLE)
 284                         c = Gauss(cback, (REAL) sqrt(cback)) / Gauss(cback, (REAL) sqrt(cback));
 285                 ind = (REAL) log(c / a);
 286         }
 287         return ind;
 288 }