src/snark/superior.cpp

   1 /*
   2 $SNARK_Header: S N A R K  1 4 - A PICTURE RECONSTRUCTION PROGRAM $
   3 $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/superior.cpp $
   4 $LastChangedRevision: 158 $
   5 $Date: 2014-07-30 16:59:01 -0400 (Wed, 30 Jul 2014) $
   6 $Author: olangthaler $
   7 */
   8
   9 #include <cstdlib>
  10 #include <iostream>
  11 #include <cmath>
  12 #include <math.h>
  13 #include <stdlib.h>
  14 #include "superior.h"
  15 #include "nav.h"
  16 #include "nav1.h"
  17 #include "nav2.h"
  18 #include "nav3.h"
  19 #include "nav4.h"
  20 #include "nav5.h"
  21 #include "sec_cri.h"
  22 #include "sec_cri1.h"
  23 #include "sec_cri2.h"
  24 #include "sec_cri3.h"
  25 #include "sec_cri4.h"
  26 #include "sec_cri5.h"
  27 #include "uiod.h"
  28 #include "int2str.h"
  29 #include "infile.h"
  30 #include "geom.h"
  31 #include "projfile.h"
  32 #include "blob.h"
  33 #include "DIGRand.h"
  34
  35 using namespace std;
  36
  37 sec_cri1_class sc1;
  38 sec_cri2_class sc2;
  39 sec_cri3_class sc3;
  40 sec_cri4_class sc4;
  41 sec_cri5_class sc5;
  42
  43 sec_cri_class* scClasses[] =
  44 {
  45         &sc1,
  46         &sc2,
  47         &sc3,
  48         &sc4,
  49         &sc5
  50 };
  51
  52 nav1_class nav1;
  53 nav2_class nav2;
  54 nav3_class nav3;
  55 nav4_class nav4;
  56 nav5_class nav5;
  57
  58 nav_class* navClasses[] =
  59 {
  60         &nav1,
  61         &nav2,
  62         &nav3,
  63         &nav4,
  64         &nav5
  65 };
  66
  67 class superior_class SuperSet =
  68 {
  69         FALSE,
  70         0,
  71         0,
  72         0,
  73         NULL,
  74         NULL,
  75         FALSE,
  76         0
  77 };
  78
  79 void initSuperiorization()
  80 {
  81         INTEGER word;
  82         REAL number;
  83         BOOLEAN eol;
  84
  85         // enable superiorization
  86         SuperSet.superiorizationEnabled = TRUE;
  87         cout << "\n         Superiorization is enabled";
  88
  89         // read N
  90         word = InFile.getint(FALSE, &eol);
  91
  92         if (word <= 0)
  93         {
  94                 cout << "\n **** N must be greater than zero";
  95                 cout << "\n **** program aborted\n";
  96                 exit(-1);
  97         }
  98         SuperSet.N = word;
  99         cout << "\n         N = " << SuperSet.N;
 100
 101         // read a
 102         number = InFile.getnum(FALSE, &eol);
 103
 104         if (number <= 0)
 105         {
 106                 cout << "\n **** a must be greater than zero";
 107                 cout << "\n **** program aborted\n";
 108                 exit(-1);
 109         }
 110         if (number >= 1)
 111         {
 112                 cout << "\n **** a must be less than one";
 113                 cout << "\n **** program aborted\n";
 114                 exit(-1);
 115         }
 116         SuperSet.a = number;
 117         cout << "\n         a = " << SuperSet.a;
 118
 119         // read b
 120         number = InFile.getnum(FALSE, &eol);
 121
 122         if (number <= 0)
 123         {
 124                 cout << "\n **** b must be greater than zero";
 125                 cout << "\n **** program aborted\n";
 126                 exit(-1);
 127         }
 128         SuperSet.b = number;
 129         cout << "\n         b = " << SuperSet.b;
 130
 131         // read secondary criterion
 132         static const INTEGER criterion[] =
 133         {
 134                 CHAR2INT('t', 'v', 'a', 'r'),
 135                 CHAR2INT('s', 'm', 'o', 'o'),
 136                 CHAR2INT('s', 'c', 'r', '3'),
 137                 CHAR2INT('s', 'c', 'r', '4'),
 138                 CHAR2INT('s', 'c', 'r', '5')
 139         };
 140
 141         word = InFile.getwrd(FALSE, &eol, criterion, 5);
 142
 143         for (INTEGER i = 0; i < 5; i++)
 144         {
 145                 if (word == criterion[i])
 146                 {
 147                         cout << "\n         secondary criterion: " << int2str(criterion[i]);
 148
 149                         SuperSet.secondaryCriterion = scClasses[i];
 150                         SuperSet.secondaryCriterion->Init();
 151
 152                         // set the default non-ascending vector for the selected secondary criterion
 153                         SuperSet.nonascendingVector = navClasses[i];
 154                         SuperSet.nonascendingVector->Init();
 155
 156                         break;
 157                 }
 158         }
 159
 160         if (SuperSet.secondaryCriterion == NULL)
 161         {
 162                 cout << "\n **** unknown secondary criterion ", int2str(word);
 163                 cout << "\n **** program aborted\n";
 164                 exit(-1);
 165         }
 166
 167         // set default values of the optional commands
 168         SuperSet.report = false;
 169         SuperSet.alternativel=0;
 170         SuperSet.positivity = false;
 171
 172         // read remaining optional commands: posi, altl (alternative l) or rprt
 173         static const INTEGER optionsf[] =
 174         {
 175                 CHAR2INT('a', 't', 'l', '1'),
 176                 CHAR2INT('a', 't', 'l', '2'),
 177                 CHAR2INT('r', 'p', 'r', 't'),
 178                 CHAR2INT('p', 'o', 's', 'i')
 179         };
 180
 181         word = InFile.getwrd(FALSE, &eol, optionsf, 4);
 182
 183         for (INTEGER i = 0; i < 4; i++)
 184         {
 185                 if (word == optionsf[i])
 186                 {
 187                         if (i==0)
 188                         {
 189                                 SuperSet.alternativel=1;
 190                                 cout << "\n         alternative l handling 1 is enabled";
 191                         }
 192                         else if (i==1)
 193                         {
 194                                 SuperSet.alternativel=2;
 195                                 cout << "\n         alternative l handling 2 is enabled";
 196                         }
 197                         else if (i==2)
 198                         {
 199                                 SuperSet.report = true;
 200                                 cout << "\n         reporting is enabled";
 201                                 cout << "\n         reporting file: RPRTsuperiorization";
 202
 203                                 SuperSet.skips = InFile.getnum(FALSE, &eol);
 204                                 if (SuperSet.skips>0)
 205                                 {
 206                                         cout << "\n         report skip factor " << SuperSet.skips;
 207                                 }
 208                                 else
 209                                 {
 210                                         cout << "\n         reporting on every iteration";
 211                                 }
 212                         }
 213                         else if (i==3)
 214                         {
 215                                 SuperSet.positivity=true;
 216                                 cout << "\n         positivity constraint is enabled";
 217                         }
 218                 }
 219         }
 220
 221         // read remaining optional commands: altl (alternative l) or rprt
 222         static const INTEGER options[] =
 223         {
 224                 CHAR2INT('a', 't', 'l', '1'),
 225                 CHAR2INT('a', 't', 'l', '2'),
 226                 CHAR2INT('r', 'p', 'r', 't')
 227         };
 228
 229         word = InFile.getwrd(FALSE, &eol, options, 3);
 230
 231         for (INTEGER i = 0; i < 3; i++)
 232         {
 233                 if (word == options[i])
 234                 {
 235                         if (i==0)
 236                         {
 237                                 SuperSet.alternativel=1;
 238                                 cout << "\n         alternative l handling 1 is enabled";
 239                         }
 240                         else if (i==1)
 241                         {
 242                                 SuperSet.alternativel=2;
 243                                 cout << "\n         alternative l handling 2 is enabled";
 244                         }
 245                         else if (i==2)
 246                         {
 247                                 SuperSet.report = true;
 248                                 cout << "\n         reporting is enabled";
 249                                 cout << "\n         reporting file: RPRTsuperiorization";
 250
 251                                 SuperSet.skips = InFile.getnum(FALSE, &eol);
 252                                 if (SuperSet.skips>0)
 253                                 {
 254                                         cout << "\n         report skip factor " << SuperSet.skips;
 255                                 }
 256                                 else
 257                                 {
 258                                         cout << "\n         reporting on every iteration";
 259                                 }
 260                         }
 261                 }
 262         }
 263
 264         // read last possible command (rprt)
 265         INTEGER rprt = CHAR2INT('r','p','r','t');
 266         if (InFile.getwrd(FALSE, &eol, &rprt, 1) == rprt)
 267         {
 268                 SuperSet.report = true;
 269                 cout << "\n         reporting is enabled";
 270                 cout << "\n         reporting file: RPRTsuperiorization";
 271
 272                 SuperSet.skips = InFile.getnum(FALSE, &eol);
 273                 if (SuperSet.skips>0)
 274                 {
 275                         cout << "\n         report skip factor " << SuperSet.skips;
 276                 }
 277                 else
 278                 {
 279                         cout << "\n         reporting on every iteration";
 280                 }
 281         }
 282
 283         // "set l = -1"
 284         SuperSet.l = -1;
 285 }
 286
 287 BOOLEAN executeSuperiorization(REAL* recon, INTEGER* list, REAL* weight, INTEGER count, alg_class* algorithm)
 288 {
 289         double min = pow(10,-15);
 290         int pixelcount = 0;
 291
 292         if(Blob.pix_basis)
 293         {
 294                 pixelcount = GeoPar.area;
 295         }
 296         else
 297         {
 298                 cout << "\nError: Superiorization does currently not support blobs, skipping Superiorization.";
 299
 300                 // execute one iteration of the selected algorithm
 301                 return algorithm->Run(recon, list, weight, count);
 302         }
 303
 304         // user defined superiorization parameters
 305         INTEGER N = SuperSet.N;
 306         REAL a = SuperSet.a;
 307         REAL b = SuperSet.b;
 308         sec_cri_class* secCrit = SuperSet.secondaryCriterion;
 309         nav_class* naV = SuperSet.nonascendingVector;
 310
 311         // xkn, vkn and z are declared here so that memory is only allocated once
 312         REAL* xkn = new REAL[pixelcount];
 313         REAL* vkn = new REAL[pixelcount];
 314         REAL* z = new REAL[pixelcount];
 315
 316         // things that should only be done on the first iteration
 317         // "set k = 0"
 318         // taken care of by loop in exalg
 319         // "set xk = xdash"
 320         // not necessary (accessing recon directly)
 321         // "set l = -1"
 322         INTEGER l;
 323         if (SuperSet.alternativel==1)
 324         {
 325                 l = count-2;
 326         }
 327         else if (SuperSet.alternativel==2 && count>1)
 328         {
 329                 l = (int)round(Rand()*(SuperSet.l-count)) + count - 1;
 330         }
 331         else if (count==1)
 332         {
 333                 SuperSet.l = l = -1;
 334         }
 335         else
 336         {
 337                 l = SuperSet.l;
 338         }
 339
 340         // "repeat"
 341         // taken care of by loop in exalg
 342         // "set n = 0"
 343         INTEGER n = 0;
 344         // "set xkn = xk"
 345         memcpy(xkn, recon, pixelcount * sizeof(REAL));
 346         // "while n < N"
 347         REAL phixk=secCrit->Run(recon); // only needs to be processed once
 348         while (n < N)
 349         {
 350                 // "set vkn to be a nonascending vector for phi at xkn"
 351                 naV->Run(xkn, vkn);
 352
 353                 // "set loop = true"
 354                 // "while (loop)"
 355
 356                 while (true)
 357                 {
 358                         // "set l = l+1"
 359                         l++;
 360                         // "set betakn = mul"
 361                         REAL betakn = b*pow(a, l);
 362                         if (betakn<min) betakn=0;
 363                         // "set z = xkn + betakn*vkn"
 364                         for (INTEGER i = 0; i < pixelcount; i++) z[i] = xkn[i] + (betakn * vkn[i]);
 365                         // "if z e delta and phi(z) <= phi(xk) then"
 366                         BOOLEAN inDelta = true;
 367                         if (SuperSet.positivity)
 368                         {
 369                                 for (INTEGER i = 0; i < pixelcount; i++) if (z[i] < 0) inDelta=false;
 370                         }
 371                         if (inDelta && secCrit->Run(z) <= phixk)
 372                         {
 373                                 // "set n = n+1"
 374                                 n++;
 375                                 // "set xkn = z"
 376                                 memcpy(xkn, z, pixelcount * sizeof(REAL));
 377                                 // "set loop = false"
 378                                 break;
 379                         }
 380                 }
 381         }
 382         // "set xk+1 = B*xkN"
 383         memcpy(recon, xkn, pixelcount * sizeof(REAL));
 384
 385         delete[] (xkn);
 386         delete[] (vkn);
 387         delete[] (z);
 388
 389         // save value of l for next iteration (checking if this is necessary will take as long as saving it)
 390         SuperSet.l = l;
 391
 392         // reporting output
 393         if (SuperSet.report && (count==1 || SuperSet.skips==0 || (count%(INTEGER)SuperSet.skips)==0))
 394         {
 395                 REAL preAlgo = secCrit->Run(recon);
 396                 BOOLEAN error = algorithm->Run(recon, list, weight, count);
 397                 cout << "\n              value of l: " << l;
 398                 cout << "\n              value of phi before algorithm operator: " << preAlgo;
 399                 REAL postAlgo = secCrit->Run(recon);
 400                 cout << "\n              value of phi after algorithm operator:  " << postAlgo;
 401
 402                 FILE* RPRTfile;
 403                 if (count==1)
 404                 {
 405                         RPRTfile = fopen("RPRTsuperiorization", "w");
 406                         fprintf(RPRTfile,"superiorization reporting output");
 407                         fprintf(RPRTfile,"\n    iter         l     phi pre-algorithm      phi post-algorithm");
 408                         fprintf(RPRTfile,"\n %5i  %10i  %20.9f  %20.9f", count, l, preAlgo, postAlgo);
 409                 }
 410                 else
 411                 {
 412                         RPRTfile = fopen("RPRTsuperiorization", "a");
 413                         fprintf(RPRTfile,"\n %5i  %10i  %20.9f  %20.9f", count, l, preAlgo, postAlgo);
 414                 }
 415                 fclose(RPRTfile);
 416
 417                 return error;
 418         }
 419
 420         // execute one iteration of the selected algorithm
 421         return algorithm->Run(recon, list, weight, count);
 422
 423         // "set k = k+1"
 424         // taken care of by loop in exalg
 425 }