src/snark/emap_CMAP.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_CMAP.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  ======================================================================
  11  one iteration of the maximum a-posteriori reconstruction algorithm
  12  for emission tomography
  13  ======================================================================
  14  */
  15
  16 #include <cstdio>
  17 #include <cmath>
  18
  19 #include "blkdta.h"
  20 #include "geom.h"
  21 #include "consts.h"
  22 #include "blob.h"
  23
  24 #include "emap.h"
  25
  26 BOOLEAN emap_class::MAP(REAL* recon)
  27 {
  28         //------ local variables
  29         INTEGER i, j, first, np, nr, idelta, usnr;
  30         REAL p, q;
  31         REAL yi;
  32         INTEGER area;
  33
  34
  35         if (Blob.pix_basis)
  36         {
  37                 area = GeoPar.area;
  38         }
  39
  40         else
  41         {
  42                 area = Blob.area;
  43         }
  44
  45         idelta = (INTEGER) (MAPEM.delta);
  46
  47         //====== step 1: compute q for iteration k+1:
  48         for (j = 0; j < area; j++)
  49         {
  50                 MAPEM.qj[j] = 0.0;
  51         }
  52
  53         for (np = 0; np < GeoPar.prjnum; np++)
  54         {
  55                 for (first = 0; first < idelta; first++)
  56                 {
  57                         for (usnr = first; usnr < GeoPar.usrays; usnr += idelta)
  58                         {
  59                                 nr = usnr + GeoPar.fusray;
  60                                 yi = MAPEM.prj[(np) * GeoPar.usrays + usnr];
  61                                 if (yi != 0.0)
  62                                 {
  63                                         if (findqpart(np, nr, yi, recon, MAPEM.qj))
  64                                                 return TRUE; //if findqpart return false, EMAP is not applicable, bywei 1/2005
  65                                 }
  66                         }
  67                 }
  68         }
  69
  70         // MAPEM.idy = MAPEM.idx = 0, see comments in emap.cpp
  71
  72 #ifdef HSTAU
  73         fprintf(output,"\n #IDY=%d #IDX=%d\n", MAPEM.idy, MAPEM.idx);fflush(output);
  74 #endif
  75
  76         if (Blob.pix_basis)
  77         {
  78                 for (i = MAPEM.idy * GeoPar.nelem;
  79                                 i < GeoPar.nelem * (GeoPar.nelem - 1) + 1;
  80                                 i += GeoPar.nelem * MAPEM.idymax)
  81                 {
  82                         for (j = i + MAPEM.idx; j < i + GeoPar.nelem; j += MAPEM.idxmax)
  83                         {
  84
  85                                 if (MAPEM.gamma1 > 0.0)
  86                                 {
  87                                         MAPEM.qj[j] *= recon[j] * MAPEM.Sjjinv[j] / MAPEM.gamma1;
  88                                 }
  89                                 else
  90                                 {
  91                                         if (MAPEM.Wjj[j] > Consts.zero)
  92                                                 MAPEM.qj[j] *= recon[j] / MAPEM.Wjj[j];
  93                                         else
  94                                                 MAPEM.qj[j] = recon[j]; //keep recon[j] unchanged when qj[j]is zero, by wei 1/2005
  95                                 }
  96                         }
  97                 }
  98         }
  99         else
 100         { //blob case
 101                 for (i = MAPEM.idy * Blob.M; i < Blob.M * (Blob.H - 1) + 1;
 102                                 i += Blob.M * MAPEM.idymax)
 103                 {
 104                         for (j = i + MAPEM.idx; j < i + Blob.M; j += MAPEM.idxmax)
 105                         {
 106
 107                                 if (j % 2 == Blob.pr)
 108                                 {
 109                                         if (MAPEM.Wjj[j] > Consts.zero)
 110                                                 MAPEM.qj[j] *= recon[j] / MAPEM.Wjj[j];
 111                                         else
 112                                                 MAPEM.qj[j] = recon[j]; //recon[j] keep unchanged when Wjj[j] is zero, by wei 1/2005
 113                                 }
 114                         }
 115                 }
 116         }
 117
 118         if (MAPEM.gamma1 > 0.0)
 119         {  // blob case will not enter here.hstau
 120                 //====== step 2: compute "smoothing component" from "recon" to "Sx"
 121                 applyK(recon, MAPEM.Kx, GeoPar.nelem);
 122
 123 //     ---Set edges to 0.0
 124                 for (i = 0; i < GeoPar.nelem; i++)
 125                 {
 126                         MAPEM.Kx[i] = 0.0;          // first row
 127                         //MAPEM.Kx[(i-1)*GeoPar.nelem+1] = 0.0;
 128                         MAPEM.Kx[(i + 1) * GeoPar.nelem - 1] = 0.0; // last column
 129                         MAPEM.Kx[i * GeoPar.nelem] = 0.0;     // first column
 130                         MAPEM.Kx[(GeoPar.nelem - 1) * GeoPar.nelem + i] = 0.0; // last row
 131                 }
 132
 133                 applyK(MAPEM.Kx, MAPEM.Sx, GeoPar.nelem);
 134
 135                 //====== step 3: compute x for iteration k+1 from "Sx" and "qj"
 136                 //cvd$     nosync
 137                 for (i = MAPEM.idy * GeoPar.nelem;
 138                                 i < GeoPar.nelem * (GeoPar.nelem - 1) + 1;
 139                                 i += GeoPar.nelem * MAPEM.idymax)
 140                 {
 141                         for (j = i + MAPEM.idx; j < i + GeoPar.nelem; j += MAPEM.idxmax)
 142                         {
 143                                 p = MAPEM.Wjj[j] / (REAL) 9.0 - recon[j]
 144                                                 + MAPEM.Sx[j] * MAPEM.Sjjinv[j] / (REAL) 9.0;
 145                                 q = MAPEM.qj[j] / (REAL) 9.0;
 146                                 if ((p <= 0.0) || (q / (p * p) >= 0.005))
 147                                 {
 148                                         MAPEM.qj[j] = (REAL) 0.5
 149                                                         * (-p + (REAL) sqrt(p * p + 4 * q));
 150                                 }
 151                                 else
 152                                 {
 153                                         MAPEM.qj[j] = q / p * (1 - (REAL) 0.5 * (q / (p * p)));
 154                                 }
 155                         }
 156                 }
 157         }    // not with the blob case .hstau
 158
 159         //====== step 4: replace k-th iteration by (k+1)-th iteration in the
 160         //       reconstruction
 161         if (Blob.pix_basis)
 162         {
 163                 for (i = MAPEM.idy * GeoPar.nelem;
 164                                 i < GeoPar.nelem * (GeoPar.nelem - 1) + 1;
 165                                 i += GeoPar.nelem * MAPEM.idymax)
 166                 {
 167                         for (j = i + MAPEM.idx; j < i + GeoPar.nelem; j += MAPEM.idxmax)
 168                         {
 169                                 if (MAPEM.qj[j] < Consts.zero)
 170                                 {
 171                                         recon[j] = Consts.zero;
 172                                 }
 173                                 else
 174                                 {
 175                                         recon[j] = MAPEM.qj[j];
 176                                 }
 177                         }
 178                 }
 179         }
 180
 181         else
 182         {  //if blobs
 183                 for (i = MAPEM.idy * Blob.M; i < Blob.M * (Blob.H - 1) + 1;
 184                                 i += Blob.M * MAPEM.idymax)
 185                 {
 186                         for (j = i + MAPEM.idx; j < i + Blob.M; j += MAPEM.idxmax)
 187                         {
 188
 189                                 if (j % 2 == Blob.pr)
 190                                 { //only the hex grid points
 191                                         if (MAPEM.qj[j] < Consts.zero)
 192                                         {
 193                                                 recon[j] = Consts.zero;
 194                                         }
 195                                         else
 196                                         {
 197                                                 recon[j] = MAPEM.qj[j];
 198                                         }
 199                                 }
 200                         }
 201                 }
 202         }
 203
 204         return FALSE;
 205 }
 206
 207 /*
 208  c======================================================================
 209  c     source & dest are nelem*nelem arrays.  source contains a picture
 210  c     to be convolved with the 3*3 kernel whose central element is 1
 211  c     and the rest are -1/8.  The result goes in dest.
 212  c======================================================================
 213  */
 214
 215 void emap_class::applyK(REAL* source, REAL* dest, INTEGER nelem)
 216 {
 217         INTEGER i, j;
 218
 219         for (i = 0; i < nelem * nelem; i++)
 220         {
 221                 dest[i] = 0.0;
 222         }
 223
 224         for (i = 0; i < nelem - 1; i++)
 225         {
 226                 for (j = 0; j < nelem; j++)
 227                 {
 228                         dest[i + (j) * nelem] += source[i + (j) * nelem + 1];
 229
 230                 }
 231         }
 232
 233         for (i = 0; i < nelem - 1; i++)
 234         {
 235                 for (j = 0; j < nelem; j++)
 236                 {
 237                         dest[i + (j) * nelem + 1] += source[i + (j) * nelem];
 238
 239                 }
 240         }
 241
 242         for (i = 0; i < nelem; i++)
 243         {
 244                 for (j = 0; j < nelem - 1; j++)
 245                 {
 246                         dest[i + (j) * nelem] += source[i + (j) * nelem + nelem];
 247                 }
 248         }
 249
 250         for (i = 0; i < nelem; i++)
 251         {
 252                 for (j = 0; j < nelem - 1; j++)
 253                 {
 254
 255                         dest[i + (j) * nelem + nelem] += source[i + (j) * nelem];
 256
 257                 }
 258         }
 259
 260         for (i = 0; i < nelem - 1; i++)
 261         {
 262                 for (j = 0; j < nelem - 1; j++)
 263                 {
 264                         dest[i + (j) * nelem] += source[i + (j) * nelem + nelem + 1];
 265                 }
 266         }
 267
 268         for (i = 0; i < nelem - 1; i++)
 269         {
 270                 for (j = 0; j < nelem - 1; j++)
 271                 {
 272                         dest[i + (j) * nelem + nelem + 1] += source[i + (j) * nelem];
 273
 274                 }
 275         }
 276
 277         for (i = 0; i < nelem - 1; i++)
 278         {
 279                 for (j = 0; j < nelem - 1; j++)
 280                 {
 281                         dest[i + (j) * nelem + nelem] += source[i + (j) * nelem + 1];
 282
 283                 }
 284         }
 285
 286         for (i = 0; i < nelem - 1; i++)
 287         {
 288                 for (j = 0; j < nelem - 1; j++)
 289                 {
 290                         dest[i + (j) * nelem + 1] += source[i + (j) * nelem + nelem];
 291
 292                 }
 293         }
 294
 295         sscal(nelem * nelem, -1 / 8., dest, 1);
 296
 297         for (i = 0; i < nelem * nelem; i++)
 298         {
 299                 dest[i] += source[i];
 300
 301         }
 302
 303         return;
 304 }