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/creapr.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 **********************************************************************
10 * THIS SUBROUTINE GENERATES PROJECTION DATA
11 * Modified to simulate emission tomography if quanin = 4 by D. Odhner
28 #include <DIGPoisson.h>
37 //jk 11/18/2007 introducing variability
42 // bug 221 - moved noise stuff to noise.cpp - swr - 3/2/07
47 REAL variability_subray_sum;
49 REAL theta, sinth, costh, raysum;
53 INTEGER i, k, n, nnr, np, l, nr, nb;
55 //jk 11/18/2007 introducing variability
56 //variables needed for computation of projections through
58 list = new INTEGER[2 * GeoPar.nelem];
60 weight_pix = new REAL[2 * GeoPar.nelem];
70 REAL emr[13][7], ax, ay, mx, my;
71 REAL* pbase = NULL; //wei 3/2005
73 // COMPUTE PROJECTION BY PROJECTION HENCE WE NEED SPACE FOR
74 // ONLY ONE PROJECTION OF LENGTH NRAYS
75 // NLO WILL DENOTE THE FIRST RAY AND NHI THE LAST RAY
78 pbase = new REAL[GeoPar.usrays];
81 //jk 09/28/2008 start the seed for noise generation here, so it is not
82 //influenced by variability number generation
84 //if any type of noise is used
85 if (NoisePar.quanin != NONOISE || NoisePar.addnfl
86 || NoisePar.ultnfl || NoisePar.sctnfl) {
88 if (Creacm.iseed < 0) Srand(time(NULL)); // bug 184 - swr - 11/04/05
89 else if (Creacm.iseed > 0) Srand(Creacm.iseed);
90 else Srand(1); //set the default seed so that the default sequence
91 //is the same regardless of presence of variability before
95 //jk 11/18/2007 introducing variability
96 //stores the value of projection, written to the FILE11 after each projection
97 //hence we need only one
98 REAL * pseudoray = NULL;
99 pseudoray = new REAL[GeoPar.usrays];
102 // SUBDIVIDE RAYS FOR AVERAGING OVER NAVE2 SUBRAYS/SUBSTRIPS
103 // CHANGE PINC, NRAYS AND MIDRAY IN /GEOM/ FOR USE BY POSIT.
104 // WHEN RDPROJ IS DONE THEESE VALUES WILL BE RESTORED
105 // SAVE PINC IN WEIGHT FOR STRIP CASE
107 weight = GeoPar.pinc;
109 GeoPar.pinc /= GeoPar.nave2;
110 GeoPar.nrays = GeoPar.usrays * GeoPar.nave2;
111 GeoPar.midray = GeoPar.nrays / 2;
113 // STORE SINE AND COSINE OF ANGLE SUBTENDED BY EACH SUBRAY AT SOURCE
114 // FOR DIVERGENT CASE ONLY FOR USE BY 'POSIT'
116 if (GeoPar.div) Anglst.genphi();
119 // SELECT PROJECTION BY PROJECTION AND CALCULATE PROJECTION DATA
121 // IF PHANTOM IS BLANK ATTENUATION IS THAT DUE TO BACKGROUND
123 for (n = 0; n < GeoPar.nave2; n++)
125 for (l = 0; l < Spctrm.nergy; l++)
127 emr[n][l] = NoisePar.raysum(Spctrm.backgr[l]);
131 for (np = 0; np < GeoPar.prjnum; np++) { //230
132 // IF CALIBRATION TYPE 1 GET THE CALIBRATION MEASUREMENT HERE
134 Anglst.getang(np, &theta, &sinth, &costh);
136 if (GeoPar.vri) strpwt = (REAL) (MAX0(fabs(sinth), fabs(costh)) * weight);
140 for (nr = 0; nr < GeoPar.usrays; nr++)
148 // CALCULATE CONTRIBUTION TO ATTENUATION BY EACH OBJECT FOR
149 // EACH SUBRAY (1 TO NRAYS)
151 for (n = 0; n < GeoPar.nave2; n++)
154 posit(np, nnr, &ax, &ay, &mx, &my);
157 for (k = 0; k < Creacm.nobj; k++)
161 REAL raylen_temp = raylen(objects[k].type, ax, ay, mx, my, objects[k].cx, objects[k].cy, objects[k].u, objects[k].v, objects[k].cosang, objects[k].sinang);
164 raylen_temp = (REAL) ((INTEGER) ((REAL) 10.0 * raylen_temp)) / (REAL) 10.0;
166 objects[k].rayl = raylen_temp;
169 objects[k].rayl = raylen(objects[k].type, ax, ay, mx, my, objects[k].cx, objects[k].cy, objects[k].u, objects[k].v, objects[k].cosang, objects[k].sinang);
174 //////////////////////////////////////////////////////////////////////////////////////////
175 //if variability is set then we need to add line integrals through variability data
176 //to calculations of projection data
180 // get information about the pixels through which the ray passes
181 // (it will be used below when the contributions are computed
182 // for each energy level):
184 // numb THE NUMBER OF PIXELS A RAY PASSES THROUGH
185 // snorm SUM OF SQUARES OF WEIGHTS
186 // weight_pix THE LENGTH OF A RAY IN A PIXEL(CORRESPONDING TO
187 // ENTRY IN THE ARRAY LIST)
188 // list THE LIST OF INDECIES INTO THE ARRAY variab_array THAT CORRESP
189 // TO THE SILECTED PIXELS
190 // values of ax, ay, mx and my are reused from a previous call to
191 // posit(np, nnr, &ax, &ay, &mx, &my);
192 wray(np, nnr, list, weight_pix, numb, snorm, &ax, &ay, &mx, &my);
196 //////////////////////////////////////////////////////////////////////////////////////////
197 //add together all the contributions to the subray
198 for (l = 0; l < Spctrm.nergy; l++)
200 //background absorption
201 raysum = Spctrm.backgr[l];
202 //variability (if present)
205 variability_subray_sum = 0.0;
206 //if the ray went through any pixels compute the contribution to the integral
209 //for each pixels through which the ray goes
210 for (nb = 0; nb < *numb; nb++)
212 //add the contribution of the pixel to subray at energy level l
213 variability_subray_sum += variab_array[l][list[nb]] * weight_pix[nb];
217 raysum += variability_subray_sum; //jk 5/21/2008 add the line ray sum for variability
219 //contributions from all the objects
220 for (k = 0; k < Creacm.nobj; k++)
222 raysum += objects[k].den1[l] * objects[k].rayl;
225 emr[n][l] = NoisePar.raysum(raysum);
228 } //if (Creacm.nobj > 0)
229 } //if (!Creacm.modif)
232 // CALCULATE NO. OF PHOTONS OR ABSORPTION WEIGHTED BY APERTURE WT
233 // AND THE SPECTRAL DISTRIBUTION
236 for (n = 0; n < GeoPar.nave2; n++)
238 for (l = 0; l < Spctrm.nergy; l++)
240 sum += Creacm.aperwt[n] * Spctrm.engwt[l] * emr[n][l];
245 // APPLY QUANTUM NOISE
246 pbase[nr] = NoisePar.applyQuantumNoise(np, nr, sum);
250 // COMPUTE ATTENUATION
251 NoisePar.scatterNoise(pbase);
254 // GET ACTUAL ATTENUATION AFTER CALIBRATION AND APPLY ADDITIVE, MULT
257 for (nr = 0; nr < GeoPar.usrays; nr++)
259 *ind = NoisePar.computeAttenuation(np, nr, *ind);
260 if (GeoPar.strip) *ind *= strpwt;
262 *ind = NoisePar.multiplicativeNoise(*ind);
263 *ind = NoisePar.additiveNoise(*ind);
269 File11.WriteProj(theta, pbase, GeoPar.usrays);
273 File11.WriteEndOfProj();
275 //Clean up all the allocated memory
276 if (pbase != NULL) delete[] pbase;
277 if (list != NULL) delete [] list;
278 if (weight_pix != NULL) delete [] weight_pix;
279 if (numb != NULL) delete numb;
280 if (snorm != NULL) delete snorm;
281 if (pseudoray != NULL) delete pseudoray;
283 for (i= 6; i >= 0; i--)
285 if (variab_array[i] != NULL) delete [] variab_array[i]; //it was allocated in creaph.cpp