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/creaph.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 **********************************************************************
10 GENERATES A TEST PHANTOM
12 * jk 11/18/2007 introducing variability
13 * jk 05/30/2008 redesigned and corrected implementation of variability
14 * the variability has to be added to each pixel at each energy
15 * level (polychromatic data) - in case of monochromatic data
16 * there is simply only one energy level
17 * each pixel value is calculated according to the formula:
18 * 1/(NAVE1^2) sum_{k=1 to NAVE1^2} sum_{j=1 to NOBJ} delta_{k,j}
19 * sum_{i=1 to NERGY} (1+g_{j,i}) d_{j,i}p_{i}
20 * (see snark14 manual for explanation of the above value)
41 //jk 11/18/2007 introducing variability
42 //following header files needed for adding tissue variability
45 //jk 05/30/08 a new header file needed
51 int i, j, k, n, ii, jj, index, ne; //loop counter variables
52 REAL g; //temp variable for Gaussian value
53 INTEGER rows; //temp variable for storing number of items in the rows
67 REAL denobj[7]; //temporarily stores the density of a pixel during computations
68 //the value for each energy level is needed (max 7)
76 REAL densit[7]; //temporarily stores the density of a pixel during computations
77 //the value for each energy level is needed (max 7)
88 // ALLOCATE WORKING SPACE - ONLY ONE ROW OF PIXELS IS NEEDED AT ONCE
90 REAL * test_aray[8]; //stores a row of pixels at a time (alocated below)
91 //each pixel needs a value for each energy level (7 of them)
92 //the last item is used to store the average value
93 REAL* coord_aray = NULL;
95 //allocate enough space for one row (up to 7 energy levels + one more for average
96 for (ne = 0; ne < 8; ne++)
97 test_aray[ne] = new REAL[GeoPar.nelem];
99 //jk 11/18/2007 introducing variability
100 //allocate memory for the array that holds variability data
101 //(for each of the energy levels)
102 //do so only if the variability is not equal to zero
103 memset(variab_array, 0, 7);
106 for (ne = 0; ne < Spctrm.nergy; ne++)
108 variab_array[ne] = new REAL[GeoPar.nelem * GeoPar.nelem];
109 memset(variab_array[ne], 0, 7);
114 //jk 09/19/2008 modifying how polychromatic pseudo projections are computed
115 //in the iterative runs of beam hardening correction
116 memset(phantom_poly_array, 0, 7);
117 for (ne = 0; ne < Spctrm.nergy; ne++)
119 phantom_poly_array[ne] = new REAL[GeoPar.nelem * GeoPar.nelem];
120 memset(phantom_poly_array[ne], 0, GeoPar.nelem * GeoPar.nelem);
123 // THIS RUN IS FIRST GENERATION OF PHANTOM
124 // ALLOCATE SPACE FOR SUBPIXELS
125 npts = GeoPar.nelem * Creacm.nave1;
127 coord_aray = new REAL[npts];
129 nmid = (REAL) (npts / 2);
131 // PRESET THE COORDINATES OF EACH SUBPIXEL CENTER
132 for (n = 0; n < npts; n++)
134 coord_aray[n] = (((REAL) (n)) - nmid) * GeoPar.pixsiz
135 / ((REAL) (Creacm.nave1));
138 // COMPUTE THE NUMBER OF SAMPLES PER PIXEL FOR AVERAGING
139 sample = Creacm.nave1 * Creacm.nave1;
141 //set the seed for variability creation
144 if (Creacm.varseed < 0)
146 else if (Creacm.varseed > 0)
147 Srand(Creacm.varseed);
149 Srand(1); //this generates the default sequence as if Srand was
150 //not called (it is called here to start this default
154 for (ne = 0; ne < Spctrm.nergy; ne++)
160 // GET EACH PIXEL DENSITY
161 for (j = 0; j < GeoPar.nelem; j++) // for each pix row
166 for (k = 0; k < GeoPar.nelem; k++)
169 test_aray[ne][k] = 0;
172 // IF NO OBJECT SPECIFIED THEN GENERATE BLANK PHANTOM
177 for (index = 0; index < Creacm.nobj; index++)
180 itype = objects[index].type;
182 // SKIP IF OBJECT WAS DELETED FROM LIST
183 co = objects[index].cosang;
184 si = objects[index].sinang;
185 if (fabs(co - 1) < 0.000001)
187 if (fabs(co + 1) < 0.000001)
189 if (fabs(co) < 0.000001)
191 if (fabs(si - 1) < 0.000001)
193 if (fabs(si + 1) < 0.000001)
195 if (fabs(si) < 0.000001)
197 unew = objects[index].u;
198 vnew = objects[index].v;
203 // 9/12/2012 jklukowska
204 //if the unew and/or vnew are smaller than the size of the pixel,
205 //use the size of the pixel as max
206 diam = (REAL) 1.5 * MAX0(MAX0(unew, vnew), GeoPar.pixsiz);
207 // diam = (REAL) 1.5 * MAX0(unew, vnew);
208 vbyusq = (vnew / unew) * (vnew / unew);
209 xcntr = objects[index].cx;
210 ycntr = objects[index].cy;
212 denobj[ne] = objects[index].den1[ne];
217 for (jj = 0; jj < Creacm.nave1; jj++)
220 // SHIFT THE ORIGIN TO THE CENTER OF THE OBJECT
221 yo = -coord_aray[ny] - ycntr;
224 // TEST IF THE POINT (X,Y) IS INSIDE THE INDEXED OBJECT
225 // IF THE ROW WE ARE AT IS OUTSIDE MAXIMUM EXTENT OF
226 // OBJECT, GET NEXT OBJECT
231 // MOVE ACROSS COLUMNS
234 // for each pixel in row
235 for (i = 0; i < GeoPar.nelem; i++)
240 for (ii = 0; ii < Creacm.nave1; ii++)
243 xo = coord_aray[nx] - xcntr;
246 // IF COLUMN IS LEFT OF EXTENT, INCREASE X
251 // IF COLUMN IS RIGHT OF EXTENT,UPDATE PIXEL ARRAY,THEN GET NEXT Y
255 test_aray[ne][i] += densit[ne];
260 // pixel within object radius
261 // ROTATE BY ANGLE OF TILT
262 xc = (REAL) fabs(xo * covru + yo * sivru);
265 if (fabs(xc) < 0.000001)
267 if (fabs(xc + 1.0) < 0.000001)
269 if (fabs(xc - 1.0) < 0.000001)
275 yc = (yo * covrv - xo * sivrv);
278 if (fabs(yc) < 0.000001)
280 if (fabs(yc - 1.0) < 0.000001)
282 if (fabs(yc + 1.0) < 0.000001)
289 case SOT_elip: // ELLIPSE
292 && ((xc * xc + yc * yc - 1.0)
295 // added tolerance - swr 11/26/04
296 densit[ne] += denobj[ne];
300 case SOT_rect: // RECTANGLE
303 densit[ne] += denobj[ne];
307 case SOT_tria: // TRIANGLE
308 // CHECK TO SEE IF POINT IS CONTAINED IN TRIANGULAR AREA
310 && (xc <= (1.0 - yc)))
312 densit[ne] += denobj[ne];
313 // NOW THE POINT IS INSIDE
314 // THE OBJECTS BOUNDARIES
318 case SOT_segm: // SEGMENT
322 // CHECK TO SEE IF POINT IS CONTAINED IN TRIANGULAR AREA
323 if (xc <= (1.0 - yc))
325 //if(xc > 1.0-yc) goto L60;
327 // TEST IF POINT IS BEYOND CIRCULAR ARC
335 // NOW THE POINT IS INSIDE THE OBJECTS BOUNDARIES
344 case SOT_sect: // SECTOR
346 // CHECK TO SEE IF POINT IS CONTAINED IN TRIANGULAR AREA
347 if (xc <= (1.0 - yc))
349 //if(xc > 1.0-yc) goto L60;
350 // TEST IF POINT IS BEYOND CIRCULAR ARC
357 // NOW THE POINT IS INSIDE THE OBJECTS BOUNDARIES
358 densit[ne] += denobj[ne];
368 test_aray[ne][i] += densit[ne];
379 rows = j * GeoPar.nelem; //this avoids repetitive multiplications below
381 for (k = 0; k < GeoPar.nelem; k++)
383 // TAKE THE AVERAGE BY DIVIDING BY THE NUMBER OF SAMPLES
384 if (Creacm.nave1 != 1)
386 test_aray[ne][k] /= sample;
388 //jk 11/18/2007 introducing variability
389 //ADD TISSUE VARIABILITY TO THE PHANTOM DATA
390 //variab_array[row][column] is assigned the value of variability for
393 //if the pixel belongs to the background (i.e. the density value
394 //is zero for each energy level) the value is set to 0
395 //(no variability is added to the background pixels)
399 if (fabs(test_aray[ne][k]) > Consts.zero)
401 g = Gauss(0.0, variability) * test_aray[ne][k];
402 test_aray[ne][k] += g;
403 variab_array[ne][k + rows] = g;
409 variab_array[ne][k + rows] = 0;
414 phantom_poly_array[ne][k + rows] = test_aray[ne][k];
416 for (k = 0; k < GeoPar.nelem; k++)
418 assert(phantom_poly_array[ne][rows + k] == test_aray[ne][k]);
421 //write row of pixels to file11 (use the weighted average
422 //of all the energy levels)
423 //File11.WriteRowOfPhanPixels(test_aray[7], GeoPar.neleM);
424 // jk 07/08/2008 use the value of the first energy level
425 // for phantom, instead of the weighted average)
427 File11.WriteRowOfPhanPixels(test_aray[0], GeoPar.nelem);
431 File11.WriteEndOfPhan();
433 //delete allocated memory
434 if (coord_aray != NULL)
437 for (ne = 0; ne < 8; ne++)
439 if (test_aray[ne] != NULL)
440 delete[] test_aray[ne];
442 //delete[] test_aray;