src/snark/creapr.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/creapr.cpp $
   5 $LastChangedRevision: 85 $
   6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7 $Author: agulati $
   8 **********************************************************************
   9
  10  * THIS SUBROUTINE GENERATES PROJECTION DATA
  11  * Modified to simulate emission tomography if quanin = 4 by D. Odhner
  12  */
  13
  14 #include <cstdlib>
  15 #include <cstdio>
  16 #include <cmath>
  17 #include <time.h>
  18
  19 #include "blkdta.h"
  20 #include "creacm.h"
  21 #include "geom.h"
  22 #include "spctrm.h"
  23 #include "noise.h"
  24 #include "consts.h"
  25 #include "uiod.h"
  26
  27 #include <DIGGauss.h>
  28 #include <DIGPoisson.h>
  29 #include "posit.h"
  30 #include "anglst.h"
  31 #include "raylen.h"
  32 #include "file11.h"
  33 #include "anglst.h"
  34
  35 #include "creapr.h"
  36
  37 //jk 11/18/2007 introducing variability
  38 #include "ray.h"
  39 #include "wray.h"
  40
  41
  42 // bug 221 - moved noise stuff to noise.cpp - swr - 3/2/07
  43
  44 void creapr(  ) {
  45
  46     REAL sum;
  47     REAL variability_subray_sum;
  48
  49     REAL theta, sinth, costh, raysum;
  50     REAL weight;
  51     REAL strpwt;
  52
  53     INTEGER i, k, n, nnr, np, l, nr, nb;
  54
  55     //jk 11/18/2007 introducing variability
  56     //variables needed for computation of projections through
  57     INTEGER * list;
  58     list = new INTEGER[2 * GeoPar.nelem];
  59     REAL * weight_pix;
  60     weight_pix = new REAL[2 * GeoPar.nelem];
  61     INTEGER * numb;
  62     numb = new INTEGER;
  63     REAL * snorm;
  64     snorm = new REAL;
  65
  66
  67     REAL* ind;
  68     REAL* nlo;
  69
  70     REAL emr[13][7], ax, ay, mx, my;
  71     REAL* pbase = NULL; //wei 3/2005
  72
  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
  76
  77
  78     pbase = new REAL[GeoPar.usrays];
  79     nlo = pbase;
  80
  81     //jk 09/28/2008 start the seed for noise generation here, so it is not
  82     //influenced by variability number generation
  83
  84     //if any type of noise is used
  85     if (NoisePar.quanin != NONOISE || NoisePar.addnfl
  86             || NoisePar.ultnfl || NoisePar.sctnfl) {
  87         //then set the seed
  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
  92     }
  93
  94
  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];
 100
 101
 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
 106
 107     weight = GeoPar.pinc;
 108     strpwt = weight;
 109     GeoPar.pinc /= GeoPar.nave2;
 110     GeoPar.nrays = GeoPar.usrays * GeoPar.nave2;
 111     GeoPar.midray = GeoPar.nrays / 2;
 112
 113     // STORE SINE AND COSINE OF ANGLE SUBTENDED BY EACH SUBRAY AT SOURCE
 114     // FOR DIVERGENT CASE ONLY FOR USE BY 'POSIT'
 115
 116     if (GeoPar.div) Anglst.genphi();
 117
 118
 119     // SELECT PROJECTION BY PROJECTION AND CALCULATE PROJECTION DATA
 120
 121     // IF PHANTOM IS BLANK ATTENUATION IS THAT DUE TO BACKGROUND
 122
 123     for (n = 0; n < GeoPar.nave2; n++)
 124         {
 125         for (l = 0; l < Spctrm.nergy; l++)
 126                 {
 127             emr[n][l] = NoisePar.raysum(Spctrm.backgr[l]);
 128         }
 129     }
 130
 131     for (np = 0; np < GeoPar.prjnum; np++) { //230
 132         // IF CALIBRATION TYPE 1   GET THE CALIBRATION MEASUREMENT HERE
 133
 134         Anglst.getang(np, &theta, &sinth, &costh);
 135
 136         if (GeoPar.vri) strpwt = (REAL) (MAX0(fabs(sinth), fabs(costh)) * weight);
 137
 138         nnr = 0;
 139
 140         for (nr = 0; nr < GeoPar.usrays; nr++)
 141                 {
 142
 143             if (!Creacm.modif)
 144                         {
 145                 if (Creacm.nobj > 0)
 146                                 {
 147
 148                     // CALCULATE CONTRIBUTION TO ATTENUATION BY EACH OBJECT FOR
 149                     // EACH SUBRAY (1 TO NRAYS)
 150
 151                     for (n = 0; n < GeoPar.nave2; n++)
 152                                         {
 153
 154                         posit(np, nnr, &ax, &ay, &mx, &my);
 155                         nnr++;
 156
 157                         for (k = 0; k < Creacm.nobj; k++)
 158                                                 {
 159 #ifdef FFCOMPARE
 160
 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);
 162
 163                             // rounding !!!
 164                             raylen_temp = (REAL) ((INTEGER) ((REAL) 10.0 * raylen_temp)) / (REAL) 10.0;
 165
 166                             objects[k].rayl = raylen_temp;
 167
 168 #else
 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);
 170 #endif
 171                         }
 172
 173
 174                         //////////////////////////////////////////////////////////////////////////////////////////
 175                         //if variability is set then we need to add line integrals through variability data
 176                         //to calculations of projection data
 177                         if (variab_set)
 178                                                 {
 179
 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):
 183                                                         //
 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);
 193
 194                         }
 195
 196                         //////////////////////////////////////////////////////////////////////////////////////////
 197                         //add together all the contributions to the subray
 198                         for (l = 0; l < Spctrm.nergy; l++)
 199                                                 {
 200                             //background absorption
 201                             raysum = Spctrm.backgr[l];
 202                             //variability (if present)
 203                             if (variab_set)
 204                                                         {
 205                                 variability_subray_sum = 0.0;
 206                                 //if the ray went through any pixels compute the contribution to the integral
 207                                 if (*numb != 0)
 208                                                                 {
 209                                     //for each pixels through which the ray goes
 210                                     for (nb = 0; nb < *numb; nb++)
 211                                                                         {
 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];
 214                                     }
 215                                 }
 216
 217                                 raysum += variability_subray_sum; //jk 5/21/2008  add the line ray sum for variability
 218                             }
 219                             //contributions from all the objects
 220                             for (k = 0; k < Creacm.nobj; k++)
 221                                                         {
 222                                 raysum += objects[k].den1[l] * objects[k].rayl;
 223                             }
 224
 225                             emr[n][l] = NoisePar.raysum(raysum);
 226                         } // for l
 227                     }  //for n
 228                 } //if (Creacm.nobj > 0)
 229             } //if (!Creacm.modif)
 230
 231
 232             // CALCULATE NO. OF PHOTONS OR ABSORPTION WEIGHTED BY APERTURE WT
 233             // AND THE SPECTRAL DISTRIBUTION
 234
 235             sum = 0.0;
 236             for (n = 0; n < GeoPar.nave2; n++)
 237                         {
 238                 for (l = 0; l < Spctrm.nergy; l++)
 239                                 {
 240                     sum += Creacm.aperwt[n] * Spctrm.engwt[l] * emr[n][l];
 241                 }
 242             }
 243
 244
 245             // APPLY QUANTUM NOISE
 246             pbase[nr] = NoisePar.applyQuantumNoise(np, nr, sum);
 247
 248         }
 249
 250         // COMPUTE ATTENUATION
 251         NoisePar.scatterNoise(pbase);
 252
 253
 254         // GET ACTUAL ATTENUATION AFTER CALIBRATION AND APPLY ADDITIVE, MULT
 255         // NOISE
 256         ind = pbase;
 257         for (nr = 0; nr < GeoPar.usrays; nr++)
 258                 {
 259             *ind = NoisePar.computeAttenuation(np, nr, *ind);
 260             if (GeoPar.strip) *ind *= strpwt;
 261
 262             *ind = NoisePar.multiplicativeNoise(*ind);
 263             *ind = NoisePar.additiveNoise(*ind);
 264
 265
 266             ind++;
 267         }
 268
 269         File11.WriteProj(theta, pbase, GeoPar.usrays);
 270
 271     }
 272
 273     File11.WriteEndOfProj();
 274
 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;
 282
 283     for (i= 6; i >= 0; i--)
 284         {
 285         if (variab_array[i] != NULL)  delete [] variab_array[i];  //it was allocated in creaph.cpp
 286     }
 287
 288     return;
 289 }