src/snark/bdhk_misl.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/bdhk_misl.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9  */
  10
  11 #include <cstdio>
  12 #include <cmath>
  13
  14 #include "blkdta.h"
  15 #include "blob.h"
  16 #include "geom.h"
  17 #include "fourie.h"
  18 #include "modefl.h"
  19 #include "consts.h"
  20 #include "uiod.h"
  21 #include "anglst.h"
  22 #include "errfac.h"
  23 #include "errpar.h"
  24 #include "pick.h"
  25 #include "pseudo.h"
  26 #include <cstdio>
  27 #include <iostream>
  28 #include "infile.h"
  29 #include "anglst.h"
  30 #include <iomanip>
  31 #include <fstream>
  32 #include "second.h"
  33 #include "wray.h"
  34 #include "raysel.h"
  35 #include "bdhk.h"
  36 using namespace std;
  37
  38 REAL bdhk_class::AVER(REAL SUM, REAL NUM)
  39 {
  40         return SUM / NUM;
  41 }
  42
  43 REAL bdhk_class::AVER3(REAL A, REAL B, REAL C)
  44 {
  45         REAL temp = (A + B + C) / 3.0;
  46         return temp;
  47 }
  48
  49 REAL bdhk_class::diffAbs(REAL A, REAL B)
  50 {
  51         REAL dif = A - B;
  52
  53         if (dif < 0)
  54                 dif *= -1;
  55
  56         return dif;
  57 }
  58
  59 REAL bdhk_class::neigSum(REAL* recon, REAL* d)
  60 {
  61         REAL sum = 0;
  62         REAL Pi1, Pi2, Pi3, Pi4, Dif1, Dif2, Dif3, Dif4;
  63         for (int i = 0; i < GeoPar.area; i++)
  64         {
  65                 if (i == 0) //upper left corner
  66                 {
  67                         Pi1 = AVER3(recon[1], recon[GeoPar.nelem], recon[GeoPar.nelem + 1]);
  68                         Dif1 = diffAbs(recon[i], Pi1);
  69                         sum += Dif1;
  70                         d[i] = Pi1 - recon[i];
  71                 }
  72                 else if (i == (GeoPar.nelem - 1)) //upper right corner
  73                 {
  74                         Pi1 = AVER3(recon[GeoPar.nelem - 2], recon[2 * GeoPar.nelem - 2],
  75                                         recon[2 * GeoPar.nelem - 1]);
  76                         Dif1 = diffAbs(recon[i], Pi1);
  77                         sum += Dif1;
  78                         d[i] = Pi1 - recon[i];
  79                 }
  80                 else if (i == (GeoPar.area - GeoPar.nelem)) //lower left corner
  81                 {
  82                         Pi1 = AVER3(recon[GeoPar.area - 2 * GeoPar.nelem],
  83                                         recon[GeoPar.area - 2 * GeoPar.nelem + 1],
  84                                         recon[GeoPar.area - GeoPar.nelem + 1]);
  85                         Dif1 = diffAbs(recon[i], Pi1);
  86                         sum += Dif1;
  87                         d[i] = Pi1 - recon[i];
  88                 }
  89                 else if (i == (GeoPar.area - 1)) //lower right  corner
  90                 {
  91                         Pi1 = AVER3(recon[GeoPar.area - 2],
  92                                         recon[GeoPar.area - GeoPar.nelem - 1],
  93                                         recon[GeoPar.area - GeoPar.nelem - 2]);
  94                         Dif1 = diffAbs(recon[i], Pi1);
  95                         sum += Dif1;
  96                         d[i] = Pi1 - recon[i];
  97                 }
  98
  99                 else if (i < GeoPar.nelem) //first row, not corners
 100                 {
 101                         //two possibilities:
 102
 103                         Pi1 = AVER3(recon[i - 1], recon[i - 1 + GeoPar.nelem],
 104                                         recon[i + GeoPar.nelem]);
 105                         Pi2 = AVER3(recon[i + 1], recon[i + 1 + GeoPar.nelem],
 106                                         recon[i + GeoPar.nelem]);
 107
 108                         Dif1 = diffAbs(recon[i], Pi1);
 109                         Dif2 = diffAbs(recon[i], Pi2);
 110
 111                         if (Dif1 < Dif2)
 112                         {
 113                                 sum += Dif1;
 114                                 d[i] = Pi1 - recon[i];
 115
 116                         }
 117                         else
 118                         {
 119                                 sum += Dif2;
 120                                 d[i] = Pi2 - recon[i];
 121                         }
 122                 }
 123
 124                 else if (i % GeoPar.nelem == 0) //first column, not corners
 125                 {
 126                         //two possibilities:
 127
 128                         Pi1 = AVER3(recon[i + 1], recon[i - GeoPar.nelem],
 129                                         recon[i - GeoPar.nelem + 1]);
 130                         Pi2 = AVER3(recon[i + 1], recon[i + GeoPar.nelem],
 131                                         recon[i + GeoPar.nelem + 1]);
 132
 133                         Dif1 = diffAbs(recon[i], Pi1);
 134                         Dif2 = diffAbs(recon[i], Pi2);
 135
 136                         if (Dif1 < Dif2)
 137                         {
 138                                 sum += Dif1;
 139                                 d[i] = Pi1 - recon[i];
 140                         }
 141                         else
 142                         {
 143                                 sum += Dif2;
 144                                 d[i] = Pi2 - recon[i];
 145                         }
 146
 147                 }
 148
 149                 else if (i % GeoPar.nelem == GeoPar.nelem - 1) //last  column, not corners
 150                 {
 151                         //two possibilities:
 152
 153                         Pi1 = AVER3(recon[i - 1], recon[i - GeoPar.nelem - 1],
 154                                         recon[i - GeoPar.nelem]);
 155                         Pi2 = AVER3(recon[i - 1], recon[i + GeoPar.nelem],
 156                                         recon[i + GeoPar.nelem - 1]);
 157
 158                         Dif1 = diffAbs(recon[i], Pi1);
 159                         Dif2 = diffAbs(recon[i], Pi2);
 160
 161                         if (Dif1 < Dif2)
 162                         {
 163                                 sum += Dif1;
 164                                 d[i] = Pi1 - recon[i];
 165                         }
 166                         else
 167                         {
 168                                 sum += Dif2;
 169                                 d[i] = Pi2 - recon[i];
 170                         }
 171
 172                 }
 173                 else if (i > GeoPar.area - GeoPar.nelem) //last  row, not corners
 174                 {
 175                         //two possibilities:
 176
 177                         Pi1 = AVER3(recon[i - 1], recon[i - GeoPar.nelem - 1],
 178                                         recon[i - GeoPar.nelem]);
 179                         Pi2 = AVER3(recon[i + 1], recon[i - GeoPar.nelem],
 180                                         recon[i - GeoPar.nelem + 1]);
 181
 182                         Dif1 = diffAbs(recon[i], Pi1);
 183                         Dif2 = diffAbs(recon[i], Pi2);
 184
 185                         if (Dif1 < Dif2)
 186                         {
 187                                 sum += Dif1;
 188                                 d[i] = Pi1 - recon[i];
 189                         }
 190                         else
 191                         {
 192                                 sum += Dif2;
 193                                 d[i] = Pi2 - recon[i];
 194                         }
 195
 196                 }
 197                 else //not in the bounderies but in the middle
 198                 {
 199                         //four  possibilities:
 200
 201                         Pi1 = AVER3(recon[i - 1], recon[i - GeoPar.nelem - 1],
 202                                         recon[i - GeoPar.nelem]);
 203                         Pi2 = AVER3(recon[i + 1], recon[i - GeoPar.nelem],
 204                                         recon[i - GeoPar.nelem + 1]);
 205                         Pi3 = AVER3(recon[i - 1], recon[i - 1 + GeoPar.nelem],
 206                                         recon[i + GeoPar.nelem]);
 207                         Pi4 = AVER3(recon[i + 1], recon[i + GeoPar.nelem],
 208                                         recon[i + GeoPar.nelem + 1]);
 209
 210                         Dif1 = diffAbs(recon[i], Pi1);
 211                         Dif2 = diffAbs(recon[i], Pi2);
 212                         Dif3 = diffAbs(recon[i], Pi3);
 213                         Dif4 = diffAbs(recon[i], Pi4);
 214
 215                         if (Dif1 <= Dif2 && Dif1 <= Dif3 && Dif1 <= Dif4)
 216                         {
 217                                 sum += Dif1;
 218                                 d[i] = Pi1 - recon[i];
 219                         }
 220                         else if (Dif2 <= Dif1 && Dif2 <= Dif3 && Dif2 <= Dif4)
 221                         {
 222                                 sum += Dif2;
 223                                 d[i] = Pi2 - recon[i];
 224                         }
 225                         else if (Dif3 <= Dif1 && Dif3 <= Dif2 && Dif3 <= Dif4)
 226                         {
 227                                 sum += Dif3;
 228                                 d[i] = Pi3 - recon[i];
 229                         }
 230                         else if (Dif4 <= Dif1 && Dif4 <= Dif2 && Dif4 <= Dif3)
 231                         {
 232                                 sum += Dif4;
 233                                 d[i] = Pi4 - recon[i];
 234                         }
 235                 }
 236
 237         }
 238         REAL norm_d = calcNorm(d);
 239
 240         if (norm_d != 0)
 241         {
 242                 //normalizing the d vector:
 243                 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 244                 for (int i1 = 0; i1 < GeoPar.area; i1++)
 245                         d[i1] /= norm_d;
 246         }
 247
 248         return sum;
 249 }
 250
 251 REAL bdhk_class::calculateResidual(REAL* recon, INTEGER* list, REAL* weight,
 252                 struct NPNR* arrayRays, INTEGER kount)
 253 {
 254         INTEGER np, nr, i, j, numb;
 255         REAL snorm, temp1, b, res = 0;
 256         for (i = 0; i < kount; i++)
 257         {
 258                 np = arrayRays[i].NP;
 259                 nr = arrayRays[i].NR;
 260                 wray(np, nr, list, weight, &numb, &snorm);
 261                 b = Anglst.prdta(np, nr);
 262                 if (numb != 0)
 263                 {
 264                         temp1 = 0;
 265                         for (j = 0; j < numb; j++)
 266                                 temp1 += weight[j] * recon[list[j]];
 267                         res += (pow((b - temp1), 2)) / snorm;
 268                 }       //end if
 269         }       //end for
 270         res = sqrt(res);
 271         return res;
 272 }
 273
 274 BOOLEAN bdhk_class::calculateGradientTV(REAL* grad, REAL* recon)
 275 {
 276         //initialiaze the Gradient vector to all zero
 277         for (int iter1 = 0; iter1 < GeoPar.area; iter1++)
 278                 grad[iter1] = 0;
 279
 280         //looping through case 1
 281         for (int k = 0; k < (GeoPar.area - GeoPar.nelem); k++) //not including the last row
 282         {
 283                 if (((k % GeoPar.nelem) - (GeoPar.nelem - 1)) != 0) //not in the last column
 284                 {
 285
 286                         REAL v1 = recon[k + 1] - recon[k];
 287                         REAL v2 = recon[k + GeoPar.nelem] - recon[k];
 288                         REAL deno = sqrt(v1 * v1 + v2 * v2);
 289                         if (deno >= Consts.zero) //continue calculating the subgradients
 290                                 //there are three subgradients to consider. They are placed in k, k+1, k+nelem
 291                                 // (1) 2x_k - x_k+1 - x_k+nelem / deno
 292                                 grad[k] += (2 * recon[k] - recon[k + 1]
 293                                                 - recon[k + GeoPar.nelem]) / deno;
 294
 295                 }        //end if
 296
 297         }        //end for k
 298
 299                          //looping through case 2
 300         for (int k = 0; k < (GeoPar.area - GeoPar.nelem); k++) //not including the last row
 301         {
 302                 if ((k % GeoPar.nelem) != 0) //not in the first column
 303                 {
 304
 305                         REAL v1 = recon[k] - recon[k - 1];
 306                         REAL v2 = recon[k + GeoPar.nelem - 1] - recon[k - 1];
 307                         REAL deno = sqrt(v1 * v1 + v2 * v2);
 308                         if (deno >= Consts.zero) //continue calculating the subgradients
 309                                 grad[k] += (recon[k] - recon[k - 1]) / deno;
 310
 311                 } //end if
 312
 313         } //end for k
 314
 315           //looping through case 3
 316         for (int k = GeoPar.nelem; k < (GeoPar.area); k++) //not including the first row
 317         {
 318                 if (((k % GeoPar.nelem) - (GeoPar.nelem - 1)) != 0) //not in the last column
 319                 {
 320
 321                         REAL v1 = recon[k + 1 - GeoPar.nelem] - recon[k - GeoPar.nelem];
 322                         REAL v2 = recon[k] - recon[k - GeoPar.nelem];
 323                         REAL deno = sqrt(v1 * v1 + v2 * v2);
 324                         if (deno >= Consts.zero) //continue calculating the subgradients
 325                                 grad[k] += (recon[k] - recon[k - GeoPar.nelem]) / deno;
 326
 327                 } //end if
 328
 329         } //end for k
 330
 331           // BOOLEAN flag=FALSE;
 332         for (int i = 0; i < GeoPar.area; i++)
 333         {
 334                 if (grad[i] != 0)
 335                         return TRUE;
 336         } //end for
 337
 338         return FALSE;
 339
 340 } //end calculateGradient
 341
 342 void bdhk_class::orderRays(struct NPNR* arrayRays, INTEGER kount)
 343 {
 344         INTEGER np, nr, i;
 345         for (i = 0; i < kount; i++)
 346         {
 347                 pick(&np, &nr);
 348                 struct NPNR temp;
 349                 temp.NP = np;
 350                 temp.NR = nr;
 351                 arrayRays[i] = temp;
 352         }
 353 }
 354
 355 void bdhk_class::MyART(REAL* recon, INTEGER* list, REAL* weight,
 356                 struct NPNR* arrayRays, INTEGER kount)
 357 {
 358         INTEGER np, nr, i, j, numb;
 359         REAL snorm, temp, b, div, dif;
 360
 361         for (i = 0; i < kount; i++)
 362         {
 363                 np = arrayRays[i].NP;
 364                 nr = arrayRays[i].NR;
 365                 wray(np, nr, list, weight, &numb, &snorm);
 366                 if (numb != 0)
 367                 {
 368                         if (snorm <= Consts.zero)
 369                         {
 370                                 cout << endl << "*** error ---> norm square is less than "
 371                                                 << Consts.zero << " ***";
 372                                 cout << endl << " *** ART CALCULATION ABORTED ***" << endl;
 373                                 return;
 374                         }
 375                         temp = 0;
 376                         for (j = 0; j < numb; j++)
 377                                 temp += weight[j] * recon[list[j]];
 378                         b = Anglst.prdta(np, nr);
 379                         dif = b - temp;
 380                         div = dif / snorm;
 381                         for (j = 0; j < numb; j++)
 382                                 recon[list[j]] += div * weight[j];
 383
 384                 } //end if
 385         } //end for
 386         return;
 387 } //end MyART
 388
 389 REAL bdhk_class::calcHalfNormSquare(REAL* recon)
 390 {
 391         REAL sum = squareNorm(recon) / 2;
 392         return sum;
 393 }
 394
 395 REAL bdhk_class::calcTV(REAL* recon)
 396 {
 397         REAL sum = 0.0;
 398
 399         for (int x = 0; x < (GeoPar.area - GeoPar.nelem - 1); x++)
 400         {
 401                 if (((x % GeoPar.nelem) - (GeoPar.nelem - 1)) != 0) //not last column
 402                         sum += sqrt(
 403                                         pow((recon[x + 1] - recon[x]), 2)
 404                                                         + pow((recon[x + GeoPar.nelem] - recon[x]), 2));
 405
 406         }
 407
 408         return sum;
 409 }
 410
 411 REAL bdhk_class::squareNorm(REAL* aVector)
 412 {
 413         REAL sum = 0;
 414         for (int i = 0; i < GeoPar.area; i++)
 415                 sum += aVector[i] * aVector[i];
 416         return sum;
 417 }
 418
 419 REAL bdhk_class::calcNorm(REAL* aVector)
 420 {
 421         return sqrt(squareNorm(aVector));
 422 }
 423