src/snark/imagewise_roi.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/imagewise_roi.c $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  Previously part of SuperSNARK
  11  */
  12
  13 #include <assert.h>
  14
  15 #include <stdlib.h>
  16 #include "experimenter.h"
  17 #include "read_eval_phantom1.h"
  18 #include "read_eval_recon1.h"
  19 #include "hit_ratio.h"
  20
  21 /* -----------------------------imagewise_roi.c -------------------------------
  22  //jk 01/11/2008     new fom introduced based on Narayan & Herman paper
  23
  24  This function computes theimagewise-roi.
  25
  26  INPUTS:
  27  itr1 - array of length niters containing the iteration numbers for the
  28  first algorithm.
  29  itr2 - array of length niters containing the iteration numbers for the
  30  second algorithm.
  31  niters - number of iterations to be compared.
  32  pairs - array of length numpairs containing the indices of the structures
  33  which are paired.
  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.
  37
  38  OUTPUTS:
  39  imagewise_roi1 - array of length niters containing the hit ratio for
  40  the first algorithm.
  41  imagewise_roi2 - array of length niters containing the hit ratio for
  42  the second algorithm.
  43
  44  */
  45
  46 void imagewise_roi(int* itr1, int* itr2, int niters, int* pairs, int numpairs, char* keywrd1, char* keywrd2, double* imagewise_roi1, double* imagewise_roi2)
  47 {
  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;
  54
  55         /*********** compute goodness of the phantom ****************************************/
  56
  57         //read in abnormality index for phantom structures
  58         read_eval_phantom1(&regions, &numstr, &phantom, &strarea); //regions is not used anywhere later on
  59
  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));
  64
  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)
  70         {
  71                 if (phantom[pairs[i] - 1] > phantom[pairs[i + 1] - 1])
  72                 {
  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]]);
  76                 }
  77                 else
  78                 {
  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]]);
  82                 }
  83
  84                 sumnsph += phantom[nonsignal[i / 2]];
  85                 sumsph += phantom[signal[i / 2]];
  86
  87         }
  88
  89         //compute the average density of non-signal containing objects
  90         avnsph = sumnsph / numpairs;
  91
  92         //compute standard deviation of density of non-signal containing objects
  93         for (i = 0; i < numpairs; i++)
  94         {
  95                 sdnsph += (phantom[nonsignal[i]] - avnsph)
  96                                 * (phantom[nonsignal[i]] - avnsph);
  97         }
  98
  99         sdnsph = sqrt(sdnsph);
 100         goodness_ph = (sumsph - sumnsph) / sdnsph;
 101
 102         /*********** compute goodness of the reconstructions ********************************/
 103
 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));
 108
 109         for (i = 0; i < niters; i++)
 110         {
 111                 /* read in abnormality indexes for structures in
 112                  the reconstructions. Initialize the arrays
 113                  which will store the structural accuracy */
 114
 115                 sumnsrecon1 = 0;
 116                 sumnsrecon2 = 0;
 117                 sumsrecon1 = 0;
 118                 sumsrecon2 = 0;
 119                 avnsrecon1 = 0;
 120                 avnsrecon2 = 0;
 121                 sdnsrecon1 = 0;
 122                 sdnsrecon2 = 0;
 123                 goodness_recon1 = 0;
 124                 goodness_recon2 = 0;
 125
 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;
 130
 131                 //compute sumnsrecon1 and sumsrecon1
 132                 //and  sumnsrecon2 and sumsrecon2
 133                 for (j = 0; j < numpairs; j++)
 134                 {
 135                         sumnsrecon1 += recon1[nonsignal[j]];
 136                         sumsrecon1 += recon1[signal[j]];
 137
 138                         sumnsrecon2 += recon2[nonsignal[j]];
 139                         sumsrecon2 += recon2[signal[j]];
 140
 141                 }
 142
 143                 //compute the average density of non-signal containing objects
 144                 avnsrecon1 = sumnsrecon1 / numpairs;
 145                 avnsrecon2 = sumnsrecon2 / numpairs;
 146
 147                 //compute standard deviation of density of non-signal containing objects
 148                 for (j = 0; j < numpairs; j++)
 149                 {
 150                         sdnsrecon1 += (recon1[nonsignal[j]] - avnsrecon1)
 151                                         * (recon1[nonsignal[j]] - avnsrecon1);
 152                         sdnsrecon2 += (recon2[nonsignal[j]] - avnsrecon2)
 153                                         * (recon2[nonsignal[j]] - avnsrecon2);
 154                 }
 155
 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;
 163         }
 164
 165         free(phantom); /* free memory for next function call */
 166         free(signal);
 167         free(nonsignal);
 168         free(recon1);
 169         free(recon2);
 170         free(strarea);
 171
 172 }