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/eval.cpp $
5 $LastChangedRevision: 179 $
6 $Date: 2015-08-03 20:13:24 -0400 (Mon, 03 Aug 2015) $
8 ***********************************************************
10 EVAL IS THE POST-PROCESSING ROUTINE WHICH DRIVES THE TWO
11 RECONSTRUCTION EVALUATION PROCEDURES, DIST AND POINT.
12 THE PARAMETER "TYPE" DETERMINES WHICH OF THESE ROUTINES TO
14 IF TYPE=1 THEN DO DIST ONLY.
15 IF TYPE=2 THEN DO POINT ONLY.
16 OTHERWISE, DO BOTH DIST AND POINT.
17 MODIFIED TO COMPUTE THE RELATIVE ERROR AS DEFINED BY RAMA+LAK.
18 ALSO, MODIFIED SO THAT THE RELATIVE ERROR AND DISTANCE
19 MEASURES ARE NOT NORMALIZED IF ABSSUM AND TSTD ARE
22 Modified 4/1/88 to enable evaluation of subregions, either
23 through boundary selection by geomtric descriptions OR (OR BOTH)
24 pixel selection through density windowing wrt to phantom.
25 Only pixels within the specified regions and density window
26 will be used in computing the metrics.
27 A switch ISW is added to indicate whether the phantom is passed through
28 the first array, a(), or the second array, b(), so that density
29 selection can always be referenced wrt the phantom.
47 void Eval_class::eval_1()
59 static const CHAR * objname[] = { "wholepic", "rect", "elip" };
77 unsigned INTEGER count;
78 unsigned INTEGER iter;
86 BOOLEAN eol, linef, linef_set = false;
87 BOOLEAN fresid, fklds, fwsqd, fulimg, flagdw;
91 static bool called = false;
93 static const INTEGER eval_codes[9] =
95 CHAR2INT('r', 'e', 's', 'o'),
96 CHAR2INT('p', 'o', 'i', 'n'),
97 CHAR2INT('l', 'i', 'n', 'e'),
98 CHAR2INT('s', 't', 'r', 'i'),
99 CHAR2INT('w', 'h', 'o', 'l'),
100 CHAR2INT('s', 'e', 'l', 'e'),
101 CHAR2INT('e', 'l', 'i', 'p'),
102 CHAR2INT('r', 'e', 'c', 't'),
103 CHAR2INT('l', 'a', 's', 't')
107 // Read evaluation type specifications
109 word = InFile.getwrd(FALSE, &eol, eval_codes, 4);
116 if (word == eval_codes[0] || word == eval_codes[1])
118 type = (word == eval_codes[0]) ? 1 : 2;
119 word = InFile.getwrd(FALSE, &eol, &(eval_codes[2]), 2);
126 if (word == eval_codes[2])
131 if (word == eval_codes[3])
138 // Read evaluation name
139 InFile.GetComment(evalnam);
140 fprintf(output, "\n %s", evalnam);
142 // Read evaluation region specifications
143 word = InFile.getwrd(TRUE, &eol, &(eval_codes[4]), 2);
147 fprintf(output, "\n **** no region (WHOLEPIC or SELECTIVE) specified");
148 fprintf(output, "\n **** evaluation aborted\n");
152 if (word == eval_codes[4])
154 // Whole picture is specified
155 ObjN = 1; // set number of objects to 1
156 Obj = new Object[ObjN];
157 getden(&flagdw, 0); // get densities
158 Obj[0].objcod = 0; // set object type
159 fulimg = TRUE; // set full immage flag
163 // Selective is specified
168 eval_Object_Node *f_object, *t_object;
169 f_object = new eval_Object_Node();
174 word = InFile.getwrd(TRUE, &eol, &(eval_codes[6]), 3);
178 fprintf(output, "\n **** missing keyword elip, rect or last");
179 fprintf(output, "\n **** program aborted\n");
183 if (word == eval_codes[8])
185 // keyword last found
188 fprintf(output, "\n **** no region specified for SELECTIVE");
189 fprintf(output, "\n **** program aborted\n");
193 t_object->data.objcod = (word == eval_codes[7]) ? 1 : 2;
196 t_object->data.cx = InFile.getnum(FALSE, &eol);
199 fprintf(output, "\n **** parameter cx of object %s is missing",int2str(word));
200 fprintf(output, "\n **** program aborted\n");
205 t_object->data.cy = InFile.getnum(FALSE, &eol);
208 fprintf(output, "\n **** parameter cy of object %s is missing",int2str(word));
209 fprintf(output, "\n **** program aborted\n");
214 t_object->data.u = InFile.getnum(FALSE, &eol);
217 fprintf(output, "\n **** parameter u of object %s is missing",int2str(word));
218 fprintf(output, "\n **** program aborted\n");
223 t_object->data.v = InFile.getnum(FALSE, &eol);
226 fprintf(output, "\n **** parameter v of object %s is missing",int2str(word));
227 fprintf(output, "\n **** program aborted\n");
232 t_object->data.ang = InFile.getnum(FALSE, &eol);
235 fprintf(output, "\n **** parameter ang of object %s is missing",int2str(word));
236 fprintf(output, "\n **** program aborted\n");
240 getdens(&flagdw, t_object);
242 t_object->next = new eval_Object_Node();
243 t_object = t_object->next;
248 // copy list elements to array Obj
254 Obj = new Object[ObjN];
255 objdon = new BOOLEAN[ObjN];
259 for (i = 0; i < ObjN; i++)
261 Obj[i].cx = t_object->data.cx;
262 Obj[i].cy = t_object->data.cy;
263 Obj[i].u = t_object->data.u;
264 Obj[i].v = t_object->data.v;
265 Obj[i].ang = t_object->data.ang;
266 Obj[i].objcod = t_object->data.objcod;
267 Obj[i].t1 = t_object->data.t1;
268 Obj[i].t2 = t_object->data.t2;
270 t_object = t_object->next;
279 ave = new REAL[ObjN];
280 dis = new REAL[ObjN];
281 rel = new REAL[ObjN];
282 var = new REAL[ObjN];
283 std = new REAL[ObjN];
284 tstd = new REAL[ObjN];
285 tvar = new REAL[ObjN];
286 abssum = new REAL[ObjN];
288 // Echo region specifications
290 fprintf(output,"\n Region cx cy u v ang t1 t2\n");
292 if (Obj[0].objcod == 0)
294 fprintf(output,"\n %s %9.2g %9.2g",objname[0], Obj[0].t1, Obj[0].t2);
298 for (j = 0; j < ObjN; j++)
300 fprintf(output,"\n %2i %s %9.2f %9.2f %9.2f %9.2g %9.2g %9.2g %9.2g", j,objname[Obj[j].objcod], Obj[j].cx, Obj[j].cy, Obj[j].u,Obj[j].v, Obj[j].ang, Obj[j].t1, Obj[j].t2);
304 // READ AND ECHO FLAGS
311 if(flag[0]==flag[i]||flag[i]==0);
313 fprintf(output, "\n **** The value of fl[0] is %d, but the value of fl[%d] is %d ",flag[0],i,flag[i]);
314 fprintf(output, "\n **** program aborted\n");
318 if(flag[0]==2){fresid = TRUE;}
319 if(flag[0]==3){fklds = TRUE; fwsqd = TRUE; fresid = TRUE;}
321 /* for (i = 0; i < 51; i++)
330 else if (flag[i] >= 2)
337 // RESIDUAL CANNOT AND WILL NOT BE CALCULATED IF SOME PIXELS
338 // ARE EXCLUDED BY DENSITY WINDOWING OR REGION SELECTION.
339 // RESET RESIDUAL FLAG
353 flag0 = (fresid) ? 2 : 1;
359 // check if called before
362 eval1 = fopen("eval", "w"); // open new file
367 eval1 = fopen("eval", "a+"); // open existing file for appending
370 // print evaluation name
371 fprintf(eval1, "\n\nevaluation name: %s\n", evalnam);
374 if (RecFile.Open("recfil") != 0)
376 fprintf(output, "\n **** unable to open recfil");
377 fprintf(output, "\n **** EVAL execution aborted\n");
382 // Get projection name
383 if (RecFile.GetProjName(prjnam) != 0)
388 // get picture dimension
389 if (RecFile.GetNelem(&(GeoPar.nelem)) != 0)
394 // get picture dimension
395 if (RecFile.GetPixSize(&(GeoPar.pixsiz)) != 0)
400 // calculate picture area
401 GeoPar.area = GeoPar.nelem * GeoPar.nelem;
402 GeoPar.midpix = (GeoPar.nelem + 1) / 2;
404 // allocate phantom and reconstruction arrays
405 Phantom = new REAL[GeoPar.area];
406 Recon = new REAL[GeoPar.area];
409 if (RecFile.ReadPhan(phnnam, Phantom) != 0)
411 fprintf(output, "\n **** phantom not present");
412 fprintf(output, "\n **** EVAL execution aborted\n");
419 // calculate resolution
421 // open prjfil only if residual is calculated
424 if (ProjFile.Open("prjfil") != 0)
426 fprintf(output, "\n **** unable to open prjfil");
427 fprintf(output, "\n **** EVAL execution aborted\n");
432 ProjFile.ReadHeader(prjnam, &NoisePar, &Spctrm, &Anglst);
434 if (!linef_set) linef = GeoPar.line;
436 bldlst(&list, &weight);
439 // COMPUTE AVERAGE, VARIANCE AND STANDARD DEVIATION OF TEST PATTERN
441 for (j = 0; j < GeoPar.area; j++)
444 Blob.pix_basis = TRUE;
445 dist(1, Recon, Phantom, linef, list, weight, flag0, ave, dis, abssum, tvar, tstd, &resid, fresid, &klds, fklds, &wsqd, fwsqd);
448 fprintf(eval1, "\n global resolution measures\n");
450 fprintf(eval1, "\n phantom name: %s", phnnam);
452 fprintf(eval1, "\n projection name: %s", prjnam);
458 fprintf(output,"\n residual calculation using line integration");
463 fprintf(output,"\n residual calculation using strip integration");
467 fprintf(eval1, "\n\n metrics for test phantom");
469 if (fresid && (fklds || fwsqd))
471 fprintf(eval1,"\n region area average distance rel err variance std dev residual klds wsqd");
475 fprintf(eval1,"\n region area average distance rel err variance std dev residual");
479 fprintf(eval1,"\n region area average distance rel err variance std dev");
482 for (j = 0; j < ObjN; j++)
484 Obj[j].avep = ave[j];
486 if (fresid && (fklds || fwsqd))
488 fprintf(eval1,"\n %6d %7d %9.4f %9.4f %9.4f %10.4f %12.4f %12.4f",j, Obj[j].npt, Obj[j].avep, tvar[j], tstd[j], resid, klds, wsqd);
492 fprintf(eval1,"\n %6d %7d %9.4f %9.4f %9.4f %10.4f",j, Obj[j].npt, Obj[j].avep, tvar[j], tstd[j], resid);
496 fprintf(eval1,"\n %6d %7d %9.4f %9.4f %9.4f",j, Obj[j].npt, Obj[j].avep, tvar[j], tstd[j]);
500 while (RecFile.ReadRec(recnam, algn, &count, &iter, Recon) == 0)
503 // if first iteration of algorithm print header
506 fprintf(eval1, "\n\n execution name: %s", recnam);
510 fprintf(eval1, "\n metrics for algorithm %s", algn);
512 if (fresid && (fklds || fwsqd))
514 fprintf(eval1,"\n iter area average distance rel err variance std dev residual klds wsqd");
518 fprintf(eval1,"\n iter area average distance rel err variance std dev residual");
522 fprintf(eval1,"\n iter area average distance rel err variance std dev");
527 // if iteration not selected continue
528 if (flag[iter - 1] == 0)
535 fprintf(eval1, "\n metrics for algorithm %4s iteration %4d",algn, count);
537 fprintf(eval1, "\n region area average distance rel err variance std dev");
539 temp = Blob.pix_basis;
540 Blob.pix_basis = TRUE;
541 dist(0, Phantom, Recon, linef, list, weight, flag[iter - 1], ave, dis, rel, var, std, &resid, fresid, &klds, fklds, &wsqd, fwsqd);
542 Blob.pix_basis = temp;
543 // Loop for each regions
544 for (j = 0; j < ObjN; j++)
546 if (abssum[j] > Consts.zero)
551 if (tstd[j] > Consts.zero)
556 if (((flag[iter - 1] == 1) || (flag[iter - 1] > 1)) && (!fresid))
560 fprintf(eval1,"\n %4d %7d %9.4f %9.4f %9.4f %9.4f %10.4f", count, Obj[j].npt, ave[j], dis[j], rel[j], var[j], std[j]);
564 fprintf(eval1,"\n %6d %7d %9.4f %9.4f %9.4f %9.4f %10.4f", j, Obj[j].npt, ave[j], dis[j], rel[j], var[j], std[j]);
568 if ((flag[iter - 1] > 1) && fresid && (fklds || fwsqd))
570 fprintf(eval1,"\n %4d %7d %9.4f %9.4f %9.4f %9.4f %9.4f %10.4f %12.4f %12.4f", count, Obj[j].npt, ave[j], dis[j], rel[j], var[j], std[j], resid, klds, wsqd);
572 else if ((flag[iter - 1] > 1) && fresid)
574 fprintf(eval1,"\n %4d %7d %9.4f %9.4f %9.4f %9.4f %9.4f %10.4f", count, Obj[j].npt, ave[j], dis[j], rel[j], var[j], std[j], resid);
580 fprintf(eval1, "\n");
590 // calculate point measurements
591 // if type was 1 we are done
592 if (type == 1 || !fulimg || (fulimg && flagdw))
600 fprintf(output,"\n **** Selective was specified, point measures will not be computed.\n");
604 fprintf(output,"\n **** Density range was specified, point measures will not be computed.\n");
609 { // if type != 1 -> type == 0 (both) or type = 2
611 fprintf(eval1, "\n point by point resolution measures\n");
613 fprintf(eval1, "\n phantom name: %s", phnnam);
615 if (fresid || (type != 2))
617 fprintf(eval1, "\n projection name: %s", prjnam);
620 // select first reconstruction in file
621 RecFile.SelectRec(0);
623 // loop though receonstructions
624 while (RecFile.ReadRec(recnam, algn, &count, &iter, Recon) == 0)
629 // if first iteration of algorithm print header
635 // if iteration not selected continue
636 if (flag[iter - 1] == 0)
641 point(Phantom, Recon, Recon, GeoPar.nelem, &nerr, epsil);
645 fprintf(eval1, "\n\n execution name: %s", recnam);
647 fprintf(eval1, "\n algn iter"); // only print out 'e(n)' as many times as there
648 for (n = 0; n < nerr; n++)
649 { // are values in 'epsil[]' (i.e. 'nerr' times)
650 fprintf(eval1, " e(%d)", n);
656 fprintf(eval1, "\n %4s %6d", algn, count);
657 for (n = 0; n < nerr; n++)
659 fprintf(eval1, " %9.4f", epsil[n]);
683 fprintf(eval1, "\n ");
692 if (list != NULL) delete[] list;
694 if (weight != NULL) delete[] weight;