src/snark/emap_CEvalMAP.cpp

   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/emap_CEvalMAP.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10
  11  ===========================================================================
  12  Evaluation of the MAP algorithms (the results are written to "user1")
  13  ===========================================================================
  14  */
  15
  16 #include <cstdio>
  17 #include <cmath>
  18
  19 #include "blkdta.h"
  20 #include "geom.h"
  21 #include "uiod.h"
  22 #include "wray.h"
  23 #include "blob.h"
  24
  25 #include "emap.h"
  26
  27 void emap_class::EvalMAP(REAL* recon, INTEGER iter)
  28 {
  29         //------ local variables
  30         INTEGER i, j, iadr, nb, np, nr, numb, list[4096], usnr;
  31         REAL snorm, weight[4096], diff1;
  32         double sum1, sum2, y, diffnbs, totalsum, sclf, smx, pictk, ssum1, ssum2;
  33
  34         //------ evaluation
  35
  36         //------ log-likelihood
  37         sum1 = 0.0;
  38         sclf = MAPEM.sclfactor;
  39         for (i = 0; i < GeoPar.prjnum; i++)
  40         {
  41                 ssum1 = 0.0;
  42                 for (j = 0; j < GeoPar.usrays; j++)
  43                 {
  44                         np = i;
  45                         usnr = j;
  46                         nr = usnr + GeoPar.fusray;
  47
  48                         y = MAPEM.prj[(np) * GeoPar.usrays + usnr] * sclf;
  49
  50                         smx = 0.0;
  51                         if (Blob.pix_basis)
  52                         {
  53                                 wray(np, nr, list, weight, &numb, &snorm);
  54                         }
  55                         else
  56                         {
  57                                 Blob.bwray(np, nr, list, weight, &numb, &snorm);
  58                         }
  59
  60                         if (numb != 0)
  61                         {
  62                                 for (nb = 0; nb < numb; nb++)
  63                                 {
  64                                         iadr = list[nb];
  65                                         pictk = recon[iadr];
  66                                         smx += pictk * weight[nb];
  67                                 }
  68                         }
  69                         smx *= sclf;
  70
  71                         if (smx > 0.0)
  72                         {
  73                                 ssum1 += y * log(smx) - smx;
  74                         }
  75                         else if (smx < 0.0)
  76                         {
  77                                 fprintf(output, "\n In EvalMAP  smx = %14.7f y = %14.7f", smx,
  78                                                 y);
  79                                 fprintf(output, "\n smx must be greater than 0.0");
  80                         }
  81
  82                 }
  83                 sum1 += ssum1;
  84         }
  85
  86         sum2 = 0.0; // moved out of the `if' block for initialization. Lajos, Feb 10, 2005
  87         if (Blob.pix_basis)
  88         { // when blob, gamma is 0.0, so this step is skept. hstau 6/25/2003
  89
  90                 for (i = GeoPar.nelem; i < GeoPar.area - 2 * GeoPar.nelem + 1; i +=
  91                                 GeoPar.nelem)
  92                 {
  93                         ssum2 = 0.0;
  94                         for (j = i + 1; j < i + GeoPar.nelem - 2; j++)
  95                         {
  96                                 diff1 = recon[j]
  97                                                 - (REAL) 0.125
  98                                                                 * (recon[j - 1 - GeoPar.nelem]
  99                                                                                 + recon[j - GeoPar.nelem]
 100                                                                                 + recon[j + 1 - GeoPar.nelem]
 101                                                                                 + recon[j - 1] + recon[j + 1]
 102                                                                                 + recon[j - 1 + GeoPar.nelem]
 103                                                                                 + recon[j + GeoPar.nelem]
 104                                                                                 + recon[j + 1 + GeoPar.nelem]);
 105                                 diffnbs = diff1 * MAPEM.sclfactor;
 106                                 ssum2 += diffnbs * diffnbs;
 107                         }
 108                         sum2 += ssum2;
 109                 }
 110
 111                 sum2 = -0.5 * MAPEM.gamma1 * sum2;
 112         }  // if pixel basis
 113
 114         totalsum = sum1 + sum2;
 115
 116         //------ write results
 117         if (iter == 1)
 118         {
 119                 if (MAPEM.gamma1 == 0.0)
 120                 {
 121                         fprintf(MAPUser1, "\n    iter     log-likelihood");
 122                 }
 123                 else
 124                 {
 125                         ///if(MAPEM.gamma1 != 0.0) {
 126                         fprintf(MAPUser1,
 127                                         "\n    iter     log-likelihood          log-posterior");
 128                 }
 129         }
 130
 131         if (MAPEM.gamma1 == 0.0)
 132         {
 133                 fprintf(MAPUser1, "\n %5i  %20.9f", iter, sum1);
 134         }
 135         else
 136         {
 137                 ///if(MAPEM.gamma1 != 0.0) {
 138                 fprintf(MAPUser1, "\n %5i  %20.9f  %20.9f", iter, sum1, totalsum);
 139         }
 140
 141         fflush(MAPUser1);
 142
 143         return;
 144 }