src/snark/DistanceMeasureWSQD.cpp

   1 /*
   2  ***********************************************************
   3  $SNARK_Header$
   4  $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/trm6.cpp $
   5  $LastChangedRevision: 122 $
   6  $Date: 2014-07-09 18:18:59 -0400 (Wed, 09 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9  */
  10
  11 #include <math.h>
  12 #include "DistanceMeasureWSQD.h"
  13 #include "anglst.h"
  14 #include "pseudo.h"
  15 #include "consts.h"
  16
  17 using namespace std;
  18
  19 DistanceMeasureWSQD::DistanceMeasureWSQD() :
  20                 DistanceMeasurePixel("weighted squared distance") {
  21 }
  22 DistanceMeasureWSQD::~DistanceMeasureWSQD() {
  23 }
  24
  25 REAL DistanceMeasureWSQD::calculateDistanceForPixel(REAL* reconPixel,
  26 INTEGER* list, REAL* weight) {
  27
  28         // variables needed during computation of the current residual
  29
  30         REAL theta = 0;
  31         REAL sinth = 0;
  32         REAL costh = 0;
  33
  34         BOOLEAN linef = true;
  35
  36         REAL raysum = 0;
  37         REAL psum = 0;
  38         INTEGER numb = 0;
  39         REAL snorm = 0;
  40         REAL raylen = 0;
  41         REAL wsqd = 0;
  42
  43         // compute current residual
  44         for (int np = 0; np < GeoPar.prjnum; np++) {
  45                 REAL rinc = GeoPar.pinc;
  46                 Anglst.getang(np, &theta, &sinth, &costh);
  47
  48                 if (GeoPar.vri) {
  49                         rinc = GeoPar.pinc * MAX0(fabs(sinth), fabs(costh));
  50                 }
  51
  52                 for (int nr = GeoPar.fusray; nr <= GeoPar.lusray; nr++) {
  53                         raysum = Anglst.prdta(np, nr);
  54                         psum = pseudo(reconPixel, np, nr, list, weight, &numb, &snorm,
  55                                         linef, TRUE);
  56
  57                         if (numb != 0) {
  58                                 // pseudo already evaluates the values for weight/number of affected pixel
  59                                 // --> can be reused
  60                                 // ATTENTION: no other operations on weight/numb are allowed between pseudo and this calculation
  61                                 raylen = 0.0;
  62                                 for (int i = 0; i < numb; i++) {
  63                                         raylen += weight[i];
  64                                 }
  65
  66                                 if (raylen > Consts.zero) {
  67                                         if (linef && !GeoPar.line)
  68                                                 raysum /= rinc;
  69                                         if (!linef && GeoPar.line)
  70                                                 raysum *= rinc;
  71                                         wsqd += pow(psum - raysum, 2) / raylen;
  72                                 }
  73                         }
  74                 }
  75         }
  76
  77         return wsqd;
  78 }