src/snark/bsmooth.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/bsmooth.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  A NINE-POINT WEIGHTED, NON-LINEAR SMOOTHING PROCEDURE.
  11  FOR DETAILS, SEE THE SNARK MANUAL.
  12  */
  13
  14 #include <cstdio>
  15 #include <cmath>
  16
  17 #include "blkdta.h"
  18 #include "consts.h"
  19 #include "uiod.h"
  20 #include "blob.h"
  21
  22 void Blob_class::bsmooth(
  23 REAL* a,      // image
  24                 INTEGER M,    // number of rows
  25                 INTEGER H,    // number of columns
  26                 REAL thresh,  // treshold
  27                 REAL w1,      // weight of center
  28                 REAL w2      // weight of horizontal and diag neighbous
  29                 )
  30 {
  31         REAL w3 = w2;       // weight of diagonal neighbous
  32         REAL* row1 = new REAL[M]; // row buffer 1
  33         REAL* row2 = new REAL[M]; // row buffer 2
  34         REAL* pres = row2;        // current row before smooth
  35         REAL* prev = row1;        // previous row before smooth
  36         INTEGER lmf, lml;
  37         INTEGER lhf, lhl;
  38         INTEGER m, h;
  39         INTEGER m1, h1;
  40         INTEGER ind;
  41
  42         BOOLEAN once = TRUE; // message flag
  43
  44         REAL centre;      // curent grid point value
  45
  46         REAL smuth;       // weighted sum of neighbours
  47         REAL divide;      // sum of weights
  48         REAL value;       // neighbour value
  49
  50         lhf = (INTEGER) (-1 * (INTEGER) (H / 2));
  51         lhl = (INTEGER) ((H - 1) / 2);
  52
  53         lmf = -1 * (INTEGER) (M / 2);
  54         lml = (INTEGER) ((M - 1) / 2);
  55
  56         for (h = lhl; h >= lhf; h--)
  57         {
  58                 REAL* itemp = prev;  // swap the buffers (necessary?)
  59                 prev = pres;
  60                 pres = itemp;
  61
  62                 for (m = lmf; m <= lml; m++)
  63                 { // m and h must have the same parity!!
  64                         if ((m + h) % 2 == 0)
  65                         {
  66                                 m1 = m + Blob.M2;
  67                                 h1 = h + Blob.H2;
  68                                 ind = h1 * M + m1;
  69
  70                                 centre = a[ind];
  71                                 pres[m1] = centre;
  72                                 smuth = w1 * centre;
  73                                 divide = w1;
  74
  75                                 if (h != lhl)
  76                                 {                   // if not first row
  77                                         if (m != lmf)
  78                                         {                 // if not first column
  79                                                 value = prev[m1 - 1];       // prev row prev column  #1
  80
  81                                                 if ((REAL) fabs(value - centre) <= thresh)
  82                                                 {
  83                                                         smuth += w3 * value;
  84                                                         divide += w3;
  85                                                 }
  86                                         }
  87
  88                                         if (m != lml)
  89                                         {              // if not last column
  90                                                 value = prev[m1 + 1];      // prev row, next column   #2
  91
  92                                                 if ((REAL) fabs(value - centre) <= thresh)
  93                                                 {
  94                                                         smuth += w3 * value;
  95                                                         divide += w3;
  96                                                 }
  97                                         }
  98                                 }
  99
 100                                 if (m >= lmf + 2)
 101                                 {                   // if not first or second column
 102                                         value = pres[m1 - 2];        // same row, two prev column #3
 103
 104                                         if ((REAL) fabs(value - centre) <= thresh)
 105                                         {
 106                                                 smuth += w2 * value;
 107                                                 divide += w2;
 108                                         }
 109                                 }
 110
 111                                 if (m <= lml - 2)
 112                                 {        // if not last column or one before last one
 113                                         value = a[ind + 2];          // same row, two next column #4
 114
 115                                         if ((REAL) fabs(value - centre) <= thresh)
 116                                         {
 117                                                 smuth += w2 * value;
 118                                                 divide += w2;
 119                                         }
 120                                 }
 121
 122                                 if (h != lhf)
 123                                 {             // if not last row
 124                                         if (m != lmf)
 125                                         {                 // if not first column
 126                                                 value = a[ind - M - 1];     // next row prev column   #5
 127
 128                                                 if ((REAL) fabs(value - centre) <= thresh)
 129                                                 {
 130                                                         smuth += w3 * value;
 131                                                         divide += w3;
 132                                                 }
 133                                         }
 134
 135                                         if (m != lml)
 136                                         {           // if not last column
 137                                                 value = a[ind - M + 1];      // next row next column  #6
 138
 139                                                 if ((REAL) fabs(value - centre) <= thresh)
 140                                                 {
 141                                                         smuth += w3 * value;
 142                                                         divide += w3;
 143                                                 }
 144                                         }
 145                                 }
 146
 147                                 if ((REAL) fabs(divide) > Consts.zero)
 148                                 {
 149                                         a[ind] = smuth / divide;
 150                                 }
 151                                 else
 152                                 {
 153                                         if (once)
 154                                         {
 155                                                 fprintf(output, "\n *** warning - sum of weights=zero");
 156                                                 once = FALSE;
 157                                         }
 158
 159                                         a[ind] =
 160                                                         ((smuth - Consts.infin) > 0) ? (REAL) 1.0 :
 161                                                         (smuth == Consts.infin) ? (REAL) 0.0 : (REAL) -1.0;
 162                                 }
 163
 164                                 //l+=2;    //blob grid
 165                         }
 166                 }
 167         }
 168
 169         delete[] row1;  // bug 92 - Lajos - 03/02/2005
 170         delete[] row2;  // bug 92 - Lajos - 03/02/2005
 171
 172         return;
 173 }