src/snark/quad_mtamx.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_mtamx.cpp $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  Y = SA*M#*A*M*X
  11
  12  X AND Y MUST BE DISTINCT
  13  */
  14
  15 #include <cstdio>
  16 #include <cmath>
  17
  18 #include "blkdta.h"
  19 #include "geom.h"
  20 #include "noise.h"
  21 #include "wray.h"
  22 #include "anglst.h"
  23 #include "ray.h"
  24 #include "errfac.h"
  25 #include "consts.h"
  26
  27 #include "blob.h"
  28 #include "quad.h"
  29
  30 void quad_class::mtamx(REAL* x, REAL* y, INTEGER* list, REAL* weight)
  31 {
  32         INTEGER i;
  33         INTEGER np;
  34         INTEGER nr;
  35         INTEGER numb;
  36         REAL snorm;
  37         REAL val;
  38         INTEGER nb;
  39         INTEGER k;
  40         REAL truval;
  41         REAL pix4;
  42         REAL theta;
  43         REAL sinth;
  44         REAL costh;
  45
  46         REAL rinc;
  47         REAL rinc2;
  48
  49         BOOLEAN alg;
  50
  51         REAL hex_voronoi;
  52         INTEGER area;
  53
  54         if (Blob.pix_basis)
  55         {
  56                 area = GeoPar.area;
  57         }
  58         else
  59         {
  60                 area = Blob.area;
  61         }
  62
  63         hex_voronoi = Consts.sqrt3 / (REAL) 2.0 * Blob.delta * Blob.delta;
  64
  65         // ZERO OUT Y
  66
  67         for (i = 0; i < area; i++)
  68         {
  69                 y[i] = 0.0;
  70         }
  71
  72         // FOR BETTER PERFORMANCE LINE AND STRIP MODES ARE SEPERATED
  73
  74         if (!GeoPar.strip)
  75         {
  76
  77                 // LINE MODE
  78
  79                 for (np = 0; np < GeoPar.prjnum; np++)
  80                 {
  81                         for (nr = 0; nr < GeoPar.nrays; nr++)
  82                         {
  83                                 if (Blob.pix_basis)
  84                                 {
  85                                         wray(np, nr, list, weight, &numb, &snorm);
  86                                 }
  87                                 else
  88                                 {
  89                                         Blob.bwray(np, nr, list, weight, &numb, &snorm);
  90                                 }
  91                                 if (numb > 0)
  92                                 {
  93                                         val = 0.0;
  94
  95
  96                                         for (nb = 0; nb < numb; nb++)
  97                                         {
  98                                                 k = list[nb];
  99                                                 val += x[k] * weight[nb];
 100                                         }
 101
 102
 103                                         if ((eopt - 2) >= 0)
 104                                         {
 105                                                 if ((eopt - 2) > 0)
 106                                                 {
 107                                                         truval = Anglst.prdta(np, nr);
 108                                                         error_1 = uerror(np, nr, truval, &alg);
 109                                                 }
 110                                                 else
 111                                                 {
 112                                                         truval = Anglst.prdta(np, nr);
 113                                                         error_1 = errfac(truval, 0.0, 0.0);
 114                                                 }
 115                                         }
 116                                         val *= sa / error_1;
 117
 118
 119                                         for (nb = 0; nb < numb; nb++)
 120                                         {
 121                                                 k = list[nb];
 122                                                 y[k] += val * weight[nb];
 123                                         }
 124                                         // NEXT NR
 125                                 }
 126                         }
 127                         // NEXT NP
 128
 129                 }
 130                 return;
 131
 132                 // STRIP MODE
 133         }
 134
 135         rinc = GeoPar.pinc;
 136         if (Blob.pix_basis)
 137         {
 138                 pix4 = sa * GeoPar.pixsiz * GeoPar.pixsiz * GeoPar.pixsiz
 139                                 * GeoPar.pixsiz;
 140         }
 141         else
 142         {
 143                 pix4 = sa * hex_voronoi * hex_voronoi;
 144         }
 145
 146         for (np = 0; np < GeoPar.prjnum; np++)
 147         {
 148                 Anglst.getang(np, &theta, &sinth, &costh);
 149                 if ((NoisePar.quanin > 0) && GeoPar.vri)
 150                 {
 151                         rinc = GeoPar.pinc * MAX0((REAL) fabs(costh), (REAL) fabs(sinth));
 152                 }
 153
 154                 rinc2 = rinc * rinc;
 155
 156                 for (nr = 0; nr < GeoPar.nrays; nr++)
 157                 {
 158                         if (Blob.pix_basis)
 159                         {
 160                                 ray(np, nr, list, weight, &numb, &snorm);
 161                         }
 162                         else
 163                         {
 164                                 Blob.bray(np, nr, list, weight, &numb, &snorm);
 165                         }
 166
 167                         if (numb > 0)
 168                         {
 169                                 val = 0.0;
 170
 171                                 for (nb = 0; nb < numb; nb++)
 172                                 {
 173                                         k = list[nb];
 174                                         val += x[k];
 175                                 }
 176
 177
 178                                 if ((eopt - 2) >= 0)
 179                                 {
 180                                         if ((eopt - 2) > 0)
 181                                         {
 182                                                 truval = Anglst.prdta(np, nr);
 183                                                 error_1 = uerror(np, nr, truval, &alg);
 184                                                 goto L16;
 185                                         }
 186                                         else
 187                                         {
 188                                                 truval = Anglst.prdta(np, nr);
 189                                                 error_1 = errfac(truval, rinc, rinc2);
 190                                         }
 191                                 }
 192
 193                                 L16: val *= pix4 / error_1;
 194
 195
 196                                 for (nb = 0; nb < numb; nb++)
 197                                 {
 198                                         k = list[nb];
 199                                         y[k] += val;
 200                                 }
 201                                 // NEXT NR
 202                         }
 203                 }
 204                 // NEXT NP
 205         }
 206         return;
 207 }