summaryrefslogtreecommitdiffstats
path: root/SD-VBS/benchmarks/sift/src/c/siftrefinemx.c
diff options
context:
space:
mode:
Diffstat (limited to 'SD-VBS/benchmarks/sift/src/c/siftrefinemx.c')
-rw-r--r--SD-VBS/benchmarks/sift/src/c/siftrefinemx.c196
1 files changed, 196 insertions, 0 deletions
<
diff --git a/SD-VBS/benchmarks/sift/src/c/siftrefinemx.c b/SD-VBS/benchmarks/sift/src/c/siftrefinemx.c
new file mode 100644
index 0000000..2b10e9d
--- /dev/null
+++ b/SD-VBS/benchmarks/sift/src/c/siftrefinemx.c
@@ -0,0 +1,196 @@
1#include "sift.h"
2
3const int max_iter = 5 ;
4
5/**
6
7 SIFTREFINEMX Subpixel localization, thresholding and on-edge test
8 Q = SIFTREFINEMX(P, OCTAVE, SMIN) refines, thresholds and performs
9 the on-edge test for the SIFT frames P extracted from the DOG
10 octave OCTAVE with parameter SMIN (see GAUSSIANSS()).
11
12 Q = SIFTREFINEMX(P, OCTAVE, SMIN, THRESH, R) specifies custom
13 values for the local maximum strength threshold THRESH and the
14 local maximum peakedeness threshold R.
15
16 OCTAVE is an octave of the Difference Of Gaussian scale space. P
17 is a 3xK matrix specifying the indexes (X,Y,S) of the points of
18 extremum of the octave OCTAVE. The spatial indexes X,Y are integer
19 with base zero. The scale index S is integer with base SMIN and
20 represents a scale sublevel in the specified octave.
21
22 The function returns a matrix Q containing the refined keypoints.
23 The matrix has the same format as P, except that the indexes are
24 now fractional. The function drops the points that do not satisfy
25 the strength and peakedness tests.
26
27**/
28
29F2D* siftrefinemx(F2D* oframes, F2D* dogss, int smin, float thresh, int r, int M, int N, int intervals)
30{
31 int S,K ;
32 F2D* out;
33
34 /* -----------------------------------------------------------------
35 ** Check the arguments
36 ** -------------------------------------------------------------- */
37
38 S = intervals;
39 K = oframes->height;
40
41 /* If the input array is empty, then output an empty array as well. */
42 if( K == 0) {
43 out = fDeepCopy(oframes);
44 return out;
45 }
46
47 /* -----------------------------------------------------------------
48 * Do the job
49 * -------------------------------------------------------------- */
50 {
51 F2D* buffer = fMallocHandle(K,3);
52 int p ;
53 const int yo = 1 ;
54 const int xo = M ;
55 const int so = M*N;
56 int buffCounter = 0;
57 int pptcount = 0;
58
59 for(p = 0 ; p < K ; ++p)
60 {
61 float tx, ty, ts;
62 int x,y,s;
63 int iter ;
64 float b[3] ;
65
66 tx = asubsref(oframes, pptcount++);
67 ty = asubsref(oframes, pptcount++);
68 ts = asubsref(oframes, pptcount++);
69
70 x = ceil(tx);
71 y = ceil(ty);
72 s = ceil(ts) - smin;
73
74
75 /* Local maxima extracted from the DOG
76 * have coorrinates 1<=x<=N-2, 1<=y<=M-2
77 * and 1<=s-mins<=S-2. This is also the range of the points
78 * that we can refine.
79 */
80
81 if(x < 1 || x > N-2 ||
82 y < 1 || y > M-2 ||
83 s < 1 || s > S-2) {
84 continue ;
85 }
86
87#define at(dx,dy,ds) asubsref(dogss, (dx)*xo + (dy)*yo + (ds)*so)
88
89 {
90 float Dx=0,Dy=0,Ds=0,Dxx=0,Dyy=0,Dss=0,Dxy=0,Dxs=0,Dys=0 ;
91 int dx = 0 ;
92 int dy = 0 ;
93
94 for(iter = 0 ; iter < max_iter ; ++iter)
95 {
96
97 float A[3*3] ;
98
99#define Aat(i,j) (A[(i)+(j)*3])
100
101 x += dx ;
102 y += dy ;
103
104 /* Compute the gradient. */
105 Dx = 0.5 * (at(x+1,y+0,s+0) - at(x-1,y+0,s+0)) ;
106 Dy = 0.5 * (at(x+0,y+1,s+0) - at(x+0,y-1,s+0));
107 Ds = 0.5 * (at(x+0,y+0,s+1) - at(x+0,y+0,s-1)) ;
108
109 /* Compute the Hessian. */
110 Dxx = (at(x+1,y,s) + at(x-1,y,s) - 2.0 * at(x,y,s)) ;
111 Dyy = (at(x,y+1,s) + at(x,y-1,s) - 2.0 * at(x,y,s)) ;
112 Dss = (at(x,y,s+1) + at(x,y,s-1) - 2.0 * at(x,y,s)) ;
113
114 Dxy = 0.25 * ( at(x+1,y+1,s) + at(x-1,y-1,s) - at(x-1,y+1,s) - at(x+1,y-1,s) ) ;
115 Dxs = 0.25 * ( at(x+1,y,s+1) + at(x-1,y,s-1) - at(x-1,y,s+1) - at(x+1,y,s-1) ) ;
116 Dys = 0.25 * ( at(x,y+1,s+1) + at(x,y-1,s-1) - at(x,y-1,s+1) - at(x,y+1,s-1) ) ;
117
118 /* Solve linear system. */
119 Aat(0,0) = Dxx ;
120 Aat(1,1) = Dyy ;
121 Aat(2,2) = Dss ;
122 Aat(0,1) = Aat(1,0) = Dxy ;
123 Aat(0,2) = Aat(2,0) = Dxs ;
124 Aat(1,2) = Aat(2,1) = Dys ;
125
126 b[0] = - Dx ;
127 b[1] = - Dy ;
128 b[2] = - Ds ;
129
130 /* If the translation of the keypoint is big, move the keypoint
131 * and re-iterate the computation. Otherwise we are all set.
132 */
133 dx= ((b[0] > 0.6 && x < N-2) ? 1 : 0 )
134 + ((b[0] < -0.6 && x > 1 ) ? -1 : 0 ) ;
135
136 dy= ((b[1] > 0.6 && y < M-2) ? 1 : 0 )
137 + ((b[1] < -0.6 && y > 1 ) ? -1 : 0 ) ;
138
139 if( dx == 0 && dy == 0 )
140 break ;
141
142 }
143
144 {
145 float val = at(x,y,s) + 0.5 * (Dx * b[0] + Dy * b[1] + Ds * b[2]) ;
146 float score = (Dxx+Dyy)*(Dxx+Dyy) / (Dxx*Dyy - Dxy*Dxy) ;
147 float xn = x + b[0] ;
148 float yn = y + b[1] ;
149 float sn = s + b[2] ;
150
151 if(fabs(val) > thresh &&
152 score < (r+1)*(r+1)/r &&
153 score >= 0 &&
154 fabs(b[0]) < 1.5 &&
155 fabs(b[1]) < 1.5 &&
156 fabs(b[2]) < 1.5 &&
157 xn >= 0 &&
158 xn <= N-1 &&
159 yn >= 0 &&
160 yn <= M-1 &&
161 sn >= 0 &&
162 sn <= S-1)
163 {
164 asubsref(buffer,buffCounter++) = xn ;
165 asubsref(buffer,buffCounter++) = yn ;
166 asubsref(buffer,buffCounter++) = sn+smin ;
167 }
168 }
169 }
170 }
171
172 /* Copy the result into an array. */
173 {
174 int i, j, k=0;
175 int NL = buffCounter/3;
176 out = fMallocHandle(3, NL);
177
178 for(i=0; i<NL; i++)
179 {
180 for(j=0; j<3; j++)
181 {
182 subsref(out,j,i) = asubsref(buffer,k);
183 k++;
184 }
185 }
186 }
187 free(buffer) ;
188 }
189
190 return out;
191}
192
193