src/snark/pix2blob.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/pix2blob.cpp $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9  */
  10
  11 /* This subroutine converts an image function in pixel basis into that in blob basis. This subroutine uses SIRT type algorithm with blob footprints for speeding up. See subroutine blob2pix() for the meaning of various variables.
  12  */
  13
  14 #include <cstdio>
  15 #include <cmath>
  16 #include <cstdlib>
  17
  18 #include "blkdta.h"
  19 #include "blob.h"
  20 #include "geom.h"
  21 #include "consts.h"
  22 #include "effpick.h"
  23 #include "uiod.h"
  24
  25 void Blob_class::pix2blob(REAL* pict, REAL* blob_array)
  26 {
  27         INTEGER m, h, m1, h1;
  28         INTEGER i, k, j, ii, kk, jj;
  29         INTEGER jc, kc;
  30         INTEGER M2, H2;
  31         INTEGER nelem2;
  32         INTEGER jx, ky;
  33         INTEGER ind, indb;
  34
  35
  36         REAL xc, yc, auy;
  37         REAL ipixsiz = (REAL) 1.0 / GeoPar.pixsiz;
  38         INTEGER count, ncount;
  39         REAL maxpict, minpict;
  40         REAL amppict;
  41         REAL nc;
  42         REAL sd;
  43         INTEGER niter;
  44         REAL width;
  45         REAL rneigh;
  46         REAL tau;
  47         REAL unif;
  48         REAL frelax;
  49
  50         frelax = 2.0; //changed on 4/8/04, from 1.0 to 2.0, as suggested by Laszlo
  51         tau = 0.3;
  52         unif = 0.0005;
  53         width = 3.5;
  54         rneigh = 0.75;
  55         niter = 5000;
  56
  57         frelax *= GeoPar.pixsiz * GeoPar.pixsiz
  58                         / (Consts.pi * Blob.support * Blob.support);
  59         frelax = MIN0(2.0, frelax);  //frelax param cannot be greater than 2 !
  60
  61         REAL *u = new REAL[GeoPar.area]; //art4 dual var
  62         REAL *Delta = new REAL[GeoPar.area]; //art4 epsilon
  63         REAL *Gamma = new REAL[GeoPar.area]; //art4 epsilon
  64         INTEGER *que = new INTEGER[GeoPar.area]; //art4
  65
  66         INTEGER pow2, num_dig;
  67         INTEGER *decomp = NULL;
  68         INTEGER *base = NULL;      //wei 3/2005
  69         INTEGER *old_numk, *new_numk;
  70         INTEGER *old_numj, *new_numj;
  71
  72         H2 = (INTEGER) (Blob.H / 2);
  73         M2 = (INTEGER) (Blob.M / 2);
  74         nelem2 = (INTEGER) (GeoPar.nelem / 2);
  75
  76         //get the decomposition of the smallest power of 2 >= nelem
  77         pow2 = 1;
  78
  79         while (pow2 < GeoPar.nelem)
  80         {
  81                 pow2 *= 2;
  82         }
  83
  84         factorization(pow2, &decomp, &num_dig);
  85         get_base(decomp, num_dig, &base);
  86         old_numk = new INTEGER[num_dig];
  87         new_numk = new INTEGER[num_dig];
  88         old_numj = new INTEGER[num_dig];
  89         new_numj = new INTEGER[num_dig];
  90
  91         // initialising the blobs with the value of the nearest pixel and computing mod2 using BLOB footprint
  92
  93         for (h = lhf; h <= lhl; h++)
  94         {
  95                 for (m = lmf; m <= lml; m++)
  96                 { // m and h must have the same parity!!
  97                         if ((m + h) % 2 == 0)
  98                         {
  99                                 m1 = m + M2;
 100                                 h1 = h + H2;
 101
 102                                 indb = h1 * Blob.M + m1;
 103
 104                                 xc = di2 * m;  // coord of the center of the blob (m,h)
 105                                 yc = di2sq3 * h;
 106
 107                                 auy = ipixsiz * yc;
 108                                 jx = (INTEGER) floor(ipixsiz * xc + 0.5); // coord of the center of
 109                                 ky = (INTEGER) floor(auy + 0.5);        // the nearest pixel
 110
 111                                 jc = jx + nelem2; // bring to the right address
 112                                 kc = (INTEGER) floor(-auy + 0.5) + nelem2;
 113
 114                                 if (jc > ljl)
 115                                         jc = ljl;
 116                                 if (jc < ljf)
 117                                         jc = ljf;
 118
 119                                 if (kc > lkl)
 120                                         kc = lkl;
 121                                 if (kc < lkf)
 122                                         kc = lkf;
 123
 124                                 blob_array[indb] = pict[kc * GeoPar.nelem + jc];
 125                         }
 126                 }
 127         }
 128
 129         //initialize dual var and finding the max amplitud between pict and 0
 130         maxpict = -Consts.infin;
 131         minpict = Consts.infin;
 132         for (k = 0; k < GeoPar.nelem; k++)
 133         {
 134                 for (j = 0; j < GeoPar.nelem; j++)
 135                 {
 136                         ind = k * GeoPar.nelem + j;
 137                         if (maxpict < pict[ind])
 138                         {
 139                                 maxpict = pict[ind];
 140                         }
 141                         if (minpict > pict[ind])
 142                         {
 143                                 minpict = pict[ind];
 144                         }
 145                         u[ind] = 0.0;
 146                 }
 147         }
 148         amppict = MAX0(maxpict,0.0) - MIN0(minpict, 0.0);
 149
 150         //computing delta and gamma of ART4
 151         for (k = 0; k < GeoPar.nelem; k++)
 152         {
 153                 for (j = 0; j < GeoPar.nelem; j++)
 154                 {
 155                         ind = k * GeoPar.nelem + j;
 156                         p2b_neighb(k, j, pict, &nc, rneigh);
 157                         if (pict[ind] >= nc)
 158                         {
 159                                 Delta[ind] = unif * amppict;
 160                                 Gamma[ind] = -(tau * fabs(pict[ind] - nc) + unif * amppict);
 161                         }
 162                         else
 163                         {
 164                                 Delta[ind] = tau * fabs(pict[ind] - nc) + unif * amppict;
 165                                 Gamma[ind] = -unif * amppict;
 166                         }
 167                 }
 168         }
 169
 170         i = 0;
 171         for (;;)
 172         {
 173                 //initializing list "que" and count
 174                 count = 0;
 175                 old_numk[0] = -1;
 176                 for (kk = 0; kk < pow2; kk++)
 177                 {
 178                         pi_map(num_dig, decomp, base, &old_numk, &new_numk, &k);
 179                         if (k < GeoPar.nelem)
 180                         {
 181                                 old_numj[0] = -1;
 182                                 for (jj = 0; jj < pow2; jj++)
 183                                 {
 184                                         pi_map(num_dig, decomp, base, &old_numj, &new_numj, &j);
 185                                         if (j < GeoPar.nelem)
 186                                         {
 187                                                 ind = k * GeoPar.nelem + j;
 188                                                 que[count] = ind;
 189                                                 count++;
 190                                         }
 191                                 } //for jj
 192                         }
 193                 } //for kk
 194                   //ART4-type algorithm
 195                 ii = 0;
 196                 for (;;)
 197                 {
 198                         Correct(pict, blob_array, u, Delta, Gamma, que, count, &ncount, &sd,
 199                                         width, frelax);
 200                         count = ncount;
 201
 202                         if (count == 0 || i >= niter)
 203                                 break;
 204                         i++;
 205                         ii++;
 206                 }
 207                 if (ii == 0 || i >= niter)
 208                         break;
 209         }
 210
 211         delete[] old_numk;  // bug 92 - Lajos - 03/02/2005
 212         delete[] new_numk;  // bug 92 - Lajos - 03/02/2005
 213         delete[] old_numj;  // bug 92 - Lajos - 03/02/2005
 214         delete[] new_numj;  // bug 92 - Lajos - 03/02/2005
 215         delete[] u;  // bug 92 - Lajos - 03/02/2005
 216         delete[] Delta;  // bug 92 - Lajos - 03/02/2005
 217         delete[] Gamma;  // bug 92 - Lajos - 03/02/2005
 218         delete[] que;  // bug 92 - Lajos - 03/02/2005
 219         if (decomp != NULL)
 220                 delete[] decomp;
 221         if (base != NULL)
 222                 delete[] base;      //wei 3/2005
 223
 224         return;
 225 }
 226
 227 void Blob_class::Correct(REAL* pict, REAL* blob_array, REAL* u, REAL* Delta,
 228                 REAL* Gamma, INTEGER* que, INTEGER count, INTEGER* ncount, REAL* sd,
 229                 REAL width, REAL frelax)
 230 {
 231
 232         INTEGER ind, indb;
 233         INTEGER c;
 234         INTEGER k, j;
 235         INTEGER numb;
 236         INTEGER multp = (INTEGER) floor(Blob.support / Blob.delta + 0.5) + 1;
 237         INTEGER npts = 4 * multp * multp;
 238         INTEGER *list = new INTEGER[npts];
 239         INTEGER g;
 240         REAL *weight = new REAL[npts];
 241         REAL sum;
 242         REAL snorm;
 243         REAL diff;
 244         REAL lowb, uppb;
 245         REAL ck;
 246
 247         *sd = 0;
 248         *ncount = 0;
 249
 250         for (c = 0; c < count; c++)
 251         {
 252                 ind = que[c];
 253                 sum = 0.0;
 254                 k = ind / GeoPar.nelem;
 255                 j = ind % GeoPar.nelem;
 256                 bpix(k, j, list, weight, &numb, &snorm);
 257                 if (numb > 0)
 258                 {
 259                         for (g = 0; g < numb; g++)
 260                         {
 261                                 indb = list[g];
 262                                 sum += blob_array[indb] * weight[g];
 263                         }
 264                         diff = pict[ind] - sum;
 265                         (*sd) += (diff * diff);
 266                         if (!(diff <= -width * Gamma[ind])
 267                                         || !(diff >= -width * Delta[ind])) //Ran bug 235 August 17, 2008
 268                         {
 269                                 que[(*ncount)] = ind;
 270                                 (*ncount)++;
 271                                 //compute the change
 272                                 lowb = frelax * (diff + Gamma[ind]) / snorm;
 273                                 uppb = frelax * (diff + Delta[ind]) / snorm;
 274                                 ck = MAX0(MIN0(uppb, u[ind]), lowb);
 275
 276                                 // updates of dual var
 277                                 u[ind] -= ck;
 278                                 // backward operation with relaxation param relax
 279                                 for (g = 0; g < numb; g++)
 280                                 {
 281                                         indb = list[g];
 282                                         blob_array[indb] += (ck * weight[g]);
 283                                 }
 284                         } // if no satisf
 285                 } //if numb
 286         } //if c
 287         delete[] list;  // bug 92 - Lajos - 03/02/2005
 288         delete[] weight;  // bug 92 - Lajos - 03/02/2005
 289         return;
 290 }
 291
 292 void Blob_class::p2b_neighb(INTEGER k, INTEGER j, REAL* pict, REAL *nc,
 293                 REAL rneigh)
 294 {
 295
 296         INTEGER m, n;
 297         INTEGER kk, jj, ind;
 298         INTEGER size;
 299         REAL mult = rneigh * Blob.support / GeoPar.pixsiz;
 300         INTEGER multn = (INTEGER) ceil(mult);
 301
 302         size = 0;
 303         (*nc) = 0.0;
 304
 305         for (m = -multn; m <= multn; m++)
 306         {
 307                 for (n = -multn; n <= multn; n++)
 308                 {
 309                         if ((n * n + m * m) <= (mult * mult))
 310                         {
 311                                 size++;
 312                         }
 313                 }
 314         }
 315         for (m = -multn; m <= multn; m++)
 316         {
 317                 kk = k + m;
 318                 if (kk >= 0 && kk < GeoPar.nelem)
 319                 {
 320                         for (n = -multn; n <= multn; n++)
 321                         {
 322                                 jj = j + n;
 323                                 if (jj >= 0 && jj < GeoPar.nelem)
 324                                 {
 325                                         if ((n * n + m * m) <= (mult * mult))
 326                                         {
 327                                                 ind = kk * GeoPar.nelem + jj;
 328                                                 (*nc) += pict[ind];
 329                                         }
 330                                 }
 331                         }
 332                 }
 333         }
 334         (*nc) /= (REAL) size;
 335
 336         return;
 337 }
 338