src/snark/quad_matrix.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_matrix.cpp $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10
  11  THIS ROUTINE CALCULATES
  12
  13  Y = D*M#*A*M*X     IF COPT .LE. 3
  14
  15  OR
  16
  17  Y = C**P*(C*(SA+M#*A*M + SB*B**C) + SC)*X   OTHERWISE
  18
  19  WHERE P = T - 1
  20  */
  21
  22 #include <cstdio>
  23
  24 #include "blkdta.h"
  25 #include "geom.h"
  26 #include "consts.h"
  27 #include "blob.h"
  28
  29 #include "quad.h"
  30
  31 void quad_class::matrix(REAL* y, REAL* x, REAL* d, INTEGER* list, REAL* weight,
  32                 REAL* s)
  33 {
  34         REAL* j;
  35         INTEGER i;
  36         INTEGER area;
  37
  38         if (Blob.pix_basis)
  39         {
  40                 area = GeoPar.area;
  41         }
  42         else
  43         {
  44                 area = Blob.area;
  45         }
  46
  47         // COMPUTE SA*M#*A*M*X STORE IN Y
  48
  49         if (sa > Consts.zero)
  50                 mtamx(x, y, list, weight);
  51
  52         if (copt > 3)
  53         {
  54
  55                 if (sb > Consts.zero)
  56                 {
  57
  58                         j = new REAL[area];
  59
  60                         for (i = 0; i < area; i++)
  61                         {
  62                                 j[i] = x[i];
  63                         }
  64
  65                         if (l > 0)
  66                         {
  67                                 if (Blob.pix_basis)
  68                                 {
  69                                         adsmos(j, l, bw1, bw2, bw3, bbcon);
  70                                 }
  71                                 else
  72                                 {
  73                                         badsmos(j, l, bw1, bw2, bbcon);
  74                                 }
  75                         }
  76                         for (i = 0; i < area; i++)
  77                         {
  78                                 y[i] += sb * j[i];
  79                         }
  80                         delete[] j;  // bug 92 - Lajos - 03/02/2005
  81                 }
  82                 if (copt != 5)
  83                 {
  84                         if (Blob.pix_basis)
  85                         {
  86                                 adsmos(y, 1, cw1, cw2, cw3, cbcon);
  87                         }
  88                         else
  89                         {
  90                                 badsmos(y, 1, cw1, cw2, cbcon);
  91                         }
  92                 }
  93
  94                 if (sc != 0.0)
  95                 {
  96
  97                         for (i = 0; i < area; i++)
  98                         {
  99                                 y[i] += sc * x[i];
 100                         }
 101                 }
 102
 103                 if ((copt != 5) && ((t - 1) > 0))
 104                 {
 105                         if (Blob.pix_basis)
 106                         {
 107                                 adsmos(y, t - 1, cw1, cw2, cw3, cbcon);
 108                         }
 109                         else
 110                         {
 111                                 badsmos(y, t - 1, cw1, cw2, cbcon);
 112                         }
 113                 }
 114
 115                 if (gopt > 0)
 116                 {
 117                         for (i = 0; i < area; i++)
 118                         {
 119                                 y[i] *= s[i];
 120                         }
 121                 }
 122
 123                 if ((copt != 5) && (t > 0))
 124                 {
 125                         if (Blob.pix_basis)
 126                         {
 127                                 adsmos(y, t, cw1, cw2, cw3, cbcon);
 128                         }
 129                         else
 130                         {
 131                                 badsmos(y, t, cw1, cw2, cbcon);
 132                         }
 133                 }
 134                 return;
 135         }
 136         // D SECTION
 137
 138         for (i = 0; i < area; i++)
 139         {
 140                 y[i] *= d[i];
 141         }
 142
 143         return;
 144 }