src/snark/emap_CInitMAP.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_CInitMAP.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  ===========================================================================
  11  initialize working space for the MAP algorithm
  12  ===========================================================================
  13  */
  14
  15 #include <cstdio>
  16
  17 #include "blkdta.h"
  18 #include "geom.h"
  19 #include "uiod.h"
  20 #include "anglst.h"
  21 #include "wray.h"
  22 #include "blob.h"
  23
  24 #include "emap.h"
  25
  26 REAL getSjj(INTEGER j, INTEGER nelem);
  27
  28 void emap_class::InitMAP(REAL* recon)
  29 {
  30
  31         INTEGER testneg;
  32         //------ local variables
  33         INTEGER i, first, np, nr, idelta, numb, list[4096], usnr;
  34         REAL snorm, weight[4096];
  35
  36         INTEGER area;
  37
  38         if (Blob.pix_basis)
  39         {
  40                 area = GeoPar.area;
  41         }
  42         else
  43         {
  44                 area = Blob.area;
  45         }
  46
  47         idelta = (INTEGER) (MAPEM.delta);
  48
  49         //------ set image to average density and initialize the work area for
  50         //       the total length of all projections through each pixel
  51         //cvd$  nodepchk
  52
  53         for (i = 0; i < area; i++)
  54         {
  55                 MAPEM.Wjj[i] = 0.0;
  56                 if (Blob.pix_basis)
  57                         MAPEM.Sjjinv[i] = 1 / getSjj(i, GeoPar.nelem); //blob case will not use Sjj Hstau 6/25/2003
  58         }
  59
  60         if (Blob.pix_basis)
  61         { //not needed in the blob case Hstau 6/3/2003
  62                 for (i = 0; i < area; i++)
  63                 {
  64                         MAPEM.Kx[i] = 0.0;
  65                         MAPEM.Sx[i] = 0.0;
  66                 }
  67         }
  68
  69         // compute the total length of all projections through each pixel
  70
  71         for (np = 0; np < GeoPar.prjnum; np++)
  72         {
  73                 for (first = 0; first < idelta; first++)
  74                 {
  75                         for (nr = first; nr < GeoPar.nrays; nr += idelta)
  76                         {
  77                                 if (Blob.pix_basis)
  78                                 {
  79                                         wray(np, nr, list, weight, &numb, &snorm);
  80                                 }
  81                                 else
  82                                 {
  83                                         Blob.bwray(np, nr, list, weight, &numb, &snorm);
  84                                 }
  85
  86                                 for (i = 0; i < numb; i++)
  87                                 {
  88
  89                                         if (MAPEM.gamma1 > 0.0)
  90                                         {  //when blobs, gamma1 = 0.0
  91                                                 MAPEM.Wjj[list[i]] += weight[i] / MAPEM.gamma1
  92                                                                 * MAPEM.Sjjinv[list[i]];
  93                                         }
  94                                         else
  95                                         {
  96                                                 MAPEM.Wjj[list[i]] += weight[i];
  97                                         }
  98                                 }
  99                         }
 100                 }
 101         }
 102
 103         //------ get projection data
 104
 105         testneg = 0;
 106         for (np = 0; np < GeoPar.prjnum; np++)
 107         {
 108                 for (usnr = 0; usnr < GeoPar.usrays; usnr++)
 109                 {
 110
 111                         //nr = usnr + GeoPar.fusray - 1;
 112                         nr = usnr + GeoPar.fusray; // mk
 113
 114                         // T.K. Narayan's modification -  24.Apr.1996
 115                         // setting negative projection data to 0.0
 116
 117                         if (Anglst.prdta(np, nr) >= 0.0)
 118                         {
 119                                 //MAPEM.prj[(np-1)*GeoPar.usrays+usnr] = prdta(np,nr);
 120                                 MAPEM.prj[(np) * GeoPar.usrays + usnr] = Anglst.prdta(np, nr);
 121                         }
 122                         else
 123                         {
 124                                 //MAPEM.prj[(np-1)*GeoPar.usrays+usnr] = 0.0;
 125                                 MAPEM.prj[(np) * GeoPar.usrays + usnr] = 0.0;
 126                                 testneg++;
 127                                 if (testneg == 1)
 128                                 {
 129                                         fprintf(output,
 130                                                         "\n     *** WARNING - There exist negative values in the projection data ***"); //bug 127, hstau 7/8/2005
 131                                         fprintf(output,
 132                                                         "\n     *** Negative values set to 0.0 ***");
 133                                 }
 134                         }
 135                 }
 136         }
 137
 138         return;
 139 }
 140
 141 REAL getSjj(INTEGER j, INTEGER nelem)
 142 // integer j, nelem
 143 {
 144         INTEGER jx, jy;
 145
 146         j++; // MK
 147
 148         if (nelem == 3)
 149         {
 150                 if (j == 5)
 151                         return 1.;
 152                 else
 153                         return 1 / 64.;
 154         }                                //add exception when nelem=3,by wei, 1/2005
 155
 156         jx = j - ((j - 1) / nelem) * nelem;
 157         jy = (j - 1) / nelem + 1;
 158         if ((jx > 2) && (jx < nelem - 1))
 159         {
 160                 if ((jy > 2) && (jy < nelem - 1))
 161                 {
 162                         return 9 / 8.;
 163                 }
 164                 else
 165                 {
 166                         if ((jy == 2) || (jy == nelem - 1))
 167                         {
 168                                 return 69 / 64.;
 169                         }
 170                         else
 171                         {
 172                                 return 3 / 64.;
 173                         }
 174                 }
 175         }
 176         else
 177         {
 178                 if ((jx == 2) || (jx == nelem - 1))
 179                 {
 180                         if ((jy > 2) && (jy < nelem - 1))
 181                         {
 182                                 return 69 / 64.;
 183                         }
 184                         else
 185                         {
 186                                 if ((jy == 2) || (jy == nelem - 1))
 187                                 {
 188                                         return 67 / 64.;
 189                                 }
 190                                 else
 191                                 {
 192                                         return 1 / 32.;
 193                                 }
 194                         }
 195                 }
 196                 else
 197                 {
 198                         if ((jy > 2) && (jy < nelem - 1))
 199                         {
 200                                 return 3 / 64.;
 201                         }
 202                         else
 203                         {
 204                                 if ((jy == 2) || (jy == nelem - 1))
 205                                 {
 206                                         return 1 / 32.;
 207                                 }
 208                                 else
 209                                 {
 210                                         return 1 / 64.;
 211                                 }
 212                         }
 213                 }
 214         }
 215 }