src/snark/quad_deset.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_deset.cpp $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  SET UP D MATRIX
  11
  12  COPT = 1 GOITEIN
  13
  14  COPT = 2 MODIFIED SIRT
  15
  16  ZERO OUT THE D VECTOR
  17  */
  18
  19 #include <cstdio>
  20 #include <cmath>
  21
  22 #include "blkdta.h"
  23 #include "geom.h"
  24 #include "consts.h"
  25 #include "uiod.h"
  26 #include "ray.h"
  27 #include "wray.h"
  28 #include "anglst.h"
  29 #include "errfac.h"
  30 #include "blob.h"
  31
  32 #include "quad.h"
  33
  34 void quad_class::deset(REAL* d, INTEGER* list, REAL* weight, BOOLEAN* alg)
  35 {
  36         INTEGER i;
  37         INTEGER np;
  38         INTEGER nr;
  39         INTEGER numb;
  40         REAL snorm;
  41         REAL totwei;
  42         INTEGER k;
  43         REAL truval;
  44         INTEGER j;
  45
  46         INTEGER area;
  47
  48         if (Blob.pix_basis)
  49         {
  50                 area = GeoPar.area;
  51         }
  52         else
  53         {
  54                 area = Blob.area;
  55         }
  56
  57         for (i = 0; i < area; i++)
  58         {
  59                 d[i] = 0.0;
  60         }
  61
  62         for (np = 0; np < GeoPar.prjnum; np++)
  63         {
  64                 for (nr = 0; nr < GeoPar.nrays; nr++)
  65                 {
  66                         if (GeoPar.line)
  67                         {
  68                                 if (Blob.pix_basis)
  69                                         wray(np, nr, list, weight, &numb, &snorm);
  70                                 else
  71                                         Blob.bwray(np, nr, list, weight, &numb, &snorm);
  72                         }
  73                         if (GeoPar.strip)
  74                         {
  75                                 if (Blob.pix_basis)
  76                                         ray(np, nr, list, weight, &numb, &snorm);
  77                                 else
  78                                         Blob.bray(np, nr, list, weight, &numb, &snorm);
  79                         }
  80                         if (numb != 0)
  81                         {
  82                                 if (copt != 1)
  83                                 {
  84
  85                                         // MODIFIED SIRT
  86
  87                                         totwei = 0.0;
  88
  89                                         for (k = 0; k < numb; k++)
  90                                         {
  91                                                 totwei += weight[k];
  92                                         }
  93
  94                                         if ((eopt - 2) >= 0)
  95                                         {
  96                                                 if ((eopt - 2) > 0)
  97                                                 {
  98                                                         truval = Anglst.prdta(np, nr);
  99                                                         error_1 = uerror(np, nr, truval, alg);
 100                                                 }
 101                                                 else
 102                                                 {
 103                                                         truval = Anglst.prdta(np, nr);
 104                                                         error_1 = errfac(truval, 0.0, 0.0);
 105                                                 }
 106                                         }
 107                                         totwei /= error_1;
 108
 109                                         for (k = 0; k < numb; k++)
 110                                         {
 111                                                 j = list[k];
 112                                                 d[j] += totwei * weight[k];
 113                                         }
 114                                 }
 115                                 else
 116                                 {
 117                                         // GOITEIN
 118
 119                                         if ((eopt - 2) >= 0)
 120                                         {
 121                                                 if ((eopt - 2) > 0)
 122                                                 {
 123                                                         truval = Anglst.prdta(np, nr);
 124                                                         error_1 = uerror(np, nr, truval, alg);
 125                                                 }
 126                                                 else
 127                                                 {
 128                                                         truval = Anglst.prdta(np, nr);
 129                                                         error_1 = errfac(truval, 0.0, 0.0);
 130                                                 }
 131                                         }
 132
 133                                         for (k = 0; k < numb; k++)
 134                                         {
 135                                                 j = list[k];
 136                                                 d[j] += (weight[k] * weight[k]) / error_1;
 137                                         }
 138                                 }
 139
 140                                 // NEXT NR
 141                         }
 142                 }
 143
 144 //     NEXT NP
 145
 146         }
 147
 148         for (i = 0; i < area; i++)
 149         {
 150                 if (fabs(d[i]) <= Consts.zero)
 151                 {
 152
 153                         if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
 154                         {
 155                                 *alg = TRUE;
 156                                 fprintf(output, "\n\n      ***a picture element is missed by all rays, can not compute matrix D, program stops***"); //hstau
 157                                 return;
 158                         }
 159                 }
 160
 161                 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
 162                 {
 163                         d[i] = (REAL) 1.0 / d[i];
 164                 }
 165         }
 166         return;
 167 }