src/snark/blob_bpix.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_bpix.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 a task similar to bwray() but pixels instead of the rays
  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::bpix(INTEGER nv, INTEGER nh, INTEGER *list, REAL* weight,
  24                 INTEGER *numb, REAL *snorm)
  25 {
  26         INTEGER multp;
  27         INTEGER m, h, m1, h1;
  28         INTEGER mc, hc;
  29         INTEGER mf, hf, ml, hl;
  30         INTEGER M2, H2;
  31         INTEGER nelem2;
  32         INTEGER mx, hy;
  33         INTEGER ind, n;
  34         INTEGER m0, h0;
  35
  36         REAL xc, yc;
  37         REAL dinv2 = (REAL) 2.0 / Blob.delta;
  38         REAL dinv2sq3 = dinv2 * Consts.isqrt3;
  39
  40         REAL dist, distx, disty; // dist is the distance betw a blob center and a pixel center. (distx, disty) is the vector distance
  41         REAL dispx, dispy; // (dispx, dispy) is the displacement from the center of a blob to the center of the nearest pixel.
  42
  43         multp = (INTEGER) floor(Blob.support / Blob.delta + 0.5) + 1;
  44
  45         H2 = (INTEGER) (Blob.H / 2);
  46         M2 = (INTEGER) (Blob.M / 2);
  47         nelem2 = (INTEGER) (GeoPar.nelem / 2);
  48
  49         REAL* footpx = new REAL[4 * multp + 1];
  50         REAL* footpy = new REAL[4 * multp + 1];
  51
  52         for (m0 = -2 * multp; m0 <= 2 * multp; m0++)
  53         {
  54                 footpx[m0 + 2 * multp] = (REAL) m0 * Blob.di2;
  55         }
  56
  57         for (h0 = -2 * multp; h0 <= 2 * multp; h0++)
  58         {
  59                 footpy[h0 + 2 * multp] = (REAL) h0 * Blob.di2sq3;
  60         }
  61
  62 // coord of the center of the pixel (nh, nv)
  63         xc = GeoPar.pixsiz * (nh - nelem2);
  64         yc = GeoPar.pixsiz * (nelem2 - nv);
  65
  66         mx = (INTEGER) floor(dinv2 * xc + 0.5); // coord of the center of
  67
  68         hy = (INTEGER) floor(dinv2sq3 * yc + 0.5);        // the nearest blob
  69
  70         mc = mx;
  71         hc = hy;
  72
  73         // displacement of the pixel footprint with resp to the nearest blob center
  74         dispx = di2 * mx - xc;
  75         dispy = di2sq3 * hy - yc;
  76
  77         mf = MAX0(lmf, mc - 2 * multp);
  78         ml = MIN0(lml, mc + 2 * multp);
  79
  80         hf = MAX0(lhf, hc - 2 * multp);
  81         hl = MIN0(lhl, hc + 2 * multp);
  82
  83         (*numb) = 0;
  84         (*snorm) = 0.0;
  85
  86         if (ml < lmf || mf > lml || hl < lhf || hf > lhl)
  87                 return;
  88         for (h = hf; h <= hl; h++)
  89         {
  90                 h0 = h - hc;
  91                 disty = footpy[h0 + 2 * multp] + dispy; //the y in the x-y coord and the h in the m-h coord are in the same direction
  92                 if (fabs(disty) > Blob.support)
  93                         continue;
  94
  95                 for (m = mf; m <= ml; m++)
  96                 {
  97                         m0 = m - mc;
  98                         distx = footpx[m0 + 2 * multp] + dispx;
  99                         if (fabs(distx) > Blob.support)
 100                                 continue;
 101
 102                         if ((m + h) % 2 == 0)
 103                         {
 104                                 m1 = m + M2;
 105                                 h1 = h + H2;
 106                                 ind = h1 * Blob.M + m1;
 107
 108                                 dist = sqrt(distx * distx + disty * disty);
 109                                 if (dist < Blob.support)
 110                                 {
 111                                         n = (INTEGER) round(dist * sampiblob); // changed "rint" to "round". Lajos, Jan 25, 2005
 112                                         list[(*numb)] = ind;
 113                                         weight[(*numb)] = Blob.vtable[n];
 114
 115                                         (*snorm) += (Blob.vtable[n] * Blob.vtable[n]);
 116                                         (*numb)++;
 117                                 }
 118                         } //for m+h
 119                 } // for m
 120         }  //for h
 121
 122         delete[] footpx;  // bug 92 - Lajos - 03/02/2005
 123         delete[] footpy;  // bug 92 - Lajos - 03/02/2005
 124         return;
 125 }
 126