src/snark/quad_costfn.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_costfn.cpp $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  THIS SUBROUTINE COMPUTES THE VALUE OF K(X) NEEDED BY THE QUADRATIC
  11  OPTIMIZATION ALGORITHM.
  12  */
  13
  14 #include <cstdio>
  15 #include <cmath>
  16
  17 #include "blkdta.h"
  18 #include "geom.h"
  19 #include "fourie.h"
  20 #include "consts.h"
  21 #include "uiod.h"
  22 #include "wray.h"
  23 #include "ray.h"
  24 #include "anglst.h"
  25 #include "errfac.h"
  26 #include "blob.h"
  27
  28 #include "quad.h"
  29
  30 void quad_class::costfn(REAL* recon, REAL* r, REAL* g, REAL* cg, REAL* w,
  31                 INTEGER* list, REAL* weight, REAL* f)
  32 {
  33
  34         static REAL ppcost; // made static to avoid garbage value for `iterx > 1'. Lajos, Feb 10, 2005
  35         INTEGER np;
  36         REAL theta;
  37         REAL sinth;
  38         REAL costh;
  39         REAL rinc2;
  40         INTEGER nr;
  41         REAL pr;
  42         BOOLEAN alg;
  43         REAL xxcost;
  44         REAL bxcost;
  45         INTEGER i;
  46         REAL xmcost;
  47         REAL pmcost;
  48         INTEGER numb;
  49         REAL snorm;
  50         REAL val;
  51         INTEGER nb;
  52         INTEGER k;
  53         REAL valder;
  54         REAL cost;
  55
  56         INTEGER area;
  57         if (Blob.pix_basis)
  58         {
  59                 area = GeoPar.area;
  60         }
  61         else
  62         {
  63                 area = Blob.area;
  64         }
  65
  66         if (copt == 4)
  67         {
  68                 fprintf(output, "\n\n      cost function can not be computed");
  69                 return;
  70         }
  71
  72         if (iterx <= 1)
  73         {
  74                 ppcost = 0.0;
  75                 Fourie.rinc = GeoPar.pinc;
  76                 for (np = 0; np < GeoPar.prjnum; np++)
  77                 {
  78                         Anglst.getang(np, &theta, &sinth, &costh);
  79                         if (GeoPar.vri)
  80                                 Fourie.rinc = GeoPar.pinc
  81                                                 * MAX0((REAL) fabs(costh), (REAL) fabs(sinth));
  82                         rinc2 = Fourie.rinc * Fourie.rinc;
  83                         for (nr = 0; nr < GeoPar.nrays; nr++)
  84                         {
  85                                 pr = Anglst.prdta(np, nr);
  86                                 if ((eopt - 2) >= 0)
  87                                 {
  88                                         if ((eopt - 2) > 0)
  89                                         {
  90                                                 error_1 = uerror(np, nr, pr, &alg);
  91                                         }
  92                                         else
  93                                         {
  94                                                 error_1 = errfac(pr, Fourie.rinc, rinc2);
  95                                         }
  96                                 }
  97                                 ppcost += pr * pr / error_1;
  98                         }
  99                 }
 100
 101                 n = w;
 102                 if (w == r)
 103                 {
 104                         n = new REAL[area];
 105                 }
 106                 m = n;
 107
 108                 if (!((sb < Consts.zero) || (l <= 0)))
 109                 {
 110                         m = cg;
 111                         if (cg == g)
 112                         {
 113                                 m = new REAL[area];
 114                         }
 115                 }
 116         }
 117
 118         xxcost = 0.0;
 119         bxcost = 0.0;
 120         if ((sb + sc) > Consts.zero)
 121         {
 122                 for (i = 0; i < area; i++)
 123                 {
 124                         if (fopt >= 3)
 125                                 xbar = f[i];
 126                         n[i] = recon[i] - xbar;
 127                         m[i] = n[i];
 128                 }
 129
 130                 if ((sb > Consts.zero) && (l > 0))
 131                 {
 132                         if (Blob.pix_basis)
 133                         {
 134                                 adsmos(n, l, bw1, bw2, bw3, bbcon);
 135                         }
 136                         else
 137                         {
 138                                 badsmos(n, l, bw1, bw2, bbcon);
 139                         }
 140                 }
 141
 142                 for (i = 0; i < area; i++)
 143                 {
 144                         if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
 145                         {
 146                                 xxcost += m[i] * m[i];
 147                                 if (sb > Consts.zero)
 148                                         bxcost += m[i] * n[i];
 149                         }
 150                 }
 151         }
 152
 153         xmcost = 0.0;
 154         pmcost = 0.0;
 155         Fourie.rinc = GeoPar.pinc;
 156         for (np = 0; np < GeoPar.prjnum; np++)
 157         {
 158                 Anglst.getang(np, &theta, &sinth, &costh);
 159                 if (GeoPar.vri)
 160                         Fourie.rinc = GeoPar.pinc * MAX0((REAL) fabs(costh), (REAL) fabs(sinth));
 161                 rinc2 = Fourie.rinc * Fourie.rinc;
 162                 for (nr = 0; nr < GeoPar.nrays; nr++)
 163                 {
 164                         pr = Anglst.prdta(np, nr);
 165                         if (GeoPar.strip)
 166                         {
 167                                 if (Blob.pix_basis)
 168                                 {
 169                                         ray(np, nr, list, weight, &numb, &snorm);
 170                                 }
 171                                 else
 172                                 {
 173                                         Blob.bray(np, nr, list, weight, &numb, &snorm);
 174                                 }
 175                         }
 176                         if (GeoPar.line)
 177                         {
 178                                 if (Blob.pix_basis)
 179                                 {
 180                                         wray(np, nr, list, weight, &numb, &snorm);
 181                                 }
 182                                 else
 183                                 {
 184                                         Blob.bwray(np, nr, list, weight, &numb, &snorm);
 185                                 }
 186                         }
 187                         if (numb != 0)
 188                         {
 189                                 val = 0.0;
 190                                 for (nb = 0; nb < numb; nb++)
 191                                 {
 192                                         k = list[nb];
 193                                         val += recon[k] * weight[nb];
 194                                 }
 195
 196                                 if ((eopt - 2) >= 0)
 197                                 {
 198                                         if ((eopt - 2) > 0)
 199                                         {
 200                                                 error_1 = uerror(np, nr, pr, &alg);
 201                                         }
 202                                         else
 203                                         {
 204                                                 error_1 = errfac(pr, Fourie.rinc, rinc2);
 205                                         }
 206                                 }
 207                                 valder = val / error_1;
 208                                 xmcost += val * valder;
 209                                 pmcost += pr * valder;
 210                         }
 211                 }
 212         }
 213
 214         cost = sa * (ppcost - (REAL) 2.0 * pmcost + xmcost) + sb * bxcost + sc * xxcost;
 215         fprintf(output, "\n\n      value of cost function %15.6e", cost);
 216
 217         return;
 218 }