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) $
8 **********************************************************************
10 Previously part of SuperSNARK
13 /*--------------------------------------------------------------------------
15 * Written by Jingsheng Zheng... November 1993
16 * Modified by Jolyon Browne ... June 1993
17 * --------------------------------------------------------------------------*/
23 #include "experimenter.h"
26 #include "create_phantom.h"
28 /*--------------------------------------------------------------------------
29 * This fuction generates the input lines that defines the elemental
30 * objects as required by the CREATE command.
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).
40 * structure - string array of length num_str containing the elemental object
42 * num_str - number of structures in the phantom.
43 * num_pairs - for PAIRed structures, this is the number of pairs of
45 * pairs - array of length 2*num_pairs containing the indicies of all PAIRed
47 * ----------------------------------------------------------------------------*/
55 char*** structure, //jk 01/01/08 modified to remove
56 // restriction for number of structures
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
65 float cx, cy, u, v, ang, den1, den2, prob, tmp;
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).
74 char shape[10], keyword[4], string[MAXLINESIZE];
80 char spect[10], spectType[15];
84 char toc[14][MAXLINESIZE];
87 char string_alt[MAXLINESIZE];
88 long phantom_seed = (long) (100000000 * drand48() );
90 if ((phanfl = fopen(phan_flname, "r")) == NULL) {
91 errorc("In create_phantom.c: error in opening file ", phan_flname);
94 /* print out beginning of CREATE input sequence */
96 fprintf(pstream, "CREATE\n");
97 fprintf(pstream, "%s, seed used in experimenter to generate this file = %ld\n", phan_flname, seed);
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);
103 while ( isSkip(string) ) //it's a comment line,
104 { //echo it to output and read next line
105 if( isComment(string) )
107 fprintf(pstream, "%s", string);
108 if ( fgets(string, sizeof (string), phanfl) == NULL) break;
111 if (strncasecmp(string, "SPEC", 4) == 0)
113 sscanf(string, "%s%s%i", spect, spectType, &erg);
115 if (strncasecmp(spectType, "mono", 4) == 0)
116 fprintf(pstream, "SPECTRUM MONOCHROMATIC %i\n", erg);
118 else if (strncasecmp(spectType, "poly", 4) == 0)
123 errorc("error: too many energy levels specified ", phan_flname);
126 fprintf(pstream, "SPECTRUM POLYCHROMATIC %i\n", nerg);
131 if (fgets(string, sizeof (string), phanfl) == NULL)
133 errorc("error: premature end of file encountered ", phan_flname);
135 while ( isBlank(string) ) //skip blank lines
136 if ( fgets(string, sizeof (string), phanfl) == NULL) break;
138 fprintf(pstream, "%s", string);
141 fprintf(pstream, "OBJECTS\n");
145 * scan the file to count the structures, then
146 * allocate two arrays: structure[][] and pairs[] to hold the info about
148 * this allows for "unlimited" number of structures in the phantom
150 while ((fgets(string, sizeof (string), phanfl)) != NULL)
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)
159 /* keyword is 'dummy','pair', 'single' or 'scale' */
160 strncpy(keyword, string, 4);
163 else if ( isSkip(string) ) //the line is a comment or blank, skip it
166 else if (strncasecmp(keyword, "pair", 4) == 0)
168 tot_nstr = tot_nstr + 2;
171 else if (strncasecmp(keyword, "sing", 4) == 0)
173 tot_nstr = tot_nstr + 1;
175 /*else do nothing - DUMMY structures do not count */
178 //allocate the arrays
179 *pairs = calloc(tot_nstr, sizeof (int));
181 *structure = calloc(tot_nstr, sizeof (char*));
183 for (i = 0; i < tot_nstr; i++)
185 (*structure)[i] = calloc(MAXLINESIZE, sizeof (char));
188 //set pointer to the file at its beginning
191 while ((fgets(string, sizeof (string), phanfl)) != NULL)
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)
200 /* keyword is 'dummy','pair', */
201 strncpy(keyword, string, 4); /* 'single' or 'scale' or 'spec'*/
204 else if ( isComment(string) ) //the line is a comment, echo it to pstream
206 fprintf(pstream, "%s", string);
209 else if ( isBlank(string) ) //the line is blank, skip it
212 else if (strncasecmp(keyword, "spec", 4) == 0)
214 //ignore the next line after spec, if any; if was written to pstream before
217 else if (strncasecmp(keyword, "pair", 4) == 0)
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.
229 if (strlen(string) - 1 != blanks(string) && isValidObject(string) ) /* non-blank line */
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
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;
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;
252 if (drand48() > prob)
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
264 if (change) fprintf(pstream, "%s", string_alt);
265 else fprintf(pstream, "%s", string);
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
272 if (change) fprintf(pstream, "%s", string);
273 else fprintf(pstream, "%s", string_alt);
276 else errorc("Wrong command string in ", phan_flname);
278 else if (strncasecmp(keyword, "sing", 4) == 0)
280 /* if a single structure then print out its description onto
281 SNARK's input stream. Also save this structure's description */
283 if (strlen(string) - 1 != blanks(string) && isValidObject(string) ) /* non-blank line */
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
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;
296 fprintf(pstream, "%s", string);
299 else errorc("Wrong command string in ", phan_flname);
300 }/* for 'scale' and 'dummy' merely echo the input */
302 else if (strncasecmp(keyword, "scal", 4) == 0)
304 tmp = sscanf(string, "%f%f", &scale, &sd);
305 if (tmp == 2) //print all: scale, phantom_seed and sd
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);
311 else if (tmp == 1) //print only scale
312 fprintf(pstream, "last %f \n", scale);
314 errorc("Invalid SCALE followup line in ", phan_flname);
317 else if (strncasecmp(keyword, "dumm", 4) == 0)
319 if (strlen(string) - 1 != blanks(string) && isValidObject(string) )
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;
329 fprintf(pstream, "%s", string);
332 else errorc("Wrong command string in ", phan_flname);
335 else errorc("Wrong command string in ", phan_flname);
340 errorc("SCALE and its followup line should appear exactly once in ", phan_flname);
344 fprintf(pstream, "phantom average %d\n", navel);
345 fprintf(pstream, "%d %f \n", nelem, pixsz);
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
354 void errorchk2(float den1, float den2, float prob, int iroi_used)
358 errorc("error: ", "paired structures with equal densities");
361 if (prob > 1 || prob < 0)
363 errorc("error: ", "invalid probability for paired structures");
367 //added a check to verify that density of one out of each
368 //paired structures is zero
375 errorc("error: ", "when IROI is used, one of the paired structures "
376 "has to have density of zero, the other non-zero");
387 errorc("error: ", "when IROI is used, one of the paired structures "
388 "has to have density of zero, the other non-zero");