src/snark/blob_footp.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/blob_footp.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  subroutine to perform the forward (1) or backward (-1) "projection" using blob footprints
  11
  12  */
  13
  14 #include <cstdio>
  15 #include <cmath>
  16
  17 #include "blkdta.h"
  18 #include "blob.h"
  19 #include "geom.h"
  20 #include "consts.h"
  21 #include "uiod.h"
  22
  23 void Blob_class::footp(REAL *blob_array, REAL *pict, INTEGER option)
  24 {
  25         INTEGER mult; //max number of pixels whose center is covered by a blob
  26         INTEGER m, h, m1, h1; // (m,h) coord of the center of a blob in the underlying sq grid. m1 * Blob.H + h1 index in the blob array.
  27         INTEGER k, j; //(j,k) coord of the center of a pixel
  28         INTEGER k0, j0;
  29         INTEGER jc, kc; // index of the nearest pixel
  30         INTEGER jf, kf, jl, kl; //upper and lower bound for j and k given a blob
  31         INTEGER nelem2;
  32         INTEGER jx, ky;
  33         INTEGER ind, n;
  34
  35         REAL xc, yc, auy;
  36         REAL ipixsiz = (REAL) 1.0 / GeoPar.pixsiz;
  37         REAL dist, distx, disty; // dist is the distance betw a blob center and a pixel center. (distx, disty) is the vector distance
  38         REAL dispx, dispy; // (dispx, dispy) is the displacement from the center of a blob to the center of the nearest pixel.
  39
  40         mult = (INTEGER) floor(Blob.support / GeoPar.pixsiz + 0.5) + 1;
  41
  42         nelem2 = (INTEGER) (GeoPar.nelem / 2);
  43
  44         REAL* footx = new REAL[2 * mult + 1];
  45         REAL* footy = new REAL[2 * mult + 1];
  46
  47         for (k0 = -mult; k0 <= mult; k0++)
  48         {
  49                 footy[k0 + mult] = (REAL) k0 * GeoPar.pixsiz;
  50         }
  51
  52         for (j0 = -mult; j0 <= mult; j0++)
  53         {
  54                 footx[j0 + mult] = (REAL) j0 * GeoPar.pixsiz;
  55         }
  56
  57         // initialization for forward operation
  58         if (option == 1)
  59         {
  60                 for (k = 0; k < GeoPar.nelem; k++)
  61                 {
  62                         for (j = 0; j < GeoPar.nelem; j++)
  63                         {
  64                                 pict[k * GeoPar.nelem + j] = 0.0;
  65                         }
  66                 }
  67         }
  68
  69         //initialization for backward operation
  70         else if (option == -1)
  71         {
  72                 for (j = 0; j <= Blob.area; j++)
  73                 {
  74                         if (j % 2 == Blob.pr)
  75                         {
  76                                 blob_array[j] = 0.0;
  77                         }
  78                 }
  79         }
  80
  81         // forward or backward operation
  82
  83         for (h = lhf; h <= lhl; h++)
  84         {
  85                 for (m = lmf; m <= lml; m++)
  86                 { // m and h must have the same parity!!
  87                         if ((m + h) % 2 == 0)
  88                         {
  89
  90                                 xc = di2 * m;  // coord of the center of the blob (m,h)
  91                                 yc = di2sq3 * h;
  92                                 auy = ipixsiz * yc;
  93                                 jx = (INTEGER) floor(ipixsiz * xc + 0.5); // coord of the center of
  94
  95                                 ky = (INTEGER) floor(auy + 0.5);        // the nearest pixel
  96
  97                                 jc = jx + nelem2; // bring to the right address
  98                                 kc = (INTEGER) floor(-auy + 0.5) + nelem2;
  99
 100                                 // displacement of the footprint with resp to the nearest pixel center
 101                                 dispx = jx * GeoPar.pixsiz - xc;
 102
 103                                 dispy = ky * GeoPar.pixsiz - yc;
 104
 105                                 jf = MAX0(ljf, jc - mult);
 106                                 jl = MIN0(ljl, jc + mult);
 107
 108                                 kf = MAX0(lkf, kc - mult);
 109                                 kl = MIN0(lkl, kc + mult);
 110
 111                                 if (jl < ljf || jf > ljl || kl < lkf || kf > lkl)
 112                                         continue;
 113
 114                                 m1 = m + Blob.M2;
 115                                 h1 = h + Blob.H2;
 116
 117                                 ind = h1 * Blob.M + m1;
 118
 119                                 for (k = kf; k <= kl; k++)
 120                                 {
 121                                         k0 = k - kc;
 122                                         disty = footy[k0 + mult] - dispy; //the y in the x-y coord and the k in the j-k coord are in the opposite directions
 123                                         if (fabs(disty) > Blob.support)
 124                                                 continue;
 125
 126                                         for (j = jf; j <= jl; j++)
 127                                         {
 128                                                 j0 = j - jc;
 129                                                 distx = footx[j0 + mult] + dispx;
 130                                                 if (fabs(distx) > Blob.support)
 131                                                         continue;
 132
 133                                                 dist = sqrt(distx * distx + disty * disty);
 134                                                 if (dist < Blob.support)
 135                                                 {
 136                                                         n = (INTEGER) round(dist * sampiblob); // changed "rint" to "round". Lajos, Jan 25, 2005
 137                                                         //forward case
 138                                                         if (option == 1)
 139                                                         {
 140                                                                 pict[k * GeoPar.nelem + j] += (blob_array[ind]
 141                                                                                 * Blob.vtable[n]);
 142                                                         }
 143
 144                                                         //backward case
 145                                                         else if (option == -1)
 146                                                         {
 147                                                                 blob_array[ind] += (pict[k * GeoPar.nelem + j]
 148                                                                                 * Blob.vtable[n]);
 149                                                         }
 150                                                 }
 151
 152                                         }
 153                                 }
 154
 155                         }
 156                 }
 157         }
 158
 159         delete[] footx;  // bug 92 - Lajos - 03/02/2005
 160         delete[] footy;  // bug 92 - Lajos - 03/02/2005
 161         return;
 162 }
 163