2 ***********************************************************
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) $
8 ***********************************************************
10 **************************************************************************************
11 Non-ascending vector for secondary optimization criterion "total variation" (sec_cri1)
12 **************************************************************************************
23 void nav1_class::Run(REAL* xkn, REAL* vkn)
25 int pixelcount = GeoPar.area;
27 REAL gradientSquareNorm = 0;
28 calcGradient(xkn, vkn);
29 for(INTEGER i = 0; i < pixelcount; i++)
30 gradientSquareNorm += pow(vkn[i], 2);
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;
37 for(INTEGER i = 0; i < pixelcount; i++)
41 void nav1_class::calcGradient(REAL* x, REAL* grad)
43 // read image dimensions
44 int nelem = GeoPar.nelem;
45 int pixelcount = GeoPar.area;
47 // define minimum value for divisions
48 const double min = pow(10,-15);
50 // set gradient starting value
51 memset(grad,0,pixelcount*sizeof(REAL));
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))
57 // (x[i]-x[i-1])/sqrt(pow(x[i-1]-x[i],2)+pow(x[i-1]-x[i-1+nelem],2))
60 // (x[i]-x[i-nelem])/sqrt(pow(x[i-nelem]-x[i+1-nelem],2)+pow(x[i-nelem]-x[i],2))
62 __________________________
64 |-|---------------------|-|
74 |_|_____________________|_|
76 -------------------------
78 // top left pixel: only formula 1 applies
79 double denx = pow(x[0]-x[1],2)+pow(x[0]-x[nelem],2);
81 grad[0]=(2*x[0]-x[1]-x[nelem])/sqrt(denx);
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);
91 grad[i] = ((x[i]-x[i-1])/sqrt(deny));
94 denx = (x[i]-x[i+1])*(x[i]-x[i+1]) + (x[i]-x[i+nelem])*(x[i]-x[i+nelem]);
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));
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));
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);
115 grad[i]=(x[i]-x[i-1])/sqrt(deny);
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);
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
135 grad[i] = ((x[i]-x[i-nelem])/sqrt(denz));
140 for(int i=0;i<(pixelcount-1);i++){