src/snark/create_phantom.c

   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/create_phantom.c $
   5 $LastChangedRevision: 85 $
   6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7 $Author: agulati $
   8 **********************************************************************
   9
  10  Previously part of SuperSNARK
  11  */
  12
  13 /*--------------------------------------------------------------------------
  14  * Super Snark V 1.0
  15  * Written by Jingsheng Zheng... November 1993
  16  * Modified by Jolyon Browne ... June 1993
  17  * --------------------------------------------------------------------------*/
  18
  19 #include <string.h>
  20 #include <strings.h>
  21 #include <stdlib.h>
  22
  23 #include "experimenter.h"
  24 #include "blanks.h"
  25 #include "errorc.h"
  26 #include "create_phantom.h"
  27
  28 /*--------------------------------------------------------------------------
  29  * This fuction generates the input lines that defines the elemental
  30  * objects as required by the CREATE command.
  31  *
  32  * INPUTS:
  33  * phan_flname - phantom filename (Note: the list of phantom filenames are
  34  * stored in the ensemble data file).
  35  * seed - seed used for the random number generator.
  36  * navel, nelem,pixsz - modifiers described by Set 4 of the SNARK CREATE
  37  * input sequence (see the SNARK manual for details).
  38  *
  39  * OUTPUTS:
  40  * structure - string array of length num_str containing the elemental object
  41  * description.
  42  * num_str - number of structures in the phantom.
  43  * num_pairs - for PAIRed structures, this is the number of pairs of
  44  * structures.
  45  * pairs - array of length 2*num_pairs containing the indicies of all PAIRed
  46  * structures.
  47  * ----------------------------------------------------------------------------*/
  48
  49 void create_phantom(
  50         char* phan_flname,
  51         long seed,
  52         int navel,
  53         int nelem,
  54         float pixsz,
  55         char*** structure, //jk 01/01/08 modified to remove
  56         // restriction for number of structures
  57         int* num_str,
  58         int* num_pairs,
  59         int** pairs, //jk 01/01/08 modified to remove
  60         // restriction for number of structures
  61         int iroi_used //jk  09/07/2008   added a check to verify that
  62         //density of one out of each paired  structures is zero
  63         )
  64 {
  65     float cx, cy, u, v, ang, den1, den2, prob, tmp;
  66     float scale, sd;
  67     int i, nstr = 0, npair = 0, tot_nstr = 0; //jk 01/01/2008 added tot_nstr
  68     //variable to keep track of the number
  69     //of structures in the phantom
  70     int scaleCount = 0; //keeps track of how many times line starting with SCALE
  71                         //appeared in the input file (it should be exactly 1
  72                         //at the end of the file construction).
  73
  74     char shape[10], keyword[4], string[MAXLINESIZE];
  75     FILE *phanfl;
  76
  77     *num_pairs = 0;
  78
  79     //jk 09/20/2008
  80     char spect[10], spectType[15];
  81     int erg, nerg;
  82     int poly = 0;
  83     int change = 0;
  84     char toc[14][MAXLINESIZE];
  85     int polyErg[7];
  86     int polyProb[7];
  87     char string_alt[MAXLINESIZE];
  88     long phantom_seed = (long) (100000000 * drand48() );
  89
  90     if ((phanfl = fopen(phan_flname, "r")) == NULL) {
  91         errorc("In create_phantom.c: error in opening file ", phan_flname);
  92     }
  93
  94     /* print out beginning of CREATE input sequence */
  95
  96     fprintf(pstream, "CREATE\n");
  97     fprintf(pstream, "%s, seed used in experimenter to generate this file = %ld\n", phan_flname, seed);
  98
  99     // jk 09/20/2008 adding polychromatic data option
 100     if (fgets(string, sizeof (string), phanfl) == NULL) {
 101         errorc("error: premature end of file encountered ", phan_flname);
 102     }
 103     while ( isSkip(string) ) //it's a comment line,
 104     {                           //echo it to output and read next line
 105         if( isComment(string) )
 106         {
 107             fprintf(pstream, "%s", string);
 108             if ( fgets(string, sizeof (string), phanfl) == NULL) break;
 109         }
 110     }
 111     if (strncasecmp(string, "SPEC", 4) == 0)
 112         {
 113         sscanf(string, "%s%s%i", spect, spectType, &erg);
 114
 115                 if (strncasecmp(spectType, "mono", 4) == 0)
 116             fprintf(pstream, "SPECTRUM MONOCHROMATIC %i\n", erg);
 117
 118                 else if (strncasecmp(spectType, "poly", 4) == 0)
 119                 {
 120             nerg = erg;
 121             if (nerg > 7)
 122                         {
 123                 errorc("error: too many energy levels specified ", phan_flname);
 124             }
 125             poly = 1;
 126             fprintf(pstream, "SPECTRUM POLYCHROMATIC %i\n", nerg);
 127         }
 128     }
 129     if (poly)
 130         {
 131         if (fgets(string, sizeof (string), phanfl) == NULL)
 132                 {
 133             errorc("error: premature end of file encountered ", phan_flname);
 134         }
 135         while ( isBlank(string) ) //skip blank lines
 136             if ( fgets(string, sizeof (string), phanfl) == NULL) break;
 137
 138                 fprintf(pstream, "%s", string);
 139     }
 140
 141     fprintf(pstream, "OBJECTS\n");
 142
 143
 144     /* jk 01/01/08
 145      * scan the file to count the structures, then
 146      * allocate two arrays: structure[][] and pairs[] to hold the info about
 147      * the structures
 148      * this allows for "unlimited" number of structures in the phantom
 149      */
 150     while ((fgets(string, sizeof (string), phanfl)) != NULL)
 151         {
 152
 153         if (    strncasecmp(string, "dumm", 4) == 0 ||
 154                 strncasecmp(string, "sing", 4) == 0 ||
 155                 strncasecmp(string, "pair", 4) == 0 ||
 156                 strncasecmp(string, "scal", 4) == 0)
 157                 {
 158
 159                         /* keyword is 'dummy','pair', 'single' or 'scale' */
 160             strncpy(keyword, string, 4);
 161         }
 162
 163                 else if ( isSkip(string) ) //the line is a comment or blank, skip it
 164                 {;}
 165
 166                 else if (strncasecmp(keyword, "pair", 4) == 0)
 167                 {
 168             tot_nstr = tot_nstr + 2;
 169         }
 170
 171                 else if (strncasecmp(keyword, "sing", 4) == 0)
 172                 {
 173             tot_nstr = tot_nstr + 1;
 174         }
 175         /*else do nothing - DUMMY structures do not count */
 176     }
 177
 178     //allocate the arrays
 179     *pairs = calloc(tot_nstr, sizeof (int));
 180
 181     *structure = calloc(tot_nstr, sizeof (char*));
 182
 183     for (i = 0; i < tot_nstr; i++)
 184         {
 185         (*structure)[i] = calloc(MAXLINESIZE, sizeof (char));
 186     }
 187
 188     //set pointer to the file at its beginning
 189     rewind(phanfl);
 190
 191     while ((fgets(string, sizeof (string), phanfl)) != NULL)
 192         {
 193
 194         if (    strncasecmp(string, "dumm", 4) == 0 ||
 195                 strncasecmp(string, "sing", 4) == 0 ||
 196                 strncasecmp(string, "pair", 4) == 0 ||
 197                 strncasecmp(string, "scal", 4) == 0 ||
 198                 strncasecmp(string, "spec", 4) == 0)
 199                 {
 200             /* keyword is 'dummy','pair', */
 201             strncpy(keyword, string, 4); /* 'single' or 'scale'  or 'spec'*/
 202         }
 203
 204                 else if ( isComment(string) ) //the line is a comment, echo it to pstream
 205                 {
 206                         fprintf(pstream, "%s", string);
 207                 }
 208
 209                 else if ( isBlank(string) ) //the line is blank, skip it
 210                 {;}
 211
 212         else if (strncasecmp(keyword, "spec", 4) == 0)
 213                 {
 214             //ignore the next line after spec, if any; if was written to pstream before
 215         }
 216
 217                 else if (strncasecmp(keyword, "pair", 4) == 0)
 218                 {
 219                          /* if a paired structure:
 220                          (a) increment number of pairs read in;
 221                          (b) read in description of pair;
 222                          (c) errorchk2: make sure paired objects have unequal densities;
 223                          (d) randomly assign densities and save first structure;
 224                          (e) print out structure description onto SNARK input stream;
 225                          (f) save second structure and also print this structure's
 226                          description onto SNARK input stream.
 227                          */
 228
 229             if (strlen(string) - 1 != blanks(string) && isValidObject(string) )         /* non-blank line */
 230             {
 231                                 (*num_pairs)++;
 232                 sscanf(string, "%s%f%f%f%f%f%f%f%f", shape, &cx, &cy, &u, &v,
 233                         &ang, &den1, &den2, &prob);
 234                 errorchk2(den1, den2, prob, iroi_used);
 235                 if (poly)  //read the next two lines with polychromatic densities
 236                 {
 237                                         fgets(string, sizeof (string), phanfl);
 238                     while ( isSkip(string) ) //it's a comment line,
 239                     {                    //echo it to output and read next line
 240                         if( isComment(string) )     fprintf(pstream, "%s", string);
 241                         if ( fgets(string, sizeof (string), phanfl) == NULL) break;
 242                     }
 243
 244                     fgets(string_alt, sizeof (string), phanfl);
 245                     while ( isSkip(string) ) //it's a comment line,
 246                     {                    //echo it to output and read next line
 247                         if( isComment(string) )     fprintf(pstream, "%s", string);
 248                         if ( fgets(string, sizeof (string), phanfl) == NULL) break;
 249                     }
 250                 }
 251                 change = 0;
 252                 if (drand48() > prob)
 253                                 {
 254                     tmp = den1;
 255                     den1 = den2;
 256                     den2 = tmp;
 257                     change = 1;
 258                 }
 259                 sprintf((*structure)[nstr], "%s %f %f %f %f %f", shape, cx, cy, u, v, ang);
 260                 (*pairs)[npair++] = nstr + 1;
 261                 fprintf(pstream, "%s %f\n", (*structure)[nstr++], den1);
 262                 if (poly)  //print the remaining densities for polychromatic
 263                 {
 264                                         if (change) fprintf(pstream, "%s", string_alt);
 265                     else fprintf(pstream, "%s", string);
 266                 }
 267                 sprintf((*structure)[nstr], "%s %f %f %f %f %f", shape, -cx, cy, u, v, -ang);
 268                 (*pairs)[npair++] = nstr + 1;
 269                 fprintf(pstream, "%s %f\n", (*structure)[nstr++], den2);
 270                 if (poly)  //print the remaining densities for polychromatic
 271                 {
 272                                         if (change) fprintf(pstream, "%s", string);
 273                     else fprintf(pstream, "%s", string_alt);
 274                 }
 275             }
 276                         else errorc("Wrong command string in ", phan_flname);
 277         }
 278                 else if (strncasecmp(keyword, "sing", 4) == 0)
 279                 {
 280                         /* if a single structure then print out its description onto
 281                         SNARK's input stream. Also save this structure's description */
 282
 283             if (strlen(string) - 1 != blanks(string) && isValidObject(string) )          /* non-blank line */
 284             {
 285                                 fprintf(pstream, "%s", string);
 286                 sscanf(string, "%s%f%f%f%f%f", shape, &cx, &cy, &u, &v, &ang);
 287                 sprintf((*structure)[nstr++], "%s %f %f %f %f %f", shape, cx, cy, u, v, ang);
 288                 if (poly)  //read and print the next  line with polychromatic densities
 289                 {
 290                                         fgets(string, sizeof (string), phanfl);
 291                     while ( isSkip(string) ) //it's a comment line,
 292                     {                    //echo it to output and read next line
 293                         if( isComment(string) )     fprintf(pstream, "%s", string);
 294                         if ( fgets(string, sizeof (string), phanfl) == NULL) break;
 295                     }
 296                     fprintf(pstream, "%s", string);
 297                 }
 298             }
 299                         else errorc("Wrong command string in ", phan_flname);
 300         }/* for 'scale' and 'dummy' merely echo the input */
 301
 302                 else if (strncasecmp(keyword, "scal", 4) == 0)
 303                 {
 304                         tmp = sscanf(string, "%f%f", &scale, &sd);
 305             if (tmp == 2)   //print all: scale, phantom_seed and sd
 306             {
 307                                 if (sd <= 0)  errorc("Value of standard deviation cannot be negative "
 308                                            "in SCALE followup line in ", phan_flname);
 309                 fprintf(pstream, "last %f %ld %f \n", scale, phantom_seed, sd);
 310             }
 311             else if (tmp == 1)   //print only scale
 312                 fprintf(pstream, "last %f \n", scale);
 313             else
 314                 errorc("Invalid SCALE followup line in ", phan_flname);
 315             scaleCount++;
 316         }
 317                 else if (strncasecmp(keyword, "dumm", 4) == 0)
 318                 {
 319             if (strlen(string) - 1 != blanks(string) && isValidObject(string) )
 320                         {
 321                 fprintf(pstream, "%s", string);
 322                 if (poly) { //read and print the next line with polychromatic densities
 323                     fgets(string, sizeof (string), phanfl);
 324                     while ( isSkip(string) ) //it's a comment line,
 325                     {                    //echo it to output and read next line
 326                         if( isComment(string) )     fprintf(pstream, "%s", string);
 327                         if ( fgets(string, sizeof (string), phanfl) == NULL) break;
 328                     }
 329                     fprintf(pstream, "%s", string);
 330                 }
 331             }
 332                         else errorc("Wrong command string in ", phan_flname);
 333
 334         }
 335                 else errorc("Wrong command string in ", phan_flname);
 336
 337     }
 338
 339     if( scaleCount != 1)
 340         errorc("SCALE and its followup line should appear exactly once in ", phan_flname);
 341
 342     *num_str = nstr;
 343     fclose(phanfl);
 344     fprintf(pstream, "phantom average %d\n", navel);
 345     fprintf(pstream, "%d %f \n", nelem, pixsz);
 346 }
 347
 348 /*
 349  * This function checks to make sure that PAIRed structures do
 350  * not have the same density. Also checks to ensure that the probability
 351  * flag is a valid number
 352  */
 353
 354 void errorchk2(float den1, float den2, float prob, int iroi_used)
 355 {
 356     if (den1 == den2)
 357         {
 358         errorc("error: ", "paired structures with equal densities");
 359     }
 360
 361     if (prob > 1 || prob < 0)
 362         {
 363         errorc("error: ", "invalid probability for paired structures");
 364     }
 365
 366     //jk 09/07/2008
 367     //added a check to verify that density of one out of each
 368     //paired  structures is zero
 369     if (iroi_used)
 370         {
 371         if (den1 != 0)
 372                 {
 373             if (den2 != 0)
 374                         {
 375                 errorc("error: ", "when IROI is used, one of the paired structures "
 376                         "has to have density of zero, the other non-zero");
 377             }
 378         }
 379                 else
 380                 {
 381             if (den2 != 0)
 382                         {
 383                         //that's fine
 384             }
 385                         else
 386                         {
 387                 errorc("error: ", "when IROI is used, one of the paired structures "
 388                         "has to have density of zero, the other non-zero");
 389             }
 390
 391         }
 392
 393     }
 394
 395 }