src/snark/emap_Cwray.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/emap_Cwray.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  */
  11
  12 #include <cstdio>
  13
  14 #include "blkdta.h"
  15 #include "wray.h"
  16 #include "blob.h"
  17 #include "consts.h"
  18 #include "uiod.h"
  19
  20 #include "emap.h"
  21
  22 BOOLEAN emap_class::findqpart(INTEGER np, INTEGER nr, REAL yi, REAL* recon,
  23                 REAL* qtilde)
  24 {
  25         //----- local variables
  26         INTEGER t, nt, bigT, n[4096];
  27         REAL l[4096], snorm, si;
  28
  29         //------ Body of outer loop in Algorithm FIND_q~ for the MAP algorithm * /
  30
  31         if (Blob.pix_basis)
  32         {
  33                 wray(np, nr, n, l, &bigT, &snorm);
  34         }
  35         else
  36         {
  37                 Blob.bwray(np, nr, n, l, &bigT, &snorm); // l and n must have enough size for blob basis. Why 4096? Hstau 6/25/2003
  38         }
  39
  40         if (bigT != 0)
  41         {
  42                 si = 0.0;
  43                 for (t = 0; t < bigT; t++)
  44                 {
  45                         si += recon[n[t]] * l[t];
  46                 }
  47
  48                 if (si < Consts.zero)
  49                 {
  50                         fprintf(output,
  51                                         "\n****The pseudo raysum for  np = %d, nr = %d is zero, EMAP is not applicable, EMAP aborted!\n",
  52                                         np, nr);
  53                         return TRUE;
  54                 }                                                  //added by wei 1/2005
  55
  56                 //cvd$     nodepchk
  57                 for (t = 0; t < bigT; t++)
  58                 {
  59                         nt = n[t];
  60                         qtilde[nt] += yi / si * l[t];
  61                 }
  62         }
  63
  64         return FALSE;
  65 }