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) $
8 ***********************************************************
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.)
21 ===========================================================================
36 INTEGER emap_class::Init()
38 MAPUser1 = fopen("MAPUser1", "r+");
44 MAPEM.prj = NULL; //wei 3/2005
49 INTEGER emap_class::Reset()
58 if (MAPEM.Wjj != NULL)
60 if (MAPEM.Sjjinv != NULL)
61 delete[] MAPEM.Sjjinv;
62 if (MAPEM.prj != NULL)
63 delete[] MAPEM.prj; //wei 3/2005
67 BOOLEAN emap_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
70 BOOLEAN flag; //added hstau 1/04
81 //------------------------ statements ----------------------------
85 MAPUser1 = fopen("MAPUser1", "w");
89 MAPEM.delta = MAX0((REAL) 1.0,
90 (REAL) 2.0 * (GeoPar.pixsiz / GeoPar.pinc));
94 MAPEM.delta = GeoPar.pixsiz / (GeoPar.pinc * (REAL) sqrt(2.) * (GeoPar.radius - (REAL) sqrt(.5) * GeoPar.pixsiz * GeoPar.nelem) / GeoPar.stod);
99 if (((INTEGER) MAPEM.delta) > GeoPar.usrays)
101 MAPEM.delta = (REAL) GeoPar.usrays;
104 // nelem not used. hstau
105 ReadInp(area, GeoPar.nelem, GeoPar.usrays, GeoPar.prjnum, &flag); //hstau 1/04
108 return TRUE; //added hstau 1/04
110 if (MAPEM.gamma1 <= Consts.zero)
115 if (GeoPar.nelem == 1)
116 MAPEM.gamma1 = 0.0; //wei 1/2005
118 if (!Blob.pix_basis && MAPEM.gamma1 > 0.0)
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
128 for (MAPEM.idy = 0; MAPEM.idy <= MAPEM.idymax - 1; MAPEM.idy++)
130 for (MAPEM.idx = 0; MAPEM.idx <= MAPEM.idxmax - 1; MAPEM.idx++)
133 return TRUE; //added by wei, 1/2005
139 EvalMAP(recon, iter);