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/imagewise_roi.c $
5 $LastChangedRevision: 80 $
6 $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
8 ***********************************************************
10 Previously part of SuperSNARK
16 #include "experimenter.h"
17 #include "read_eval_phantom1.h"
18 #include "read_eval_recon1.h"
19 #include "hit_ratio.h"
21 /* -----------------------------imagewise_roi.c -------------------------------
22 //jk 01/11/2008 new fom introduced based on Narayan & Herman paper
24 This function computes theimagewise-roi.
27 itr1 - array of length niters containing the iteration numbers for the
29 itr2 - array of length niters containing the iteration numbers for the
31 niters - number of iterations to be compared.
32 pairs - array of length numpairs containing the indices of the structures
34 numpairs - number of pairs of structures.
35 keywrd1 - string containing the keyword which defines the first algorithm.
36 keywrd2 - string containing the keyword which defines the second algorithm.
39 imagewise_roi1 - array of length niters containing the hit ratio for
41 imagewise_roi2 - array of length niters containing the hit ratio for
46 void imagewise_roi(int* itr1, int* itr2, int niters, int* pairs, int numpairs, char* keywrd1, char* keywrd2, double* imagewise_roi1, double* imagewise_roi2)
48 double *phantom, *recon1, *recon2;
49 int *regions, numstr, i, j, *strarea, *signal, *nonsignal;
50 double sumnsph = 0, sumnsrecon1 = 0, sumnsrecon2 = 0, sumsph = 0,
51 sumsrecon1 = 0, sumsrecon2 = 0, avnsph = 0, avnsrecon1 = 0,
52 avnsrecon2 = 0, sdnsph = 0, sdnsrecon1 = 0, sdnsrecon2 = 0,
53 goodness_ph = 0, goodness_recon1 = 0, goodness_recon2 = 0;
55 /*********** compute goodness of the phantom ****************************************/
57 //read in abnormality index for phantom structures
58 read_eval_phantom1(®ions, &numstr, &phantom, &strarea); //regions is not used anywhere later on
60 //allocate memmory to store indecis of structures that contain signal and
61 //the ones that do not contain signal
62 signal = (int*) malloc(numpairs * sizeof(int));
63 nonsignal = (int*) malloc(numpairs * sizeof(int));
65 //scan the phantom array to identify the structures that contain signal
66 //ASSUMPTION: the higher density is associated with the signal.
67 //the lower with non-signal structure
68 //compute sumnsph and sumsph
69 for (i = 0; i < numpairs * 2; i += 2)
71 if (phantom[pairs[i] - 1] > phantom[pairs[i + 1] - 1])
73 signal[i / 2] = pairs[i] - 1;
74 nonsignal[i / 2] = pairs[i + 1] - 1;
75 //assert(phantom[signal[i/2]] > phantom[nonsignal[i/2]]);
79 signal[i / 2] = pairs[i + 1] - 1;
80 nonsignal[i / 2] = pairs[i] - 1;
81 //assert(phantom[signal[i/2]] > phantom[nonsignal[i/2]]);
84 sumnsph += phantom[nonsignal[i / 2]];
85 sumsph += phantom[signal[i / 2]];
89 //compute the average density of non-signal containing objects
90 avnsph = sumnsph / numpairs;
92 //compute standard deviation of density of non-signal containing objects
93 for (i = 0; i < numpairs; i++)
95 sdnsph += (phantom[nonsignal[i]] - avnsph)
96 * (phantom[nonsignal[i]] - avnsph);
99 sdnsph = sqrt(sdnsph);
100 goodness_ph = (sumsph - sumnsph) / sdnsph;
102 /*********** compute goodness of the reconstructions ********************************/
104 //allocate memory to store abnormality indexes for
105 //structures in the reconstructions
106 recon1 = (double *) malloc(numstr * sizeof(double));
107 recon2 = (double *) malloc(numstr * sizeof(double));
109 for (i = 0; i < niters; i++)
111 /* read in abnormality indexes for structures in
112 the reconstructions. Initialize the arrays
113 which will store the structural accuracy */
126 read_eval_recon1(&itr1[i], keywrd1, numstr, recon1);
127 read_eval_recon1(&itr2[i], keywrd2, numstr, recon2);
128 imagewise_roi1[i] = 0.0;
129 imagewise_roi2[i] = 0.0;
131 //compute sumnsrecon1 and sumsrecon1
132 //and sumnsrecon2 and sumsrecon2
133 for (j = 0; j < numpairs; j++)
135 sumnsrecon1 += recon1[nonsignal[j]];
136 sumsrecon1 += recon1[signal[j]];
138 sumnsrecon2 += recon2[nonsignal[j]];
139 sumsrecon2 += recon2[signal[j]];
143 //compute the average density of non-signal containing objects
144 avnsrecon1 = sumnsrecon1 / numpairs;
145 avnsrecon2 = sumnsrecon2 / numpairs;
147 //compute standard deviation of density of non-signal containing objects
148 for (j = 0; j < numpairs; j++)
150 sdnsrecon1 += (recon1[nonsignal[j]] - avnsrecon1)
151 * (recon1[nonsignal[j]] - avnsrecon1);
152 sdnsrecon2 += (recon2[nonsignal[j]] - avnsrecon2)
153 * (recon2[nonsignal[j]] - avnsrecon2);
156 sdnsrecon1 = sqrt(sdnsrecon1);
157 sdnsrecon2 = sqrt(sdnsrecon2);
158 goodness_recon1 = ((sumsrecon1 - sumnsrecon1) / sdnsrecon1);
159 goodness_recon2 = ((sumsrecon2 - sumnsrecon2) / sdnsrecon2);
160 //output the value of the normalized FOM
161 imagewise_roi1[i] = goodness_recon1 / goodness_ph;
162 imagewise_roi2[i] = goodness_recon2 / goodness_ph;
165 free(phantom); /* free memory for next function call */