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) $
8 ***********************************************************
15 #include <DIGPoisson.h>
22 // bug221 - moved noise code to this class - swr - 03/03/07
24 Noise_class NoisePar = Noise_class();
26 Noise_class::Noise_class()
46 void Noise_class::resetNoise()
48 // jklukowska, bug 263, in case the noise was set
49 // previously, we need to reset it to initial values
68 Noise_class::~Noise_class()
78 void Noise_class::setAdditiveNoise(REAL mean, REAL std)
85 REAL Noise_class::additiveNoise(REAL value)
87 return addnfl ? value + Gauss(addnmn, addnsd) : value;
90 void Noise_class::setMultiplicativeNoise(REAL mean, REAL std)
97 REAL Noise_class::multiplicativeNoise(REAL value)
99 return ultnfl ? value * Gauss(ultnmn, ultnsd) : value;
102 void Noise_class::setScatter(REAL peak, REAL 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
111 sctcnt = new REAL[GeoPar.usrays];
112 nsctr = (int) (width / GeoPar.pinc + 1.0);
114 sctwts = new REAL[nsctr];
115 for (INTEGER n = 0; n < nsctr; n++)
117 sctwts[n] = (REAL) (peak * (1.0 - n * GeoPar.pinc / width));
122 void Noise_class::scatterNoise(REAL* projection)
127 // CONVOLUTE WITH THE SCATTER FUNCTION
129 for (INTEGER nr = 0; nr < GeoPar.usrays; nr++)
131 sctcnt[nr] = projection[nr];
134 for (INTEGER nr = 0; nr < GeoPar.usrays; nr++)
136 REAL sum = sctcnt[nr] * sctwts[0];
137 REAL vsum = sctwts[0];
139 for (INTEGER nrsct = 1; nrsct < nsctr; nrsct++)
141 INTEGER nshift = nrsct;
143 // CONVOLVE CAREFULLY
145 if ((nr - nshift) >= 0)
147 vsum += sctwts[nrsct];
148 sum += sctcnt[nr - nshift] * sctwts[nrsct];
151 if ((nr + nshift) < GeoPar.usrays)
153 vsum += sctwts[nrsct];
154 sum += sctcnt[nr + nshift] * sctwts[nrsct];
158 projection[nr] = sum / vsum;
164 void Noise_class::setQuantumNoise(INTEGER type, REAL mean, REAL calibration)
168 quancm = calibration;
172 void Noise_class::setQuantumNoise()
174 // GET AVERAGE BACKGROUND ABSORPTION
178 for (INTEGER l = 0; l < Spctrm.nergy; l++)
180 back += Spctrm.engwt[l] * Spctrm.backgr[l];
185 for (INTEGER l = 0; l < Spctrm.nergy; l++)
187 back += (REAL) (Spctrm.engwt[l] * exp(-Spctrm.backgr[l]));
191 // COMPUTE THE MEAN NUMBER OF PHOTONS EXITING FOR ACTUAL AND
192 // CALIBRATION MEASUREMENTS
193 aback = quanmn * back;
196 cback = quancm * aback;
198 // IF CALIBRATION TYPE 2 GET THE USRAYS CALIBRATION VALUES AND STORE
200 if (quanin == CONSTANT_RAY)
203 calray = new REAL[GeoPar.nrays];
205 for (INTEGER nr = 0; nr < GeoPar.nrays; nr++)
207 calray[nr] = Gauss(cback, (REAL) sqrt(cback))
208 / Gauss(cback, (REAL) sqrt(cback));
211 fprintf(output,"\n calray= %36.30f", calray[nr]);
219 REAL Noise_class::applyQuantumNoise(INTEGER np, INTEGER nr, REAL sum)
221 if (quanin == NONOISE)
226 sum = (REAL) Poisson(sum);
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));
237 // CHECK IF THE OBJECT IS OPAQUE FOR THIS RAY; TOO MUCH ABSORPTION
239 if ((sum <= Consts.zero * aback) && (quanin != PET))
242 "\n **** in projection %5i ray %5i no photons get through", np,
245 "\n **** in order to do tomography, more photons are needed");
246 fprintf(output, "\n **** program aborted\n");
253 REAL Noise_class::raysum(REAL raysum)
255 return quanin == PET ? raysum : exp(-raysum);
258 REAL Noise_class::computeAttenuation(INTEGER np, INTEGER nr, REAL ind)
260 if (np != oldnp && quanin == CONSTANT_PROJECTION)
262 c = Gauss(cback, (REAL) sqrt(cback)) / Gauss(cback, (REAL) sqrt(cback));
273 ind += Poisson(aback);
279 if (quanin > NONOISE)
280 a /= Gauss(aback, (REAL) sqrt(aback));
281 if (quanin == CONSTANT_RAY)
283 if (quanin == VARIABLE)
284 c = Gauss(cback, (REAL) sqrt(cback)) / Gauss(cback, (REAL) sqrt(cback));
285 ind = (REAL) log(c / a);