src/snark/dist.cpp

   1 /*
   2  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
   3  *                                                                           *
   4  *                              S N A R K   1 4                              *
   5  *                                                                           *
   6  *                     A PICTURE RECONSTRUCTION PROGRAM                      *
   7  *                                                                           *
   8  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
   9
  10  dist.cpp,v 1.2 2008/08/25 16:11:08 jklukowska Exp
  11
  12  COMPARE THE RECONSTRUCTION WITH THE TEST PICTURE
  13  ***** INPUT VARIABLES *****
  14  A,B:     TWO NELEM*NELEM ARRAYS CONTAINING THE TWO PICTURES
  15  TO BE COMPARED.
  16  LINEF:    A LOGICAL VARIABLE WHICH DETERMINES IF THE RESIDUAL
  17  IS TO BE COMPUTED IN LINE OR STRIP MODE.
  18  LIST,WEIGHT: ARRAYS NECESSARY FOR RAY AND WRAY.
  19  FLAG:     INDICATOR TO CONTROL RESIDUAL COMPUTATION.
  20  FLAG.NE.1 CAUSES RESIDUAL TO BE COMPUTED.
  21  ***** OUTPUT VARIABLES *****
  22  AVE(J):      THE AVERAGE DENSITY IN REGION J OF PICTURE B.
  23  DIS(J):      THE TWO NORM DISTANCE IN REGION J BETWEEN A AND B.
  24  REL(J):      THE ONE NORM DISTANCE IN REGION J BETWEEN A AND B.
  25  VAR(J):      THE VARIANCE IN REGION J OF PICTURE B.
  26  STD(J):      THE STANDARD DEVIATION IN REGION J OF B (=SQRT(VAR)).
  27  RESID:    THE RESIDUAL- I.E. THE TWO NORM OF THE DIFFERENCE
  28  BETWEEN THE PSEUDO PROJECTION OF B AND THE GIVEN
  29  PROJECTION DATA FOR THE WHOLE PICTURE.
  30
  31  ---  Modified 4/1/88 to enable evaluation of subregions, either
  32  through boundary selection by geomtric descriptions OR (OR BOTH)
  33  pixel selection through density windowing wrt to phantom.
  34  ---  Only pixels within the specified regions and density window
  35  will be used in computing the metrics
  36  ---  A switch ISW is added to indicate whether the phantom is passed through
  37  the first array, a(), or the second array, b(), so that density
  38  selection can always be referenced wrt the phantom
  39  */
  40
  41 #include <cstdlib>
  42 #include <cstdio>
  43 #include <cmath>
  44 #include <stdio.h>
  45 #include <iostream>
  46 #include <limits>
  47 #include "math.h"
  48 #include "blkdta.h"
  49 #include "geom.h"
  50 #include "region.h"
  51 #include "anglst.h"
  52 #include "pseudo.h"
  53 #include "blob.h"
  54 #include "consts.h"
  55 #include "anglst.h"
  56 #include "creacm.h"
  57 #include "eval.h"
  58 #include "DistanceMeasureWSQD.h"
  59
  60 void Eval_class::dist(
  61 INTEGER isw,
  62 REAL* a,
  63 REAL* b,
  64 BOOLEAN linef,
  65 INTEGER* list,
  66 REAL* weight,
  67 INTEGER flag,
  68 REAL* ave,
  69 REAL* dis,
  70 REAL* rel,
  71 REAL* var,
  72 REAL* std,
  73 REAL* resid,
  74 BOOLEAN fresid,
  75 REAL* klds,
  76 BOOLEAN fklds,
  77 REAL* wsqd,
  78 BOOLEAN fwsqd)
  79 {
  80         REAL xd;
  81         REAL yd;
  82         REAL arg;
  83
  84         REAL theta, sinth, costh;
  85
  86         REAL raysum;
  87         REAL psum;
  88         INTEGER numb;
  89         REAL snorm;
  90         REAL rinc;
  91         REAL drecon;
  92         REAL dfantom;
  93         REAL dtmp;
  94         REAL divis;
  95         REAL t1;
  96         REAL t2;
  97         REAL sumx;
  98         REAL rel1;
  99         REAL sumx2;
 100
 101         INTEGER npt1;
 102         INTEGER rarea;
 103         INTEGER i1;
 104         INTEGER idx0;
 105         INTEGER objtyp;
 106         REAL rdis;
 107         REAL wsqds;
 108
 109         // The near0 parameter is added to allow resetting non-zero variances
 110         // (presumably arising from rounding error), computed variance < near0
 111         // will be set to zero
 112
 113         REAL near0 = (REAL) 1e-5;
 114
 115         //  Repeat for other specified regions until all are considered
 116         for (INTEGER ObjJ = 0; ObjJ < ObjN; ObjJ++)
 117         {
 118
 119                 t1 = Obj[ObjJ].t1;
 120                 t2 = Obj[ObjJ].t2;
 121                 sumx = 0.0;
 122                 rel1 = 0.0;
 123                 sumx2 = 0.0;
 124                 rdis = 0;
 125                 wsqds = 0.0;
 126                 divis = 0.0;
 127                 npt1 = 0;
 128
 129                 objtyp = Obj[ObjJ].objcod;
 130
 131                 // Specified region is the whole picture
 132                 if (objtyp == 0)
 133                 {
 134
 135                         // Test if pixel density (of phantom) is within selected range,
 136                         // pixel outside density range are excluded
 137                         for (int i = 0; i < GeoPar.area; i++)
 138                         {
 139
 140                                 if (isw == 0)
 141                                 {
 142                                         dtmp = a[i];
 143                                 }
 144                                 else
 145                                 {
 146                                         dtmp = b[i];
 147                                 }
 148
 149                                 if ((dtmp >= t1) && (dtmp <= t2))
 150                                 {
 151                                         dfantom = a[i];
 152                                         drecon = b[i];
 153                                         sumx += drecon;
 154                                         sumx2 += SQR(drecon);
 155                                         drecon = (REAL) fabs(dfantom - drecon);
 156                                         rel1 += drecon;
 157                                         divis += drecon * drecon;
 158                                         npt1++;
 159                                 }
 160                         }
 161                 }
 162                 else
 163                 {
 164
 165                         // Call subroutine "subreg" to compute parameters of specified region.
 166                         subreg(ObjJ);
 167
 168                         idx0 = (Region.iy1 - 2) * GeoPar.nelem;
 169
 170                         for (int j = Region.iy1; j <= Region.iy2; j++)
 171                         {
 172                                 xd = Region.xd0;
 173                                 yd = Region.yd0;
 174                                 idx0 += GeoPar.nelem;
 175
 176                                 // Test if current pixel is in specified region
 177                                 for (int i = Region.ix1; i <= Region.ix2; i++)
 178                                 {
 179                                         if (objtyp == 1)
 180                                         {
 181                                                 // Rectangular subregion
 182                                                 if ((fabs(xd) > Region.aplus)
 183                                                                 || (fabs(yd) > Region.bplus))
 184                                                         goto L222;
 185                                         }
 186                                         else
 187                                         {
 188                                                 // Elliptical subregion
 189                                                 arg = (xd * xd) / Region.asq + (yd * yd) / Region.bsq;
 190                                                 if (arg > 1.0001)
 191                                                         goto L222;
 192                                         }
 193
 194                                         // Test if pixel density is within selected range
 195                                         i1 = idx0 + i;
 196
 197                                         if (isw == 0)
 198                                         {
 199                                                 dtmp = a[i1];
 200                                         }
 201                                         else
 202                                         {
 203                                                 dtmp = b[i1];
 204                                         }
 205
 206                                         if ((dtmp >= t1) && (dtmp <= t2))
 207                                         {
 208                                                 dfantom = a[i1];
 209                                                 drecon = b[i1];
 210                                                 sumx += drecon;
 211                                                 sumx2 += SQR(drecon); // sumx2 += SQR(drecon - ave[ObjJ]);
 212                                                 drecon = (REAL) fabs(dfantom - drecon);
 213                                                 rel1 += drecon;
 214                                                 divis += drecon * drecon;
 215                                                 npt1++;
 216                                         }
 217
 218                                         L222: xd += Region.d1;
 219                                         yd -= Region.d2;
 220                                 }
 221
 222                                 Region.xd0 -= Region.d2;
 223                                 Region.yd0 -= Region.d1;
 224                         }
 225                 }
 226
 227                 rel[ObjJ] = rel1;
 228                 Obj[ObjJ].npt = npt1;
 229                 rarea = npt1;
 230                 ave[ObjJ] = sumx / rarea;
 231                 dis[ObjJ] = (REAL) sqrt(divis / rarea);
 232                 var[ObjJ] = MAX0((REAL) 0.0, sumx2 / rarea - SQR(ave[ObjJ]));
 233
 234                 // Reset very small non-zero variance (presumably due to rounding
 235                 // error) to zero
 236
 237                 if (var[ObjJ] <= near0)
 238                 {
 239                         var[ObjJ] = 0;
 240                 }
 241
 242                 std[ObjJ] = (REAL) sqrt(var[ObjJ]);
 243         }
 244
 245         // ALLOW USER TO SKIP EXPENSIVE STATISTIC COMPUTATION
 246         if ((fresid || fklds || fwsqd) && flag != 1)
 247         {
 248                 rdis = 0;
 249                 REAL kldis = 0;
 250                 bool kldserrorreported=false;
 251                 for (int np = 0; np < GeoPar.prjnum; np++)
 252                 {
 253                         rinc = GeoPar.pinc;
 254                         Anglst.getang(np, &theta, &sinth, &costh);
 255
 256                         if (GeoPar.vri)
 257                         {
 258                                 rinc = GeoPar.pinc * MAX0((REAL) fabs(sinth), (REAL) fabs(costh));
 259                         }
 260
 261                         for (int nr = GeoPar.fusray; nr <= GeoPar.lusray; nr++)
 262                         {
 263                                 raysum = Anglst.prdta(np, nr);
 264                                 psum = pseudo(b, np, nr, list, weight, &numb, &snorm, linef, TRUE);
 265
 266                                 if (numb != 0)
 267                                 {
 268                                         if (linef && !GeoPar.line) raysum /= rinc;
 269                                         if (!linef && GeoPar.line) raysum *= rinc;
 270                                         rdis += SQR(psum - raysum);
 271
 272                                         // calculate kullback-leibler distance
 273                                         if (!fklds)
 274                                         {
 275                                                 // if it wasn't requested, set it to nan
 276                                                 kldis=nan("");
 277                                         }
 278                                         else if (raysum<0 && !kldserrorreported)
 279                                         {
 280                                                 // report on non-applicability
 281                                                 std::cout << "\n****Warning: Negative raysum in prjfil detected, KLDS is not applicable and will be nan!";
 282                                                 kldis=nan("");
 283                                                 kldserrorreported=true;
 284                                         }
 285                                         // if it was requested and no errors occured, perform calculation
 286                                         else if (!kldserrorreported && !isnan(kldis) && kldis!=std::numeric_limits<double>::infinity() && kldis!=std::numeric_limits<double>::max())
 287                                         {
 288                                                 if (raysum<=Consts.zero && psum<=Consts.zero)
 289                                                 {
 290                                                         // do nothing
 291                                                 }
 292                                                 else if (raysum<=Consts.zero)
 293                                                 {
 294                                                         kldis += psum;
 295                                                 }
 296                                                 else if (psum<=Consts.zero)
 297                                                 {
 298                                                         // set to infinity
 299                                                         if (std::numeric_limits<double>::has_infinity) kldis = std::numeric_limits<double>::infinity();
 300                                                         else kldis = std::numeric_limits<double>::max();
 301                                                 }
 302                                                 else
 303                                                 {
 304                                                         if ((raysum / psum) > Consts.zero) kldis += (raysum * log(raysum / psum));
 305                                                         else kldis += (raysum * log(Consts.zero));
 306
 307                                                         kldis += psum - raysum;
 308                                                 }
 309                                         }
 310                                 }
 311                         }
 312                 }
 313
 314                 *resid = (REAL) sqrt(rdis);
 315                 *klds = kldis;
 316
 317                 if (fwsqd) {
 318                         DistanceMeasureWSQD* distanceMeasureWSQD = new DistanceMeasureWSQD();
 319                         *wsqd = distanceMeasureWSQD->calculateDistance(b, list, weight);
 320                         delete distanceMeasureWSQD;
 321                 }
 322         }
 323         return;
 324 }