src/snark/nav2.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/nav2.cpp $
   5 $LastChangedRevision: 145 $
   6 $Date: 2014-07-17 19:49:25 -0400 (Thu, 17 Jul 2014) $
   7 $Author: olangthaler $
   8 ***********************************************************
   9
  10 *********************************************************************************
  11 Non-ascending vector for secondary optimization criterion "smoothness" (sec_cri2)
  12 *********************************************************************************
  13 */
  14
  15 #include <cstdio>
  16 #include <stdlib.h>
  17 #include <string.h>
  18 #include <iostream>
  19 #include "nav2.h"
  20 #include "math.h"
  21 #include "anglst.h"
  22 #include "blob.h"
  23
  24 void nav2_class::Run(REAL* xkn, REAL* vkn)
  25 {
  26         int pixelcount = GeoPar.area;
  27
  28         REAL gradientSquareNorm = 0;
  29         calcGradient(xkn, vkn);
  30         for (INTEGER i = 0; i < pixelcount; i++) gradientSquareNorm += pow(vkn[i], 2);
  31         REAL gradientNorm = sqrt(gradientSquareNorm);
  32         if (gradientNorm > pow(10,-15)) for (INTEGER i = 0; i < pixelcount; i++) vkn[i] = (-1.0) * vkn[i] / gradientNorm;
  33         else for (INTEGER i = 0; i < pixelcount; i++) vkn[i] = 0;
  34 }
  35
  36 void nav2_class::calcGradient(REAL* x, REAL* grad)
  37 {
  38         // read image dimensions
  39         int nelem = GeoPar.nelem;
  40         int pixelcount = GeoPar.area;
  41         double temp=0;
  42
  43         // define minimum value for divisions
  44         double min = pow(10,-15);
  45
  46         // set gradient starting value
  47         memset(grad,0,pixelcount*sizeof(REAL));
  48
  49         // pre-processing
  50         REAL* pre = new REAL[pixelcount];
  51         memset(pre,0,pixelcount*sizeof(REAL));
  52
  53         // iterate through all pixels that are not on the border
  54         for (int i=nelem+1; i<pixelcount-nelem-1; i++)
  55         {
  56                 // skip left and right column
  57                 if (!(i%nelem) || !((i+1)%nelem)) continue;
  58
  59                 pre[i]=x[i]-(0.125*(x[i-nelem-1]+x[i-nelem]+x[i-nelem+1]+x[i-1]+x[i+1]+x[i+nelem-1]+x[i+nelem]+x[i+nelem+1]));
  60         }
  61
  62         // actual processing
  63         // the four corners
  64         temp=-0.25*pre[nelem+1];
  65         if (fabs(temp)>min) grad[0]=temp;
  66         temp=-0.25*pre[2*nelem-2];
  67         if (fabs(temp)>min) grad[nelem-1]=temp;
  68         temp=-0.25*pre[pixelcount-2*nelem+1];
  69         if (fabs(temp)>min) grad[pixelcount-nelem]=temp;
  70         temp=-0.25*pre[pixelcount-nelem-2];
  71         if (fabs(temp)>min) grad[pixelcount-1]=temp;
  72
  73         // the 8 pixels adjoining the corners
  74         temp=-0.25*(pre[nelem+1]+pre[nelem+2]);
  75         if (fabs(temp)>min) grad[1]=temp;
  76         temp=-0.25*(pre[nelem+1]+pre[2*nelem+1]);
  77         if (fabs(temp)>min) grad[nelem]=temp;
  78
  79         temp=-0.25*(pre[2*nelem-2]+pre[2*nelem-3]);
  80         if (fabs(temp)>min) grad[nelem-2]=temp;
  81         temp=-0.25*(pre[2*nelem-2]+pre[3*nelem-2]);
  82         if (fabs(temp)>min) grad[2*nelem-1]=temp;
  83
  84         temp=-0.25*(pre[pixelcount-2*nelem+1]+pre[pixelcount-2*nelem+2]);
  85         if (fabs(temp)>min) grad[pixelcount-nelem+1]=temp;
  86         temp=-0.25*(pre[pixelcount-2*nelem+1]+pre[pixelcount-3*nelem+1]);
  87         if (fabs(temp)>min) grad[pixelcount-2*nelem]=temp;
  88
  89         temp=-0.25*(pre[pixelcount-nelem-2]+pre[pixelcount-nelem-3]);
  90         if (fabs(temp)>min) grad[pixelcount-2]=temp;
  91         temp=-0.25*(pre[pixelcount-nelem-2]+pre[pixelcount-2*nelem-2]);
  92         if (fabs(temp)>min) grad[pixelcount-nelem-1]=temp;
  93
  94         // the remainder of the top and bottom rows
  95         int t=0;
  96         for (int i=2; i<nelem-2; i++)
  97         {
  98                 temp=-0.25*(pre[i+nelem-1]+pre[i+nelem]+pre[i+nelem+1]);
  99                 if (fabs(temp)>min) grad[i]=temp;
 100
 101                 t = i+(nelem*(nelem-1));
 102                 temp=-0.25*(pre[t-nelem-1]+pre[t-nelem]+pre[t-nelem+1]);
 103                 if (fabs(temp)>min) grad[t]=temp;
 104         }
 105
 106         // the remainder of the left and right columns
 107         for (int i=2*nelem; i<nelem*(nelem-3)+1; i+=nelem)
 108         {
 109                 temp=-0.25*(pre[i-nelem+1]+pre[i+1]+pre[i+nelem+1]);
 110                 if (fabs(temp)>min) grad[i]=temp;
 111
 112                 t = i+nelem-1;
 113                 temp=-0.25*(pre[t-nelem-1]+pre[t-1]+pre[t+nelem-1]);
 114                 if (fabs(temp)>min) grad[t]=temp;
 115         }
 116
 117         // the remainder of the image: all pixels that are not on the border
 118         // note: since the border pixels of the preprocessed array are 0, the
 119         // pixels of the "inner border" do not need to be handled separately
 120         for (int i=nelem+1; i<pixelcount-nelem-1; i++)
 121         {
 122                 // skip left and right column
 123                 if (!((i)%nelem) || !((i+1)%nelem)) continue;
 124
 125                 temp=2*pre[i]-0.25*(pre[i-nelem-1]+pre[i-nelem]+pre[i-nelem+1]+pre[i-1]+pre[i+1]+pre[i+nelem-1]+pre[i+nelem]+pre[i+nelem+1]);
 126
 127                 if (fabs(temp)>min) grad[i]=temp;
 128         }
 129
 130         delete[] (pre);
 131
 132         // failsafe
 133         for (int i=0;i<pixelcount;i++)
 134         {
 135                 if (isnan(grad[i]))
 136                 {
 137                         grad[i]=0;
 138                 }
 139         }
 140 }