src/snark/creaph.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/creaph.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  **********************************************************************
   9
  10  GENERATES A TEST PHANTOM
  11  *
  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)
  21  */
  22
  23 #include <cstdlib>
  24 #include <cstdio>
  25 #include <cmath>
  26 #include <string.h>
  27 #include <time.h>
  28
  29 #include "blkdta.h"
  30 #include "creacm.h"
  31 #include "geom.h"
  32 #include "objects.h"
  33 #include "uiod.h"
  34
  35 #include "file11.h"
  36
  37 #include "creaph.h"
  38
  39 #include "assert.h"
  40
  41 //jk 11/18/2007 introducing variability
  42 //following header files needed for adding tissue variability
  43 #include "DIGGauss.h"
  44 #include "consts.h"
  45 //jk  05/30/08   a new header file needed
  46 #include "spctrm.h"
  47
  48 void creaph()
  49 {
  50
  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
  54
  55         REAL nmid;
  56
  57         REAL unew;
  58         REAL vnew;
  59         REAL covru;
  60         REAL covrv;
  61         REAL sivru;
  62         REAL sivrv;
  63         REAL diam;
  64         REAL vbyusq;
  65         REAL xcntr;
  66         REAL ycntr;
  67         REAL denobj[7]; //temporarily stores the density of a pixel during computations
  68         //the value for each energy level is needed (max 7)
  69
  70         REAL co;
  71         REAL si;
  72
  73         INTEGER ny;
  74         REAL yo;
  75         INTEGER nx;
  76         REAL densit[7]; //temporarily stores the density of a pixel during computations
  77         //the value for each energy level is needed (max 7)
  78
  79         INTEGER ny0;
  80         INTEGER nx0;
  81
  82         INTEGER npts;
  83         INTEGER sample;
  84         INTEGER itype;
  85
  86         REAL yc, xc, xo;
  87
  88         // ALLOCATE WORKING SPACE - ONLY ONE ROW OF PIXELS IS NEEDED AT ONCE
  89
  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;
  94
  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];
  98
  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);
 104         if (variab_set)
 105         {
 106                 for (ne = 0; ne < Spctrm.nergy; ne++)
 107                 {
 108                         variab_array[ne] = new REAL[GeoPar.nelem * GeoPar.nelem];
 109                         memset(variab_array[ne], 0, 7);
 110                 }
 111
 112         }
 113
 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++)
 118         {
 119                 phantom_poly_array[ne] = new REAL[GeoPar.nelem * GeoPar.nelem];
 120                 memset(phantom_poly_array[ne], 0, GeoPar.nelem * GeoPar.nelem);
 121         }
 122
 123         // THIS RUN IS FIRST GENERATION OF PHANTOM
 124         // ALLOCATE SPACE FOR SUBPIXELS
 125         npts = GeoPar.nelem * Creacm.nave1;
 126
 127         coord_aray = new REAL[npts];
 128
 129         nmid = (REAL) (npts / 2);
 130
 131         // PRESET THE COORDINATES OF EACH SUBPIXEL CENTER
 132         for (n = 0; n < npts; n++)
 133         {
 134                 coord_aray[n] = (((REAL) (n)) - nmid) * GeoPar.pixsiz
 135                                 / ((REAL) (Creacm.nave1));
 136         }
 137
 138         // COMPUTE THE NUMBER OF SAMPLES PER PIXEL FOR AVERAGING
 139         sample = Creacm.nave1 * Creacm.nave1;
 140
 141         //set the seed for variability creation
 142         if (variab_set)
 143         {
 144                 if (Creacm.varseed < 0)
 145                         Srand(time(NULL));
 146                 else if (Creacm.varseed > 0)
 147                         Srand(Creacm.varseed);
 148                 else
 149                         Srand(1); //this generates the default sequence as if Srand was
 150                 //not called (it is called here to start this default
 151                 //sequence over
 152         }
 153
 154         for (ne = 0; ne < Spctrm.nergy; ne++)
 155         {
 156
 157                 ny0 = -Creacm.nave1;
 158                 nx0 = 0;
 159
 160                 // GET EACH PIXEL DENSITY
 161                 for (j = 0; j < GeoPar.nelem; j++)              // for each pix row
 162                 {
 163
 164                         ny0 += Creacm.nave1;
 165
 166                         for (k = 0; k < GeoPar.nelem; k++)
 167                         {
 168                                 // clear row buffer
 169                                 test_aray[ne][k] = 0;
 170                         }
 171
 172                         // IF NO OBJECT SPECIFIED THEN GENERATE BLANK PHANTOM
 173                         if (Creacm.nobj > 0)
 174                         {
 175
 176                                 // for each object
 177                                 for (index = 0; index < Creacm.nobj; index++)
 178                                 {
 179
 180                                         itype = objects[index].type;
 181
 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)
 186                                                 co = 1.;
 187                                         if (fabs(co + 1) < 0.000001)
 188                                                 co = -1.;
 189                                         if (fabs(co) < 0.000001)
 190                                                 co = 0.;
 191                                         if (fabs(si - 1) < 0.000001)
 192                                                 si = 1.;
 193                                         if (fabs(si + 1) < 0.000001)
 194                                                 si = -1.;
 195                                         if (fabs(si) < 0.000001)
 196                                                 si = 0.;
 197                                         unew = objects[index].u;
 198                                         vnew = objects[index].v;
 199                                         covru = co / unew;
 200                                         covrv = co / vnew;
 201                                         sivru = si / unew;
 202                                         sivrv = si / vnew;
 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;
 211
 212                                         denobj[ne] = objects[index].den1[ne];
 213
 214                                         ny = ny0;
 215
 216                                         // for each subpix
 217                                         for (jj = 0; jj < Creacm.nave1; jj++)
 218                                         {
 219
 220                                                 // SHIFT THE ORIGIN TO THE CENTER OF THE OBJECT
 221                                                 yo = -coord_aray[ny] - ycntr;
 222                                                 ny++;
 223
 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
 227                                                 if (fabs(yo) > diam)
 228                                                         goto L95;
 229                                                 // goto nex object
 230
 231                                                 // MOVE ACROSS COLUMNS
 232
 233                                                 nx = nx0;
 234                                                 // for each pixel in row
 235                                                 for (i = 0; i < GeoPar.nelem; i++)
 236                                                 {
 237                                                         densit[ne] = 0.0;
 238
 239                                                         // for each subpixel
 240                                                         for (ii = 0; ii < Creacm.nave1; ii++)
 241                                                         {
 242
 243                                                                 xo = coord_aray[nx] - xcntr;
 244                                                                 nx++;
 245
 246                                                                 // IF COLUMN IS LEFT OF EXTENT, INCREASE X
 247                                                                 if (xo < -diam)
 248                                                                         goto L70;
 249                                                                 // goto next column
 250
 251                                                                 // IF COLUMN IS RIGHT OF EXTENT,UPDATE PIXEL ARRAY,THEN GET NEXT Y
 252                                                                 if (xo > diam)
 253                                                                 {
 254                                                                         // HERE xo .gt. diam
 255                                                                         test_aray[ne][i] += densit[ne];
 256
 257                                                                         goto L90;
 258                                                                         // goto next row
 259                                                                 }
 260                                                                 // pixel within object radius
 261                                                                 // ROTATE BY ANGLE OF TILT
 262                                                                 xc = (REAL) fabs(xo * covru + yo * sivru);
 263
 264                                                                 // CHECK XC FOR FUZZ
 265                                                                 if (fabs(xc) < 0.000001)
 266                                                                         xc = 0.0;
 267                                                                 if (fabs(xc + 1.0) < 0.000001)
 268                                                                         xc = -1.0;
 269                                                                 if (fabs(xc - 1.0) < 0.000001)
 270                                                                         xc = 1.0;
 271
 272                                                                 if (xc <= 1.0)
 273                                                                 {
 274
 275                                                                         yc = (yo * covrv - xo * sivrv);
 276
 277                                                                         // CHECK YC FOR FUZZ
 278                                                                         if (fabs(yc) < 0.000001)
 279                                                                                 yc = 0.0;
 280                                                                         if (fabs(yc - 1.0) < 0.000001)
 281                                                                                 yc = 1.0;
 282                                                                         if (fabs(yc + 1.0) < 0.000001)
 283                                                                                 yc = -1.0;
 284
 285                                                                         if (yc <= 1.0)
 286                                                                         {
 287                                                                                 switch (itype)
 288                                                                                 {
 289                                                                                 case SOT_elip: // ELLIPSE
 290
 291                                                                                         if ((yc >= -1.0)
 292                                                                                                         && ((xc * xc + yc * yc - 1.0)
 293                                                                                                                         < 0.000001))
 294                                                                                         {
 295                                                                                                 // added tolerance - swr 11/26/04
 296                                                                                                 densit[ne] += denobj[ne];
 297                                                                                         }
 298                                                                                         break;
 299
 300                                                                                 case SOT_rect: // RECTANGLE
 301                                                                                         if (yc >= -1.0)
 302                                                                                         {
 303                                                                                                 densit[ne] += denobj[ne];
 304                                                                                         }
 305                                                                                         break;
 306
 307                                                                                 case SOT_tria: // TRIANGLE
 308                                                                                         // CHECK TO SEE IF POINT IS CONTAINED IN TRIANGULAR AREA
 309                                                                                         if ((yc >= 0.0)
 310                                                                                                         && (xc <= (1.0 - yc)))
 311                                                                                         {
 312                                                                                                 densit[ne] += denobj[ne];
 313                                                                                                 // NOW THE POINT IS INSIDE
 314                                                                                                 // THE OBJECTS BOUNDARIES
 315                                                                                         }
 316                                                                                         break;
 317
 318                                                                                 case SOT_segm: // SEGMENT
 319                                                                                         if (yc <= 0.0)
 320                                                                                         {
 321
 322                                                                                                 // CHECK TO SEE IF POINT IS CONTAINED IN TRIANGULAR AREA
 323                                                                                                 if (xc <= (1.0 - yc))
 324                                                                                                 {
 325                                                                                                         //if(xc > 1.0-yc) goto L60;
 326
 327                                                                                                         // TEST IF POINT IS BEYOND CIRCULAR ARC
 328                                                                                                         if ((xc * xc)
 329                                                                                                                         <= (1.0
 330                                                                                                                                         + yc
 331                                                                                                                                                         * vbyusq
 332                                                                                                                                                         * (2.0
 333                                                                                                                                                                         - yc)))
 334                                                                                                         {
 335                                                                                                                 // NOW THE POINT IS INSIDE THE OBJECTS BOUNDARIES
 336                                                                                                                 densit[ne] +=
 337                                                                                                                                 denobj[ne];
 338
 339                                                                                                         }
 340                                                                                                 }
 341                                                                                         }
 342                                                                                         break;
 343
 344                                                                                 case SOT_sect: // SECTOR
 345
 346                                                                                         // CHECK TO SEE IF POINT IS CONTAINED IN TRIANGULAR AREA
 347                                                                                         if (xc <= (1.0 - yc))
 348                                                                                         {
 349                                                                                                 //if(xc > 1.0-yc) goto L60;
 350                                                                                                 // TEST IF POINT IS BEYOND CIRCULAR ARC
 351                                                                                                 if ((xc * xc)
 352                                                                                                                 <= (1.0
 353                                                                                                                                 + yc * vbyusq
 354                                                                                                                                                 * (2.0
 355                                                                                                                                                                 - yc)))
 356                                                                                                 {
 357                                                                                                         // NOW THE POINT IS INSIDE THE OBJECTS BOUNDARIES
 358                                                                                                         densit[ne] += denobj[ne];
 359                                                                                                 }
 360                                                                                         }
 361                                                                                         break;
 362
 363                                                                                 }
 364                                                                         }
 365                                                                 }
 366                                                                 L70: ;
 367                                                         } // for(ii
 368                                                         test_aray[ne][i] += densit[ne];
 369
 370                                                 } //for(i
 371                                                 L90: ;
 372                                         } //for(jj
 373
 374                                         L95: ;
 375                                 } //for(index
 376
 377                         }
 378
 379                         rows = j * GeoPar.nelem; //this avoids repetitive multiplications below
 380
 381                         for (k = 0; k < GeoPar.nelem; k++)
 382                         {
 383                                 // TAKE THE AVERAGE BY DIVIDING BY THE NUMBER OF SAMPLES
 384                                 if (Creacm.nave1 != 1)
 385                                 {
 386                                         test_aray[ne][k] /= sample;
 387                                 }
 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
 391                                 //each pixel
 392
 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)
 396                                 if (variab_set)
 397                                 {
 398
 399                                         if (fabs(test_aray[ne][k]) > Consts.zero)
 400                                         {
 401                                                 g = Gauss(0.0, variability) * test_aray[ne][k];
 402                                                 test_aray[ne][k] += g;
 403                                                 variab_array[ne][k + rows] = g;
 404
 405                                         }
 406
 407                                         else
 408                                         {
 409                                                 variab_array[ne][k + rows] = 0;
 410                                         }
 411
 412                                 }
 413
 414                                 phantom_poly_array[ne][k + rows] = test_aray[ne][k];
 415                         }
 416                         for (k = 0; k < GeoPar.nelem; k++)
 417                         {
 418                                 assert(phantom_poly_array[ne][rows + k] == test_aray[ne][k]);
 419                         }
 420
 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)
 426                         if (ne == 0)
 427                                 File11.WriteRowOfPhanPixels(test_aray[0], GeoPar.nelem);
 428
 429                 }  //for j
 430         } //for ne
 431         File11.WriteEndOfPhan();
 432
 433         //delete allocated memory
 434         if (coord_aray != NULL)
 435                 delete[] coord_aray;
 436
 437         for (ne = 0; ne < 8; ne++)
 438         {
 439                 if (test_aray[ne] != NULL)
 440                         delete[] test_aray[ne];
 441         }
 442         //delete[] test_aray;
 443
 444         return;
 445 }