src/snark/trm3.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/trm3.cpp $
   5  $LastChangedRevision: 122 $
   6  $Date: 2014-07-09 18:18:59 -0400 (Wed, 09 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  trm3.cpp,v 1.5 2009/06/01 03:33:26 jklukowska Exp
  11
  12  */
  13
  14 #include <math.h>
  15
  16 #include "blkdta.h"
  17 #include "geom.h"
  18 #include "uiod.h"
  19 #include "consts.h"  // added. hstau
  20 #include "trm3.h"
  21 #include "blob.h"   //wei, 9/05
  22
  23 // bug220 - allow user args for termination tests - swr - 1/30/07
  24
  25 void trm3_class::Init()
  26 {
  27         BOOLEAN eol;
  28
  29         epsiln = InFile.getnum(FALSE, &eol);
  30         // bug 128 - added check for epsilon - swr - 7/6/05
  31         if (epsiln <= Consts_class::zero)
  32         {
  33                 fprintf(output, "\n **** epsilon must be specified and greater than ZERO");
  34                 fprintf(output, "\n **** program aborted\n");
  35                 exit(-1);
  36         }
  37         fprintf(output, "\n         epsilon = %10.6f", epsiln);
  38 }
  39
  40 BOOLEAN trm3_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
  41 {
  42
  43         static REAL mean;
  44         static REAL var1;
  45         static REAL var2;
  46         static REAL chgvar;
  47         REAL area;                       //wei, 9/05
  48
  49         int i;
  50
  51         REAL sum = 0.0;
  52         REAL sqsum = 0.0;
  53
  54         INTEGER lhf, lhl, lmf, lml, h, m, h1, m1, ind;    //wei, 9/05
  55
  56         if (!Blob.pix_basis)
  57         {
  58
  59                 lhf = (INTEGER) (-1 * (INTEGER) (Blob.H / 2));
  60                 lhl = (INTEGER) ((Blob.H - 1) / 2);
  61
  62                 lmf = -1 * (INTEGER) (Blob.M / 2);
  63                 lml = (INTEGER) ((Blob.M - 1) / 2);
  64
  65                 area = Blob.area / 2;
  66
  67                 for (h = lhl; h >= lhf; h--)
  68                 {
  69                         for (m = lmf; m <= lml; m++)
  70                         { // m and h must have the same parity!!
  71                                 if ((m + h) % 2 == 0)
  72                                 {
  73                                         m1 = m + Blob.M2;
  74                                         h1 = h + Blob.H2;
  75                                         ind = h1 * Blob.M + m1;
  76                                         sum += recon[ind];
  77                                         sqsum += recon[ind] * recon[ind];
  78                                 }
  79                         }
  80                 }
  81
  82         } //wei, 9/05
  83
  84         else
  85         {
  86                 area = GeoPar.area;
  87
  88                 for (i = 0; i < area; i++)
  89                 {
  90                         sum += recon[i];
  91                         sqsum += recon[i] * recon[i];
  92                 }
  93         }     //wei, 9/05
  94
  95         mean = sum / area;
  96         var2 = sqsum / area - mean * mean;   //wei, 9/05
  97
  98         if (!(iter == 1))
  99         {
 100                 if (fabs(var1) < Consts.zero)
 101                 {  // added for robustness. hstau
 102                         fprintf(output,
 103                                         "\n          the variance in the previous iteration is less than ZERO");
 104                         return TRUE;
 105                 }
 106                 else
 107                 {             //added for robustness. hsta
 108                         chgvar = ((REAL) fabs(var2 - var1)) / var1;
 109                         if (chgvar < epsiln)
 110                         {
 111                                 fprintf(output,
 112                                                 "\n          iterative process stops at iteration %5i",
 113                                                 iter);
 114                                 fprintf(output,
 115                                                 "\n          the change in variance is less than %10.6f of the variance",
 116                                                 chgvar);
 117                                 return TRUE;
 118                         }
 119                 }
 120         }
 121         var1 = var2;
 122
 123         return FALSE;
 124 }