src/snark/quad_adsmos.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_adsmos.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.
  12  */
  13
  14 #include <cstdio>
  15
  16 #include "blkdta.h"
  17 #include "geom.h"
  18
  19 #include "blob.h"
  20 #include "quad.h"
  21
  22 //pragma message( "to fix!" )
  23 /// fix indexes 1 -> 0
  24
  25 void quad_class::adsmos(REAL* x, INTEGER m, REAL w1, REAL w2, REAL w3,
  26                 INTEGER bcon)
  27 {
  28
  29         BOOLEAN snar, sajc;
  30
  31         static const INTEGER hsnar = CHAR2INT('s', 'n', 'a', 'r');
  32         static const INTEGER hsajc = CHAR2INT('s', 'a', 'j', 'c');
  33         static const INTEGER hkami = CHAR2INT('k', 'a', 'm', 'i');
  34
  35         INTEGER kami;
  36         INTEGER n;
  37         INTEGER n1;
  38         INTEGER n2;
  39         REAL weit4;
  40         REAL weit6;
  41         REAL weit9;
  42         REAL* ip = NULL;
  43         REAL* in = NULL;   //wei  3/2005
  44         // REAL* ipf;          //wei 3/2005
  45         INTEGER j;
  46         INTEGER i;
  47         INTEGER l1;
  48         INTEGER ln;
  49         INTEGER k;
  50         REAL* ipin;
  51         INTEGER li;
  52         INTEGER l;
  53
  54         // FOR THE KAMI METHOD THE ROUTINE IS EXECUTED TWICE
  55
  56         kami = 1;
  57
  58         // INDEX PARAMETERS
  59
  60         n = GeoPar.nelem;
  61         n1 = n - 1;
  62         n2 = n + 1;
  63         snar = FALSE;
  64         sajc = FALSE;
  65         if (bcon == hsnar)
  66                 snar = TRUE;
  67         if (bcon == hsajc)
  68                 sajc = TRUE;
  69         if (snar)
  70         {
  71                 weit4 = (REAL) (w1 + 2.0 * w2 + w3);
  72                 weit6 = weit4 + w2 + w3;
  73                 weit9 = (REAL) (weit6 + w2 + 2.0 * w3);
  74
  75         }
  76         // ALLOCATE WORKING SPACE
  77
  78         ip = new REAL[n];
  79
  80         in = new REAL[n];
  81
  82         // IPF IS  USE FOR SUB FREE ONLY
  83
  84         //ipf = ip;         wei, 3/2005
  85
  86         // SMOOTHING IS DONE M TIMES
  87
  88         for (;;)
  89         {
  90                 for (j = 0; j < m; j++)
  91                 {
  92
  93                         // SMOOTH FIRST ROW OF PICTURE
  94
  95                         in[0] = x[0];
  96                         in[n1] = x[n1];
  97
  98                         if (!sajc)
  99                         {
 100
 101                                 x[0] = w1 * x[0] + w2 * (x[1] + x[n]) + w3 * x[n2];
 102                                 if (snar)
 103                                         x[0] /= weit4;
 104
 105                                 for (i = 1; i < n1; i++)
 106                                 {
 107
 108                                         in[i] = x[i];
 109
 110                                         x[i] = w1 * x[i] + w2 * (in[i - 1] + x[i + 1] + x[n + i])
 111                                                         + w3 * (x[n1 + i] + x[n2 + i]);
 112
 113                                         if (snar)
 114                                                 x[i] /= weit6;
 115                                 }
 116                                 x[n1] = w1 * x[n1] + w2 * (in[n1 - 1] + x[n1 + n])
 117                                                 + w3 * x[n + n1 - 1];
 118
 119                                 if (snar)
 120                                         x[n - 1] /= weit4;
 121                         }
 122                         else
 123                         {
 124                                 x[0] = (REAL) (w1 * x[0] + w2 * (2.0 * x[0] + x[1] + x[n])
 125                                                 + w3 * (x[0] + x[1] + x[n2] + x[n]));
 126
 127                                 for (i = 1; i < n1; i++)
 128                                 {
 129                                         in[i] = x[i];
 130
 131                                         x[i] = w1 * x[i] + (w2 + w3) * (in[i - 1] + x[i + 1])
 132                                                         + w2 * (x[i] + x[n + i])
 133                                                         + w3 * (x[n1 + i] + x[n2 + i]);
 134
 135                                 }
 136
 137                                 x[n1] = (REAL) (w1 * x[n1]
 138                                                 + w2 * (in[n1 - 1] + 2. * x[n1] + x[n + n1])
 139                                                 + w3 * (x[n1] + in[n1 - 1] + x[n1 + n] + x[n1 + n1]));
 140                         }
 141
 142                         // ROW 2 TO ROW N - 1
 143
 144                         l = n;
 145                         //l1 = n2;
 146                         l1 = n;
 147                         //ln = 2*n;
 148                         ln = 2 * n - 1;
 149
 150                         for (k = 1; k < n1; k++)
 151                         {
 152                                 ipin = in;
 153                                 in = ip;
 154                                 ip = ipin;
 155                                 in[0] = x[l1];
 156                                 in[n1] = x[ln];
 157                                 x[l1] = w1 * x[l1] + w2 * (ip[0] + x[l1 + 1] + x[n + l1])
 158                                                 + w3 * (ip[1] + x[l1 + n2]);
 159
 160                                 if (sajc)
 161                                         x[l1] += w2 * (in[0]) + w3 * (ip[0] + x[l1 + n1]);
 162                                 if (snar)
 163                                         x[l1] /= weit6;
 164
 165                                 for (i = 1; i < n1; i++)
 166                                 {
 167
 168                                         li = l + i;
 169                                         in[i] = x[li];
 170
 171                                         x[li] = w1 * x[li]
 172                                                         + w2 * (in[i - 1] + x[li + 1] + ip[i] + x[li + n])
 173                                                         + w3
 174                                                                         * (ip[i - 1] + ip[i + 1] + x[li + n1]
 175                                                                                         + x[li + n2]);
 176
 177                                         if (snar)
 178                                                 x[li] /= weit9;
 179                                 }
 180
 181                                 x[ln] = w1 * x[ln] + w2 * (in[n - 2] + ip[n - 1] + x[ln + n])
 182                                                 + w3 * (ip[n - 2] + x[ln + n1]);
 183
 184                                 if (snar)
 185                                         x[ln] /= weit6;
 186                                 if (sajc)
 187                                         x[ln] += w2 * in[n1] + w3 * (ip[n1] + x[ln + n]);
 188
 189                                 l = l + n;
 190                                 l1 = l1 + n;
 191                                 ln = ln + n;
 192                         }
 193
 194                         // SMOOTH LAST ROW
 195
 196                         // IN LAST ROW CHANGE ROLE OF IN AND IP
 197
 198                         //ip[1] = x[l1];
 199                         ip[0] = x[l1];
 200
 201                         if (!sajc)
 202                         {
 203
 204                                 x[l1] = w1 * x[l1] + w2 * (in[0] + x[l1 + 1]) + w3 * (in[1]);
 205
 206                                 if (snar)
 207                                         x[l1] /= weit4;
 208
 209                                 for (i = 1; i < n1; i++)
 210                                 {
 211                                         li = l + i;
 212                                         ip[i] = x[li];
 213
 214                                         x[li] = w1 * x[li] + w2 * (in[i] + ip[i - 1] + x[li + 1])
 215                                                         + w3 * (in[i - 1] + in[i + 1]);
 216
 217                                         if (snar)
 218                                                 x[li] /= weit6;
 219                                 }
 220
 221                                 x[ln] = w1 * x[ln] + w2 * (ip[n1 - 1] + in[n1])
 222                                                 + w3 * in[n1 - 1];
 223
 224                                 if (snar)
 225                                         x[ln] /= weit4;
 226                         }
 227                         else
 228                         {
 229                                 x[l1] = (REAL) (w1 * x[l1]
 230                                                 + w2 * (2.0 * x[l1] + in[0] + x[l1 + 1])
 231                                                 + w3 * (x[l1] + in[1] + in[0] + x[l1 + 1]));
 232
 233                                 for (i = 1; i < n1; i++)
 234                                 {
 235                                         li = l + i;
 236
 237                                         ip[i] = x[li];
 238
 239                                         x[li] = w1 * x[li] + (w2 + w3) * (ip[i - 1] + x[li + 1])
 240                                                         + w2 * (in[i] + x[li])
 241                                                         + w3 * (in[i - 1] + in[i + 1]);
 242
 243                                 }
 244
 245                                 x[ln] = (REAL) (w1 * x[ln]
 246                                                 + w2 * (in[n1] + ip[n1 - 1] + 2.0 * x[ln])
 247                                                 + w3 * (x[ln] + in[n - 1] + in[n1 - 1] + ip[n1 - 1]));
 248                         }
 249                 }
 250
 251                 if ((bcon != hkami) || (kami == 2))
 252                         goto L17;
 253
 254                 // PREPARE FOR SECOND ROUND FOR KAMI
 255
 256                 kami = 2;
 257
 258                 for (i = 0; i < n; i++)
 259                 {
 260                         x[i] = 0.0;
 261                         x[l + i] = 0.0;
 262                 }
 263
 264                 l1 = 0;
 265                 ln = n;
 266                 for (i = 2; i < n; i++)
 267                 {
 268                         l1 = l1 + n;
 269                         ln = ln + n;
 270                         x[l1] = 0.0;
 271                         x[ln] = 0.0;
 272                 }
 273         }
 274
 275         L17:
 276         //delete[] ipf;  // bug 92 - Lajos - 03/02/2005   wei,3/2005
 277         delete[] ip;
 278         delete[] in;  // wei 3/2005
 279
 280         return;
 281 }