summaryrefslogtreecommitdiffstats
path: root/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.c
diff options
context:
space:
mode:
Diffstat (limited to 'SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.c')
-rw-r--r--SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.c342
1 files changed, 0 insertions, 342 deletions
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.c b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.c
deleted file mode 100644
index a1ef8c7..0000000
--- a/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.c
+++ /dev/null
@@ -1,342 +0,0 @@
1/* file: siftrefinemx.c
2
3Benchmark - sift
4Data set - sqcif
5
6 < M A T L A B >
7 Copyright 1984-2006 The MathWorks, Inc.
8 Version 7.3.0.298 (R2006b)
9 August 03, 2006
10
11
12 To get started, type one of these: helpwin, helpdesk, or demo.
13 For product information, visit www.mathworks.com.
14
15Warning: Function /h/g2/kvs/checkParallel/sdvbs-svn/common/matlab/randn.m has the same name as a MATLAB builtin. We suggest you rename the function to avoid a potential name conflict.
16> In path at 113
17 In script_run_profile at 3
18Warning: You are using gcc version "4.1.1". The earliest gcc version supported
19with mex is "3.4.0". The latest version tested for use with mex is "3.4.5".
20To download a different version of gcc, visit http://gcc.gnu.org
21Warning: You are using gcc version "4.1.1". The earliest gcc version supported
22with mex is "3.4.0". The latest version tested for use with mex is "3.4.5".
23To download a different version of gcc, visit http://gcc.gnu.org
24Warning: You are using gcc version "4.1.1". The earliest gcc version supported
25with mex is "3.4.0". The latest version tested for use with mex is "3.4.5".
26To download a different version of gcc, visit http://gcc.gnu.org
27Warning: You are using gcc version "4.1.1". The earliest gcc version supported
28with mex is "3.4.0". The latest version tested for use with mex is "3.4.5".
29To download a different version of gcc, visit http://gcc.gnu.org
30Warning: You are using gcc version "4.1.1". The earliest gcc version supported
31with mex is "3.4.0". The latest version tested for use with mex is "3.4.5".
32To download a different version of gcc, visit http://gcc.gnu.org
33Warning: You are using gcc version "4.1.1". The earliest gcc version supported
34with mex is "3.4.0". The latest version tested for use with mex is "3.4.5".
35To download a different version of gcc, visit http://gcc.gnu.org
36Input size - (96x96)
37** author: Andrea Vedaldi
38** description: Subpixel localization, thresholding and off-edge test.
39**/
40
41/* AUTORIGHTS
42Copyright (c) 2006 The Regents of the University of California.
43All Rights Reserved.
44
45Created by Andrea Vedaldi
46UCLA Vision Lab - Department of Computer Science
47
48Permission to use, copy, modify, and distribute this software and its
49documentation for educational, research and non-profit purposes,
50without fee, and without a written agreement is hereby granted,
51provided that the above copyright notice, this paragraph and the
52following three paragraphs appear in all copies.
53
54This software program and documentation are copyrighted by The Regents
55of the University of California. The software program and
56documentation are supplied "as is", without any accompanying services
57from The Regents. The Regents does not warrant that the operation of
58the program will be uninterrupted or error-free. The end-user
59understands that the program was developed for research purposes and
60is advised not to rely exclusively on the program for any reason.
61
62This software embodies a method for which the following patent has
63been issued: "Method and apparatus for identifying scale invariant
64features in an image and use of same for locating an object in an
65image," David G. Lowe, US Patent 6,711,293 (March 23,
662004). Provisional application filed March 8, 1999. Asignee: The
67University of British Columbia.
68
69IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
70FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
71INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
72ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
73ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
74CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
75LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
76A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
77BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
78MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
79*/
80
81#include"mex.h"
82
83#include<mexutils.c>
84
85#include<stdlib.h>
86#include<string.h>
87
88/* Prototype of DGESV LAPACK function for the solution of a linear system. */
89#ifdef WINDOWS
90#define DGESV dgesv
91#undef min
92#undef max
93#else
94#define DGESV dgesv_
95#endif
96
97#ifdef WINDOWS
98#ifdef __cplusplus__
99extern "C" {
100 extern int DGESV(int *n, int *nrhs, double *a, int *lda,
101 int *ipiv, double *b, int *ldb, int *info) ;
102}
103#else
104 extern int DGESV(int *n, int *nrhs, double *a, int *lda,
105 int *ipiv, double *b, int *ldb, int *info) ;
106#define sqrtf(x) ((float)sqrt((double)x)
107#define powf(x) ((float)pow((double)x)
108#define fabsf(x) ((float)fabs((double)x)
109#endif
110#else
111extern int DGESV(int *n, int *nrhs, double *a, int *lda,
112 int *ipiv, double *b, int *ldb, int *info) ;
113#endif
114
115#define greater(a,b) ((a) > (b))
116#define min(a,b) (((a)<(b))?(a):(b))
117#define max(a,b) (((a)>(b))?(a):(b))
118
119const int max_iter = 5 ;
120
121void
122mexFunction(int nout, mxArray *out[],
123 int nin, const mxArray *in[])
124{
125 int M,N,S,smin,K ;
126 const int* dimensions ;
127 const double* P_pt ;
128 const double* D_pt ;
129 double threshold = 0.01 ; /*0.02 ;*/
130 double r = 10.0 ;
131 double* result ;
132 enum {IN_P=0,IN_D,IN_SMIN,IN_THRESHOLD,IN_R} ;
133 enum {OUT_Q=0} ;
134
135 /* -----------------------------------------------------------------
136 ** Check the arguments
137 ** -------------------------------------------------------------- */
138 if (nin < 3) {
139 mexErrMsgTxt("At least three input arguments required.");
140 } else if (nout > 1) {
141 mexErrMsgTxt("Too many output arguments.");
142 }
143
144 if( !uIsRealMatrix(in[IN_P],3,-1) ) {
145 mexErrMsgTxt("P must be a 3xK real matrix") ;
146 }
147
148 if( !mxIsDouble(in[IN_D]) || mxGetNumberOfDimensions(in[IN_D]) != 3) {
149 mexErrMsgTxt("G must be a three dimensional real array.") ;
150 }
151
152 if( !uIsRealScalar(in[IN_SMIN]) ) {
153 mexErrMsgTxt("SMIN must be a real scalar.") ;
154 }
155
156 if(nin >= 4) {
157 if(!uIsRealScalar(in[IN_THRESHOLD])) {
158 mexErrMsgTxt("THRESHOLD must be a real scalar.") ;
159 }
160 threshold = *mxGetPr(in[IN_THRESHOLD]) ;
161 }
162
163 if(nin >= 5) {
164 if(!uIsRealScalar(in[IN_R])) {
165 mexErrMsgTxt("R must be a real scalar.") ;
166 }
167 r = *mxGetPr(in[IN_R]) ;
168 }
169
170 dimensions = mxGetDimensions(in[IN_D]) ;
171 M = dimensions[0] ;
172 N = dimensions[1] ;
173 S = dimensions[2] ;
174 smin = (int)(*mxGetPr(in[IN_SMIN])) ;
175
176 if(S < 3 || M < 3 || N < 3) {
177 mexErrMsgTxt("All dimensions of DOG must be not less than 3.") ;
178 }
179
180 K = mxGetN(in[IN_P]) ;
181 P_pt = mxGetPr(in[IN_P]) ;
182 D_pt = mxGetPr(in[IN_D]) ;
183
184 /* If the input array is empty, then output an empty array as well. */
185 if( K == 0) {
186 out[OUT_Q] = mxDuplicateArray(in[IN_P]) ;
187 return ;
188 }
189
190 /* -----------------------------------------------------------------
191 * Do the job
192 * -------------------------------------------------------------- */
193 {
194 double* buffer = (double*) mxMalloc(K*3*sizeof(double)) ;
195 double* buffer_iterator = buffer ;
196 int p ;
197 const int yo = 1 ;
198 const int xo = M ;
199 const int so = M*N ;
200
201/* printf("Actual values = %d\n\n", K);
202*/
203 for(p = 0 ; p < K ; ++p) {
204 int x = ((int)*P_pt++) ;
205 int y = ((int)*P_pt++) ;
206/* printf("%d\t%d\n", ((int)*P_pt), smin);
207*/
208 int s = ((int)*P_pt++) - smin ;
209 int iter ;
210 double b[3] ;
211
212 /* Local maxima extracted from the DOG
213 * have coorrinates 1<=x<=N-2, 1<=y<=M-2
214 * and 1<=s-mins<=S-2. This is also the range of the points
215 * that we can refine.
216 */
217/* printf("%d\t%d\t%d\t%d\t%d\t%d\n", x, N-2,y,M-2, s, S-2);
218*/
219 if(x < 1 || x > N-2 ||
220 y < 1 || y > M-2 ||
221 s < 1 || s > S-2) {
222 continue ;
223 }
224
225#define at(dx,dy,ds) (*(pt + (dx)*xo + (dy)*yo + (ds)*so))
226
227 {
228 const double* pt = D_pt + y*yo + x*xo + s*so ;
229
230 double Dx=0,Dy=0,Ds=0,Dxx=0,Dyy=0,Dss=0,Dxy=0,Dxs=0,Dys=0 ;
231 int dx = 0 ;
232 int dy = 0 ;
233/* printf("%d\t%d\t%d\t%d\t%d\t%d\t%d\t%f\t%f\n",S, y, yo, x, xo, s, so, *D_pt, *pt);
234*/
235 for(iter = 0 ; iter < max_iter ; ++iter) {
236
237 double A[3*3] ;
238 int ipiv[3] ;
239 int n = 3 ;
240 int one = 1 ;
241 int info = 0 ;
242
243#define Aat(i,j) (A[(i)+(j)*3])
244
245 x += dx ;
246 y += dy ;
247 pt = D_pt + y*yo + x*xo + s*so ;
248
249 /* Compute the gradient. */
250 Dx = 0.5 * (at(+1,0,0) - at(-1,0,0)) ;
251 Dy = 0.5 * (at(0,+1,0) - at(0,-1,0));
252 Ds = 0.5 * (at(0,0,+1) - at(0,0,-1)) ;
253
254 /* Compute the Hessian. */
255 Dxx = (at(+1,0,0) + at(-1,0,0) - 2.0 * at(0,0,0)) ;
256 Dyy = (at(0,+1,0) + at(0,-1,0) - 2.0 * at(0,0,0)) ;
257 Dss = (at(0,0,+1) + at(0,0,-1) - 2.0 * at(0,0,0)) ;
258
259 Dxy = 0.25 * ( at(+1,+1,0) + at(-1,-1,0) - at(-1,+1,0) - at(+1,-1,0) ) ;
260 Dxs = 0.25 * ( at(+1,0,+1) + at(-1,0,-1) - at(-1,0,+1) - at(+1,0,-1) ) ;
261 Dys = 0.25 * ( at(0,+1,+1) + at(0,-1,-1) - at(0,-1,+1) - at(0,+1,-1) ) ;
262
263 /* Solve linear system. */
264 Aat(0,0) = Dxx ;
265 Aat(1,1) = Dyy ;
266 Aat(2,2) = Dss ;
267 Aat(0,1) = Aat(1,0) = Dxy ;
268 Aat(0,2) = Aat(2,0) = Dxs ;
269 Aat(1,2) = Aat(2,1) = Dys ;
270
271 b[0] = - Dx ;
272 b[1] = - Dy ;
273 b[2] = - Ds ;
274
275 /* DGESV (&n, &one, A, &n, ipiv, b, &n, &info) ;
276 */
277 /* If the translation of the keypoint is big, move the keypoint
278 * and re-iterate the computation. Otherwise we are all set.
279 */
280 dx= ((b[0] > 0.6 && x < N-2) ? 1 : 0 )
281 + ((b[0] < -0.6 && x > 1 ) ? -1 : 0 ) ;
282
283 dy= ((b[1] > 0.6 && y < M-2) ? 1 : 0 )
284 + ((b[1] < -0.6 && y > 1 ) ? -1 : 0 ) ;
285
286 if( dx == 0 && dy == 0 ) break ;
287
288 }
289
290 {
291 double val = at(0,0,0) + 0.5 * (Dx * b[0] + Dy * b[1] + Ds * b[2]) ;
292 double score = (Dxx+Dyy)*(Dxx+Dyy) / (Dxx*Dyy - Dxy*Dxy) ;
293 double xn = x + b[0] ;
294 double yn = y + b[1] ;
295 double sn = s + b[2] ;
296
297/* printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", fabs(val),threshold,score,(r+1)*(r+1)/r,fabs(b[0]), fabs(b[1]), fabs(b[2]),xn,yn,sn,r);
298*/
299 if(fabs(val) > threshold &&
300 score < (r+1)*(r+1)/r &&
301 score >= 0 &&
302 fabs(b[0]) < 1.5 &&
303 fabs(b[1]) < 1.5 &&
304 fabs(b[2]) < 1.5 &&
305 xn >= 0 &&
306 xn <= N-1 &&
307 yn >= 0 &&
308 yn <= M-1 &&
309 sn >= 0 &&
310 sn <= S-1) {
311 *buffer_iterator++ = xn ;
312 *buffer_iterator++ = yn ;
313 *buffer_iterator++ = sn+smin ;
314 }
315 }
316 }
317 }
318
319 /* Copy the result into an array. */
320 {
321 int i;
322 int NL = (buffer_iterator - buffer)/3 ;
323/* printf("%NL VALUE = %d\t%d\t%d\n", NL, buffer_iterator, buffer);
324*/
325 out[OUT_Q] = mxCreateDoubleMatrix(3, NL, mxREAL) ;
326 result = mxGetPr(out[OUT_Q]);
327 for(i=0; i<(3*NL); i++)
328 {
329 result[i] = buffer[i];
330/* printf("%f\t", buffer[i]);
331*/
332 }
333/* printf("\n");
334 memcpy(result, buffer, sizeof(double) * 3 * NL) ;
335*/
336 }
337 mxFree(buffer) ;
338 }
339
340}
341
342
generated by cgit v1.2.2 (git 2.25.0) at 2025-04-30 17:55:46 -0400