src/snark/quad_mtamxp.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_mtamxp.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 - P)
  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 "fourie.h"
  21 #include "noise.h"
  22 #include "anglst.h"
  23 #include "wray.h"
  24 #include "ray.h"
  25 #include "errfac.h"
  26 #include "consts.h"
  27
  28 #include "blob.h"
  29 #include "quad.h"
  30
  31 void quad_class::mtamxp(REAL* x, REAL* y, INTEGER* list, REAL* weight)
  32 {
  33         BOOLEAN alg;
  34
  35         INTEGER i;
  36         INTEGER np;
  37         INTEGER nr;
  38         INTEGER numb;
  39         REAL snorm;
  40         REAL val;
  41         INTEGER nb;
  42         INTEGER k;
  43         REAL truval;
  44         REAL pix2;
  45         REAL pix3;
  46         REAL pix4;
  47         REAL theta;
  48         REAL sinth;
  49         REAL costh;
  50         REAL rinc2;
  51         REAL hex_voronoi;
  52         hex_voronoi = Consts.sqrt3 / (REAL) 2.0 * Blob.delta * Blob.delta;
  53
  54         INTEGER area;
  55
  56         if (Blob.pix_basis)
  57         {
  58
  59                 area = GeoPar.area;
  60         }
  61         else
  62         {
  63                 area = Blob.area;
  64         }
  65
  66         ///tracer
  67         //for(i = 0; i < area; i++) {
  68         //   fprintf(output,"\nmtamxpx=%f", x[i]);
  69         // }
  70         ////tracer
  71
  72         // ZERO OUT Y
  73
  74         for (i = 0; i < area; i++)
  75         {
  76                 y[i] = 0.0;
  77         }
  78
  79         // FOR BETTER PERFORMANCE LINE AND STRIP MODES ARE SEPERATED
  80
  81         if (!GeoPar.strip)
  82         {
  83
  84                 // LINE MODE
  85
  86                 for (np = 0; np < GeoPar.prjnum; np++)
  87                 {
  88                         for (nr = 0; nr < GeoPar.nrays; nr++)
  89                         {
  90
  91                                 if (Blob.pix_basis)
  92                                 {
  93                                         wray(np, nr, list, weight, &numb, &snorm);
  94                                 }
  95                                 else
  96                                 {
  97                                         Blob.bwray(np, nr, list, weight, &numb, &snorm);
  98                                 }
  99                                 if (numb > 0)
 100                                 {
 101                                         val = 0.0;
 102
 103                                         // VAL = M*X
 104
 105                                         for (nb = 0; nb < numb; nb++)
 106                                         {
 107                                                 k = list[nb];
 108                                                 val += x[k] * weight[nb];
 109                                         }
 110
 111                                         // VAL = M*X - P
 112
 113                                         truval = Anglst.prdta(np, nr);
 114                                         val -= truval;
 115
 116                                         // VAL = SA*A*(M*X - P)
 117
 118                                         if ((eopt - 2) >= 0)
 119                                         {
 120                                                 if ((eopt - 2) > 0)
 121                                                 {
 122                                                         error_1 = uerror(np, nr, truval, &alg);
 123                                                 }
 124                                                 else
 125                                                 {
 126                                                         error_1 = errfac(truval, 0.0, 0.0);
 127                                                 }
 128                                         }
 129
 130                                         val *= sa / error_1;
 131
 132                                         // Y = SA*M#*A*(M*X - P)
 133
 134                                         for (nb = 0; nb < numb; nb++)
 135                                         {
 136                                                 k = list[nb];
 137                                                 y[k] += val * weight[nb];
 138                                         }
 139
 140 //    NEXT NR
 141                                 }
 142                         }
 143
 144 //    NEXT NP
 145
 146                 }
 147
 148                 return;
 149         }
 150 //    STRIP MODE
 151
 152         Fourie.rinc = GeoPar.pinc;  //think rinc and these pix* for blob case. hstau
 153         if (Blob.pix_basis)
 154         {
 155                 pix2 = GeoPar.pixsiz * GeoPar.pixsiz;
 156                 pix3 = pix2 * sa;
 157                 pix4 = pix3;   // bug in Fortran?. NO hstau
 158         }
 159         else
 160         {
 161                 pix2 = hex_voronoi;
 162                 pix3 = pix2 * sa;
 163                 pix4 = pix3;
 164         }
 165
 166         for (np = 0; np < GeoPar.prjnum; np++)
 167         {
 168                 Anglst.getang(np, &theta, &sinth, &costh);
 169                 if ((NoisePar.quanin > 0) && GeoPar.vri)
 170                         Fourie.rinc = GeoPar.pinc
 171                                         * MAX0((REAL) fabs(costh), (REAL) fabs(sinth));
 172                 rinc2 = Fourie.rinc * Fourie.rinc;
 173
 174                 for (nr = 0; nr < GeoPar.nrays; nr++)
 175                 {
 176                         if (Blob.pix_basis)
 177                         {
 178                                 ray(np, nr, list, weight, &numb, &snorm);
 179                         }
 180                         else
 181                         {
 182                                 Blob.bray(np, nr, list, weight, &numb, &snorm);
 183                         }
 184
 185                         if (numb > 0)
 186                         {
 187                                 val = 0.0;
 188
 189                                 // VAL = M*X  (NOTE THAT M IS A 0-1 MATRIX)
 190
 191                                 for (nb = 0; nb < numb; nb++)
 192                                 {
 193                                         k = list[nb];
 194                                         val += x[k];
 195                                 }
 196
 197                                 val *= pix2;
 198
 199                                 // VAL = M*X - P
 200
 201                                 truval = Anglst.prdta(np, nr);
 202
 203                                 val -= truval;
 204
 205                                 // VAL = SA*M#*A*(M*X - P)
 206
 207                                 if ((eopt - 2) >= 0)
 208                                 {
 209                                         if ((eopt - 2) > 0)
 210                                         {
 211                                                 error_1 = uerror(np, nr, truval, &alg);
 212                                         }
 213                                         else
 214                                         {
 215                                                 error_1 = errfac(truval, Fourie.rinc, rinc2);
 216                                         }
 217                                 }
 218                                 val *= pix4 / error_1;
 219
 220                                 // Y = SA*M#*A*(M*X - P)
 221
 222                                 for (nb = 0; nb < numb; nb++)
 223                                 {
 224
 225                                         k = list[nb];
 226                                         y[k] += val;
 227                                 }
 228
 229                                 // NEXT NR
 230                         }
 231
 232                 }
 233
 234                 // NEXT NP
 235
 236         }
 237         return;
 238 }