src/snark/quad_deltad.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_deltad.cpp $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10                  R#*D*D*R
  11       DELTA = ---------------------
  12               G#*D*(D*M#*A*M*D)*D*G
  13
  14     NOTE: R = M#*A*(M*X-P)
  15
  16     AOPT=1            AOPT=2           AOPT=3             AOPT=4
  17
  18     (R=G=V)           (R=G)            (R=V)
  19
  20 1.                                     PRTN=TWONRM        SAME
  21 2.  TWONRM=D**2*R**2  SAME             SAME               SAME
  22 3.                                     GMA=TWONRM/PRTN    SAME
  23 4.                                     G=R+GMA*G          SAME
  24 5.  CG=D**2*G         SAME             SAME               SAME
  25 6.  V=M#*A*M*CG       SAME             SAME               SAME
  26 7.  GTQG=CG*V         SAME             SAME               SAME
  27           TWONRM
  28 8.  DELTA=------      SAME             SAME               SAME
  29            GTQG
  30
  31                       NOTE:            NOTE:              NOTE:
  32                       TWONRM AND       FOR CONSTANT
  33                       DELTA ARE NOT    DELTA ONLY         SAME
  34                       COMPUTED FOR     V IS COMPUTED
  35                       CONSTANT DELTA
  36 */
  37
  38 #include <cstdio>
  39 #include <cmath>
  40
  41 #include "blkdta.h"
  42 #include "geom.h"
  43 #include "consts.h"
  44 #include "uiod.h"
  45 #include "blob.h"
  46
  47 #include "quad.h"
  48
  49 void quad_class::deltad(
  50 REAL* r,
  51 REAL* g,
  52 REAL* v,
  53 REAL* cg,
  54 REAL* d,
  55 REAL* delta,  // corr. hstau
  56                 REAL* gma,     //corr. hstau
  57                 INTEGER* list,
  58                 REAL* weight,
  59                 BOOLEAN *alg)
  60 {
  61         static REAL twonrm;     // set static. hstau
  62         static REAL prtn;     // set static. hstau
  63         REAL gtqg;
  64         INTEGER i;
  65
  66         INTEGER area;
  67
  68         if (Blob.pix_basis)
  69         {
  70                 area = GeoPar.area;
  71         }
  72         else
  73         {
  74                 area = Blob.area;
  75         }
  76
  77         if (iterx == 1)
  78                 twonrm = 0.0;  //corr? hstau
  79
  80 //C1.
  81
  82         if (!((bopt <= 2) && (popt <= 0)))
  83         {
  84                 prtn = twonrm;
  85
  86 //C2.
  87
  88                 twonrm = 0.0;
  89                 for (i = 0; i < area; i++)
  90                 {
  91                         if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
  92                         {
  93                                 twonrm += d[i] * r[i] * r[i];
  94                         }
  95                 }
  96
  97                 if (popt >= 1)
  98                 {
  99                         fprintf(output, "\n\n      squared norm of the residual is %15.6e",
 100                                         twonrm);     //wei, bug 191, 12/2005
 101                 }
 102
 103                 if (twonrm <= Consts.zero)
 104                 {
 105                         fprintf(output,
 106                                         "\n\n      norm of the residual is zero, iteration stops"); //wei, bug 191, 12/2005
 107                         *alg = TRUE;
 108                         return;
 109                 }
 110         }
 111
 112         if (aopt > 2)
 113         {
 114                 if (iterx == 1)
 115                 {
 116
 117 //C3.
 118                         for (i = 0; i < area; i++)
 119                         {
 120                                 g[i] = r[i];
 121                         }
 122                         goto L4;
 123                         //corr.hstau
 124                 }
 125                 if (bopt == 3)
 126                         *gma = twonrm / prtn; //gma = twonrm/prtn;
 127
 128 //C4.
 129
 130                 for (i = 0; i < area; i++)
 131                 {
 132                         g[i] = r[i] + (*gma) * g[i];  // corr. hstau
 133 ;
 134                 }
 135
 136 //C5.
 137         }
 138         L4:  // corr.hstau
 139         for (i = 0; i < area; i++)
 140         {
 141                 cg[i] = d[i] * g[i];
 142
 143
 144         }
 145
 146         // RETURN FOR AOPT .LE. 2 AND BOPT .LE. 2
 147
 148         if ((bopt <= 2) && (aopt <= 2))
 149                 return; // corr.hstau
 150 //C6.
 151
 152         mtamx(cg, v, list, weight);
 153
 154         if (bopt <= 2)
 155                 return;
 156
 157 //C7.
 158
 159         gtqg = 0.0;
 160         for (i = 0; i < area; i++)
 161         {
 162                 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
 163                 {
 164                         gtqg += cg[i] * v[i];
 165                 }
 166         }
 167
 168
 169         // CHECK FOR ZERO DENOMINATOR
 170
 171         if (fabs(gtqg) > Consts.zero)
 172         {
 173
 174 //C8.
 175                 *delta = twonrm / gtqg;
 176
 177                 return;
 178         }
 179
 180         fprintf(output, "\n***zero denominator when computing delta iteratively, program stops***"); //not selfexplained hstau
 181         *alg = TRUE;
 182         return;
 183 }