src/snark/nav1.cpp

   1 /*
   2  ***********************************************************
   3  $SNARK_Header$
   4  $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/nav1.cpp $
   5  $LastChangedRevision: 153 $
   6  $Date: 2014-07-27 21:27:06 -0500 (Sun, 27 Jul 2014) $
   7  $Author: olangthaler $
   8  ***********************************************************
   9
  10 **************************************************************************************
  11 Non-ascending vector for secondary optimization criterion "total variation" (sec_cri1)
  12 **************************************************************************************
  13 */
  14
  15 #include <stdio.h>
  16 #include <string.h>
  17 #include <iostream>
  18 #include "nav1.h"
  19 #include "math.h"
  20 #include "anglst.h"
  21 #include "blob.h"
  22
  23 void nav1_class::Run(REAL* xkn, REAL* vkn)
  24 {
  25  int pixelcount = GeoPar.area;
  26
  27  REAL gradientSquareNorm = 0;
  28  calcGradient(xkn, vkn);
  29  for(INTEGER i = 0; i < pixelcount; i++)
  30      gradientSquareNorm += pow(vkn[i], 2);
  31
  32  REAL gradientNorm = sqrt(gradientSquareNorm);
  33  if(gradientNorm > pow(10,-15))
  34     for(INTEGER i = 0; i < pixelcount; i++)
  35         vkn[i] = (-1.0) * vkn[i] / gradientNorm;
  36  else
  37     for(INTEGER i = 0; i < pixelcount; i++)
  38         vkn[i] = 0.0;
  39 }
  40
  41 void nav1_class::calcGradient(REAL* x, REAL* grad)
  42 {
  43  // read image dimensions
  44  int nelem = GeoPar.nelem;
  45  int pixelcount = GeoPar.area;
  46
  47  // define minimum value for divisions
  48  const double min = pow(10,-15);
  49
  50  // set gradient starting value
  51  memset(grad,0,pixelcount*sizeof(REAL));
  52
  53  // formula 1
  54  // (2*x[i]-x[i+1]-x[i+nelem])/sqrt(pow(x[i]-x[i+1],2)+pow(x[i]-x[i+nelem],2))
  55
  56  // formula 2
  57  // (x[i]-x[i-1])/sqrt(pow(x[i-1]-x[i],2)+pow(x[i-1]-x[i-1+nelem],2))
  58
  59  // formula 3
  60  // (x[i]-x[i-nelem])/sqrt(pow(x[i-nelem]-x[i+1-nelem],2)+pow(x[i-nelem]-x[i],2))
  61  /*
  62   __________________________
  63  |1|       1 & 2         |2|
  64  |-|---------------------|-|
  65  | |                     | |
  66  | |                     | |
  67  |1|                     | |
  68  | |                     | |
  69  |&|      1 & 2 & 3      |2|
  70  | |                     | |
  71  |3|                     | |
  72  | |                     | |
  73  | |                     | |
  74  |_|_____________________|_|
  75  |3|          3          |0|
  76   -------------------------
  77  */
  78  // top left pixel: only formula 1 applies
  79  double denx = pow(x[0]-x[1],2)+pow(x[0]-x[nelem],2);
  80  if(denx>min){
  81     grad[0]=(2*x[0]-x[1]-x[nelem])/sqrt(denx);
  82    }
  83
  84  double deny;
  85  // (rest of) top row: formulas 1 and 2 apply
  86  for(int i=1; i<nelem; i++){
  87      deny = pow(x[i-1]-x[i],2)+pow(x[i-1]-x[i-1+nelem],2);
  88
  89      if((i+1)==nelem){
  90         if(deny>min)
  91            grad[i] = ((x[i]-x[i-1])/sqrt(deny));
  92        }
  93      else{
  94         denx = (x[i]-x[i+1])*(x[i]-x[i+1]) + (x[i]-x[i+nelem])*(x[i]-x[i+nelem]);
  95
  96         if(denx>min && deny>min)
  97            grad[i]=((2*x[i]-x[i+1]-x[i+nelem])/sqrt(denx)) + ((x[i]-x[i-1])/sqrt(deny));
  98        }
  99  }
 100
 101  double denz;
 102  // rest of image except left and right column and bottom row: all 3 formulas apply
 103  for(int i=nelem; i<pixelcount-nelem; i++){
 104      if(!(i%nelem)){ // first column: only cases 1 & 3
 105         denx = pow(x[i]-x[i+1],2) + pow(x[i]-x[i+nelem],2);
 106         denz = pow(x[i-nelem]-x[i-nelem+1],2) + pow(x[i-nelem]-x[i],2);
 107         if(denx>min && denz>min){
 108            grad[i]=((2*x[i]-x[i+1]-x[i+nelem])/sqrt(denx)) + ((x[i]-x[i-nelem])/sqrt(denz));
 109           }
 110        }
 111      else{
 112         if(!((i+1)%nelem)){ // last column: only case 2 applies
 113            deny = pow(x[i-1]-x[i],2) + pow(x[i-1]-x[i-1+nelem],2);
 114            if(deny>min){
 115               grad[i]=(x[i]-x[i-1])/sqrt(deny);
 116              }
 117           }
 118         else{
 119            denx = pow(x[i]-x[i+1],2)+pow(x[i]-x[i+nelem],2);
 120            deny = pow(x[i-1]-x[i],2)+pow(x[i-1]-x[i-1+nelem],2);
 121            denz = pow(x[i-nelem]-x[i-nelem+1],2)+pow(x[i-nelem]-x[i],2);
 122            if(denx>min && deny>min && denz>min){
 123               grad[i]=(2*x[i]-x[i+1]-x[i+nelem])/sqrt(denx) +
 124                       (x[i]-x[i-1])/sqrt(deny) +
 125                       (x[i]-x[i-nelem])/sqrt(denz);
 126              }
 127           }
 128        }
 129     }
 130  // Bottom row without last pixel: only third case applies
 131  for(int i=pixelcount-nelem; i<(pixelcount-1); i++){
 132      denz = pow(x[i-nelem]-x[i-nelem+1],2)+pow(x[i-nelem]-x[i],2);
 133      // left column: formulas 1 and 3 apply
 134      if(denz>min){
 135         grad[i] = ((x[i]-x[i-nelem])/sqrt(denz));
 136        }
 137     }
 138
 139  // failsafe
 140  for(int i=0;i<(pixelcount-1);i++){
 141      if(isnan(grad[i]))
 142         grad[i]=0.0;
 143     }
 144 }