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) $
8 ***********************************************************
10 *********************************************************************************
11 Non-ascending vector for secondary optimization criterion "smoothness" (sec_cri2)
12 *********************************************************************************
24 void nav2_class::Run(REAL* xkn, REAL* vkn)
26 int pixelcount = GeoPar.area;
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;
36 void nav2_class::calcGradient(REAL* x, REAL* grad)
38 // read image dimensions
39 int nelem = GeoPar.nelem;
40 int pixelcount = GeoPar.area;
43 // define minimum value for divisions
44 double min = pow(10,-15);
46 // set gradient starting value
47 memset(grad,0,pixelcount*sizeof(REAL));
50 REAL* pre = new REAL[pixelcount];
51 memset(pre,0,pixelcount*sizeof(REAL));
53 // iterate through all pixels that are not on the border
54 for (int i=nelem+1; i<pixelcount-nelem-1; i++)
56 // skip left and right column
57 if (!(i%nelem) || !((i+1)%nelem)) continue;
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]));
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;
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;
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;
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;
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;
94 // the remainder of the top and bottom rows
96 for (int i=2; i<nelem-2; i++)
98 temp=-0.25*(pre[i+nelem-1]+pre[i+nelem]+pre[i+nelem+1]);
99 if (fabs(temp)>min) grad[i]=temp;
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;
106 // the remainder of the left and right columns
107 for (int i=2*nelem; i<nelem*(nelem-3)+1; i+=nelem)
109 temp=-0.25*(pre[i-nelem+1]+pre[i+1]+pre[i+nelem+1]);
110 if (fabs(temp)>min) grad[i]=temp;
113 temp=-0.25*(pre[t-nelem-1]+pre[t-1]+pre[t+nelem-1]);
114 if (fabs(temp)>min) grad[t]=temp;
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++)
122 // skip left and right column
123 if (!((i)%nelem) || !((i+1)%nelem)) continue;
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]);
127 if (fabs(temp)>min) grad[i]=temp;
133 for (int i=0;i<pixelcount;i++)