src/snark/quad_badsmos.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/quad_badsmos.cpp $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  THIS ROUTINE PERFORMS THE FOUR DIFFERENT SMOOTHING PROCEDURES
  11  USED IN THE QUADRATIC OPTIMIZATION ALGORITHM FOR BLOB RECONSTRUCTIONS.//changed . hstau
  12
  13  */
  14
  15 #include <cstdio>
  16
  17 #include "blkdta.h"
  18 #include "geom.h"
  19
  20 #include "blob.h"
  21 #include "quad.h"
  22
  23 //pragma message( "to fix!" )
  24 /// fix indexes 1 -> 0
  25
  26 void quad_class::badsmos(REAL* x, INTEGER m, REAL w1, REAL w2, INTEGER bcon)
  27 {
  28         BOOLEAN snar, sajc;
  29
  30         static const INTEGER hsnar = CHAR2INT('s', 'n', 'a', 'r');
  31         static const INTEGER hsajc = CHAR2INT('s', 'a', 'j', 'c');
  32         static const INTEGER hkami = CHAR2INT('k', 'a', 'm', 'i');
  33
  34         INTEGER kami;
  35         //INTEGER n,    // number of rows
  36         INTEGER H;    // number of rows
  37         INTEGER M;   // number of columns
  38         INTEGER m1;  // Blob.M-1
  39         INTEGER m2;  //Blob.M+1
  40         INTEGER h1;  //Blob.H-1
  41
  42         REAL weit4;
  43         REAL weit6;
  44         //REAL weit9;
  45         REAL weit3;
  46         REAL weit5;
  47         REAL weit7;
  48
  49         REAL w3 = w2;
  50         REAL* ip = NULL;
  51         REAL* in = NULL;  //wei 3/2005
  52         //REAL* ipf;         //wei 3/2005
  53         //REAL lb, rb;
  54         REAL lt, rt;
  55         INTEGER j;
  56         INTEGER i;
  57         INTEGER l1;
  58         INTEGER ln;
  59         INTEGER k;
  60         REAL* ipin;
  61         INTEGER li;
  62         INTEGER l;
  63
  64         INTEGER area;
  65         INTEGER kpr;
  66
  67         area = Blob.area;
  68
  69         // FOR THE KAMI METHOD THE ROUTINE IS EXECUTED TWICE
  70
  71         kami = 1;
  72
  73         // INDEX PARAMETERS
  74
  75         //n = GeoPar.nelem;
  76         H = Blob.H;
  77         M = Blob.M;
  78         //n1 = n - 1;
  79         //n2 = n + 1;
  80         h1 = H - 1;
  81         m1 = M - 1;
  82         m2 = M + 1;
  83
  84         snar = FALSE;
  85         sajc = FALSE;
  86         if (bcon == hsnar)
  87                 snar = TRUE;
  88         if (bcon == hsajc)
  89                 sajc = TRUE;
  90         if (snar)
  91         {
  92                 weit4 = (REAL) (w1 + w2 + 2.0 * w3);
  93                 weit6 = weit4 + w2 + w3;
  94                 //weit9 = (REAL) (weit6 + w2 + 2.0 * w3);
  95                 weit3 = (REAL) (w1 + w2 + w3);
  96                 weit5 = weit3 + w2 + w3;
  97                 weit7 = (REAL) (weit5 + 2.0 * w3);
  98
  99         }
 100         // ALLOCATE WORKING SPACE
 101
 102         //ip = new REAL[n];
 103         ip = new REAL[M];
 104
 105         //in = new REAL[n];
 106         in = new REAL[M];
 107
 108         for (;;)
 109         {
 110                 for (j = 0; j < m; j++)
 111                 {
 112
 113                         // SMOOTH FIRST ROW OF PICTURE
 114
 115                         in[Blob.pr] = x[Blob.pr];
 116                         in[m1 - Blob.pr] = x[m1 - Blob.pr];
 117
 118                         if (!sajc)
 119                         {
 120                                 if (Blob.pr == 0)
 121                                 {
 122                                         x[0] = w1 * x[0] + w2 * x[2] + w3 * x[m2];
 123                                         if (snar)
 124                                                 x[Blob.pr] /= weit3; //x[0] /= weit4;
 125                                 }
 126
 127                                 else
 128                                 {
 129                                         x[1] = w1 * x[1] + w2 * x[3] + w3 * x[Blob.pr + m2]
 130                                                         + w3 * x[M];
 131                                         if (snar)
 132                                                 x[Blob.pr] /= weit4; //x[0] /= weit4;
 133                                 }
 134                                 for (i = Blob.pr + 2; i < m1 - Blob.pr; i += 2)
 135                                 {
 136                                         in[i] = x[i];
 137                                         x[i] = w1 * x[i] + w2 * (in[i - 2] + x[i + 2])
 138                                                         + w3 * (x[m1 + i] + x[m2 + i]);
 139
 140                                         if (snar)
 141                                                 x[i] /= weit5; //x[i] /= weit6;
 142                                 }
 143                                 if (Blob.pr == 0)
 144                                 {
 145                                         x[m1] = w1 * x[m1] + w2 * in[m1 - 2] + w3 * x[m1 + m1];
 146                                         //if(snar) x[n-1] /= weit4;
 147                                         if (snar)
 148                                                 x[m1] /= weit3;
 149                                 }
 150                                 else
 151                                 {
 152                                         x[m1 - 1] = w1 * x[m1 - 1] + w2 * in[m1 - 3]
 153                                                         + w3 * x[m1 - 1 + m1] + w3 * x[m1 - 1 + m2];
 154                                         if (snar)
 155                                                 x[m1 - 1] /= weit4;
 156                                 }
 157
 158                         }
 159                         else
 160                         {
 161                                 k = 0;
 162                                 if (Blob.pr == 0)
 163                                 {
 164                                         REAL a1 = x[0];
 165                                         REAL a2 = (x[0] + x[m2] + x[2 * M]) / 3.0;
 166                                         x[0] = (REAL) (w1 * x[0] + w2 * (a1 + x[2])
 167                                                         + w3 * (x[0] + (x[0] + x[2]) / 2.0 + a2 + x[m2]));
 168                                 }
 169                                 else
 170                                 {
 171                                         REAL a1 = (x[1] + x[M]) / 2.0;
 172                                         x[1] = (REAL) (w1 * x[1] + w2 * (a1 + x[3])
 173                                                         + w3
 174                                                                         * (x[1] + (x[1] + x[3]) / 2.0 + x[M]
 175                                                                                         + x[M + 2]));
 176                                 }
 177
 178                                 for (i = Blob.pr + 2; i < m1 - Blob.pr; i += 2)
 179                                 {
 180                                         in[i] = x[i];
 181                                         x[i] = w1 * x[i] + w2 * (in[i - 2] + x[i + 2]) + w3 * (x[m1 + i] + x[m2 + i] + x[i] + (in[i - 2] + x[i + 2]) / 2.0);
 182
 183                                 }
 184
 185                                 // bug 191 - SAJC smoothing - swr - 12/17/05
 186                                 if (Blob.pr == 0)
 187                                 {
 188                                         REAL a1 = x[m1];
 189                                         REAL a2 = (x[m1] + x[m1 + m1] + x[2 * M + m1]) / 3.0;
 190                                         x[m1] = (REAL) (w1 * x[m1] + w2 * (in[m1 - 2] + a1) + w3 * (x[m1] + (x[m1] + in[m1 - 2]) / 2.0 + x[m1 + m1] + a2));
 191                                 }
 192                                 else
 193                                 {
 194                                         REAL a1 = (x[m1 - 1] + x[m1 - 1 + m2]) / 2.0;
 195                                         x[m1 - 1] = (REAL) (w1 * x[m1 - 1] + w2 * (in[m1 - 3] + a1) + w3 * ((in[m1 - 3] + x[m1 - 1]) / 2.0 + x[m1 - 1] + x[m1 + m1 - 1] + x[m1 + m1 + 1]));
 196                                 }
 197
 198                         }
 199
 200                         for (k = 1; k < h1; k++)
 201                         {
 202                                 ipin = in;
 203                                 in = ip;
 204                                 ip = ipin;
 205                                 kpr = (k - Blob.pr) % 2;
 206                                 l = k * M;
 207                                 l1 = l + kpr;
 208                                 ln = l1 + m1 - 2 * kpr;
 209
 210                                 lt = in[kpr];
 211                                 rt = in[m1 - kpr];
 212                                 in[kpr] = x[l1];
 213                                 in[m1 - kpr] = x[ln];
 214
 215                                 x[l1] = w1 * x[l1] + w2 * x[l1 + 2];
 216
 217                                 if (kpr == 0)
 218                                 {
 219                                         x[l1] += w3 * (ip[1] + x[l1 + m2]);
 220                                 }
 221                                 else
 222                                 {
 223                                         x[l1] += w3 * (ip[0] + ip[2] + x[l1 + m2] + x[l1 + m1]);
 224                                 }
 225                                 if (sajc)
 226                                 {
 227                                         // bug 191 - SAJC smoothing - swr - 12/17/05
 228                                         if (kpr == 0)
 229                                         {
 230                                                 REAL a1;
 231                                                 if (k > 1)
 232                                                         a1 = (lt + ip[1] + in[0]) / 3.0;
 233                                                 else
 234                                                         a1 = (ip[1] + in[0]) / 2.0;
 235                                                 REAL a2 = in[0];
 236                                                 REAL a3;
 237                                                 if (k < h1 - 1)
 238                                                         a3 = (in[0] + x[l1 + m2] + x[l1 + 2 * M]) / 3.0;
 239                                                 else
 240                                                         a3 = (in[0] + x[l1 + m2]) / 2.0;
 241                                                 x[l1] += w2 * a2 + w3 * (a1 + a3);
 242                                         }
 243                                         else
 244                                         {
 245                                                 REAL a1 = (ip[0] + in[1] + x[l1 + m1]) / 3.0;
 246                                                 x[l1] += w2 * a1;
 247                                         }
 248                                 }
 249
 250                                 if (snar)
 251                                 {
 252                                         if (kpr == 0)
 253                                                 x[l1] /= weit5;
 254                                         else
 255                                                 x[l1] /= weit6;
 256                                 }
 257                                 for (i = 2 + kpr; i < m1 - kpr; i += 2)
 258                                 {    //next point on the same row
 259                                         li = k * M + i;
 260                                         in[i] = x[li];
 261                                         x[li] = w1 * x[li] + w2 * (in[i - 2] + x[li + 2]) + w3 * (ip[i - 1] + ip[i + 1] + x[li + m1] + x[li + m2]);
 262
 263                                         if (snar)
 264                                                 x[li] /= weit7;
 265                                 }
 266                                 x[ln] = w1 * x[ln] + w2 * in[m1 - 2 - kpr];
 267
 268                                 if (kpr % 2 == 0)
 269                                 {
 270                                         x[ln] += w3 * (ip[m1 - 1] + x[ln + m1]);
 271                                 }
 272                                 else
 273                                 {
 274                                         x[ln] += w3
 275                                                         * (ip[m1 - 2] + x[ln + m1] + x[ln + m2] + ip[m1]);
 276                                 }
 277
 278                                 if (snar)
 279                                 {
 280                                         if (kpr == 0)
 281                                                 x[ln] /= weit4;
 282                                         else
 283                                                 x[ln] /= weit6;
 284                                 }
 285
 286                                 if (sajc)
 287                                 {
 288                                         // bug 191 - SAJC smoothing - swr - 12/17/05
 289                                         if (kpr == 0)
 290                                         {
 291                                                 REAL a1;
 292                                                 if (k > 1)
 293                                                         a1 = (rt + in[m1] + ip[m1 - 1]) / 3.0;
 294                                                 else
 295                                                         a1 = (in[m1] + ip[m1 - 1]) / 2.0;
 296                                                 REAL a2 = in[m1];
 297                                                 REAL a3;
 298                                                 if (k < h1 - 1)
 299                                                         a3 = (in[m1] + x[ln + m1] + x[ln + 2 * M]) / 3.0;
 300                                                 else
 301                                                         a3 = (in[m1] + x[ln + m1]) / 2.0;
 302                                                 x[ln] += w2 * a2 + w3 * (a1 + a3);
 303                                         }
 304                                         else
 305                                         {
 306                                                 REAL a1 = (in[m1 - 1] + ip[m1] + x[ln + m2]) / 3.0;
 307                                                 x[ln] += w2 * a1;
 308                                         }
 309                                 }
 310                         }
 311
 312                         // SMOOTH LAST ROW
 313
 314                         // IN LAST ROW CHANGE ROLE OF IN AND IP
 315
 316                         l1 = h1 * M + Blob.pr;
 317                         ln = l1 + m1 - 2 * Blob.pr;
 318                         l = h1 * M;
 319                         lt = ip[Blob.pr];
 320                         rt = ip[m1 - Blob.pr];
 321
 322                         ip[Blob.pr] = x[l1];
 323
 324                         if (!sajc)
 325                         {
 326                                 if (Blob.pr == 0)
 327                                 {
 328                                         x[l1] = w1 * x[l1] + w2 * x[l1 + 2] + w3 * (in[1]);
 329                                 }
 330                                 else
 331                                 {
 332                                         x[l1] = w1 * x[l1] + w2 * x[l1 + 2] + w3 * (in[0] + in[2]);
 333                                 }
 334
 335                                 if (snar)
 336                                 {
 337                                         if (Blob.pr == 0)
 338                                                 x[l1] /= weit3;
 339                                         else
 340                                                 x[l1] /= weit4;
 341                                 }
 342
 343                                 for (i = 2 + Blob.pr; i < m1 - Blob.pr; i += 2)
 344                                 {
 345                                         li = l + i;
 346                                         ip[i] = x[li];
 347
 348                                         x[li] = w1 * x[li] + w2 * (ip[i - 2] + x[li + 2])
 349                                                         + w3 * (in[i - 1] + in[i + 1]);
 350
 351                                         if (snar)
 352                                                 x[li] /= weit5;
 353                                 }
 354                                 if (Blob.pr == 0)
 355                                 {
 356                                         x[ln] = w1 * x[ln] + w2 * ip[m1 - 2] + w3 * in[m1 - 1];
 357                                 }
 358                                 else
 359                                 {
 360                                         x[ln] = w1 * x[ln] + w2 * ip[m1 - 3]
 361                                                         + w3 * (in[m1 - 2] + in[m1]);
 362                                 }
 363
 364                                 if (snar)
 365                                 {
 366                                         if (Blob.pr == 0)
 367                                                 x[ln] /= weit3;
 368                                         else
 369                                                 x[ln] /= weit4;
 370                                 }
 371                         }
 372                         else
 373                         {
 374                                 // bug 191 - SAJC smoothing - swr - 12/17/05
 375                                 if (Blob.pr == 0)
 376                                 {
 377                                         REAL a1 = x[l1];
 378                                         REAL a2 = (x[l1] + lt + in[1]) / 3.0;
 379                                         x[l1] = w1 * x[l1] + w2 * (x[l1 + 2] + a1)
 380                                                         + w3
 381                                                                         * (a2 + in[1] + x[l1]
 382                                                                                         + (x[l1] + x[l1 + 2]) / 2.0);
 383                                 }
 384                                 else
 385                                 {
 386                                         REAL a1 = (x[l1] + in[0]) / 2.0;
 387                                         x[l1] = w1 * x[l1] + w2 * (a1 + x[l1 + 2]) + w3 * (in[0] + in[2] + x[l1] + (x[l1] + x[l1 + 2]) / 2.0);
 388                                 }
 389
 390                                 for (i = 2 + Blob.pr; i < m1 - Blob.pr; i += 2)
 391                                 {
 392                                         li = l + i;
 393                                         ip[i] = x[li];
 394
 395                                         x[li] = w1 * x[li] + w2 * (ip[i - 2] + x[li + 2])
 396                                                         + w3
 397                                                                         * (x[li] + 0.5 * ip[i - 2] + 0.5 * x[li + 2]
 398                                                                                         + in[i - 1] + in[i + 1]);
 399
 400                                 }
 401
 402                                 // bug 191 - SAJC smoothing - swr - 12/17/05
 403                                 if (Blob.pr == 0)
 404                                 {
 405                                         REAL a1 = x[ln];
 406                                         REAL a2 = (x[ln] + in[m1 - 1] + rt) / 3.0;
 407                                         x[ln] = w1 * x[ln] + w2 * (ip[m1 - 2] + a1) + w3 * (a2 + in[m1 - 1] + (ip[m1 - 2] + x[ln]) / 2.0 + x[ln]);
 408                                 }
 409                                 else
 410                                 {
 411                                         REAL a1 = (x[ln] + in[m1]) / 2.0;
 412                                         x[ln] = w1 * x[ln] + w2 * (ip[m1 - 3] + a1) + w3 * (in[m1 - 2] + in[m1] + (ip[m1 - 3] + x[ln]) / 2.0 + x[ln]);
 413                                 }
 414                         }
 415                 }
 416
 417                 if ((bcon != hkami) || (kami == 2))
 418                         goto L17;
 419
 420                 // PREPARE FOR SECOND ROUND FOR KAMI
 421
 422                 kami = 2;
 423
 424                 for (i = 0; i < M; i++)
 425                 {
 426                         x[i] = 0.0;
 427                         x[l + i] = 0.0;
 428                 }
 429
 430                 l1 = 0;
 431                 ln = m1;
 432
 433                 for (i = 1; i < h1; i++)
 434                 {
 435                         l1 = l1 + M;
 436                         ln = ln + m1;
 437                         x[l1] = 0.0;
 438                         x[ln] = 0.0;
 439                 }
 440         }
 441
 442         L17:
 443         // delete[] ipf;  // bug 92 - Lajos - 03/02/2005      wei/2005
 444         delete[] ip;
 445         delete[] in;   // wei 3/2005
 446
 447         return;
 448 }