src/snark/emap.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.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  ===========================================================================
  11  maximum a-posteriori probability expectation maximization algorithm
  12  for the reconstruction in emission tomography (Herman, De Pierro, and
  13  Gai, "On Methods for Maximum A Posteriori Image Reconstruction with
  14  a Normal Prior," Technical Report MIPG181, Medical Image Processing
  15  Group, Dept. of Radiology, Univ. of Pennsylvania, 1992.)
  16
  17  input line in snark:
  18  gamma  [eval sy1 sy2]]
  19  e.g.:  1.0    eval 10  1
  20
  21  ===========================================================================
  22  */
  23
  24 #include <cstdio>
  25 #include <cstdlib>
  26 #include <cmath>
  27
  28 #include "blkdta.h"
  29 #include "geom.h"
  30 #include "consts.h"
  31 #include "blob.h"
  32 #include "uiod.h"
  33
  34 #include "emap.h"
  35
  36 INTEGER emap_class::Init()
  37 {
  38         MAPUser1 = fopen("MAPUser1", "r+");
  39         MAPEM.qj = NULL;
  40         MAPEM.Kx = NULL;
  41         MAPEM.Sx = NULL;
  42         MAPEM.Wjj = NULL;
  43         MAPEM.Sjjinv = NULL;
  44         MAPEM.prj = NULL;                       //wei 3/2005
  45
  46         return 0;
  47 }
  48
  49 INTEGER emap_class::Reset()
  50 {
  51         fclose(MAPUser1);
  52         if (MAPEM.qj != NULL)
  53                 delete[] MAPEM.qj;
  54         if (MAPEM.Kx != NULL)
  55                 delete[] MAPEM.Kx;
  56         if (MAPEM.Sx != NULL)
  57                 delete[] MAPEM.Sx;
  58         if (MAPEM.Wjj != NULL)
  59                 delete[] MAPEM.Wjj;
  60         if (MAPEM.Sjjinv != NULL)
  61                 delete[] MAPEM.Sjjinv;
  62         if (MAPEM.prj != NULL)
  63                 delete[] MAPEM.prj;                //wei 3/2005
  64         return 0;
  65 }
  66
  67 BOOLEAN emap_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
  68 {
  69         INTEGER area;
  70         BOOLEAN flag;  //added hstau 1/04
  71
  72         if (Blob.pix_basis)
  73         {
  74                 area = GeoPar.area;
  75         }
  76         else
  77         {
  78                 area = Blob.area;
  79         }
  80
  81 //------------------------   statements    ----------------------------
  82         if (iter == 1)
  83         {
  84
  85                 MAPUser1 = fopen("MAPUser1", "w");
  86
  87                 if (GeoPar.par)
  88                 {
  89                         MAPEM.delta = MAX0((REAL) 1.0,
  90                                         (REAL) 2.0 * (GeoPar.pixsiz / GeoPar.pinc));
  91                 }
  92                 else
  93                 {
  94                         MAPEM.delta = GeoPar.pixsiz / (GeoPar.pinc * (REAL) sqrt(2.) * (GeoPar.radius - (REAL) sqrt(.5) * GeoPar.pixsiz * GeoPar.nelem) / GeoPar.stod);
  95                 }
  96
  97                 MAPEM.delta++;
  98
  99                 if (((INTEGER) MAPEM.delta) > GeoPar.usrays)
 100                 {
 101                         MAPEM.delta = (REAL) GeoPar.usrays;
 102                 }
 103
 104                 // nelem not used. hstau
 105                 ReadInp(area, GeoPar.nelem, GeoPar.usrays, GeoPar.prjnum, &flag); //hstau 1/04
 106
 107                 if (flag)
 108                         return TRUE; //added hstau 1/04
 109
 110                 if (MAPEM.gamma1 <= Consts.zero)
 111                 {
 112                         MAPEM.gamma1 = 0.0;
 113                 }
 114
 115                 if (GeoPar.nelem == 1)
 116                         MAPEM.gamma1 = 0.0;                                  //wei 1/2005
 117
 118                 if (!Blob.pix_basis && MAPEM.gamma1 > 0.0)
 119                 {
 120                         fprintf(output,
 121                                         "\n *** error: gamma must be 0.0 when using blob basis");
 122                         return TRUE;  //it was exit(1); bug 123, hstau 7/8/05
 123                 }
 124
 125                 InitMAP(recon);
 126         }  // iter ==1
 127
 128         for (MAPEM.idy = 0; MAPEM.idy <= MAPEM.idymax - 1; MAPEM.idy++)
 129         {
 130                 for (MAPEM.idx = 0; MAPEM.idx <= MAPEM.idxmax - 1; MAPEM.idx++)
 131                 {
 132                         if (MAP(recon))
 133                                 return TRUE;            //added by wei, 1/2005
 134                 }
 135         }
 136
 137         if (MAPEM.mapevl)
 138         {
 139                 EvalMAP(recon, iter);
 140         }
 141
 142         return FALSE;
 143 }