summaryrefslogtreecommitdiffstats
path: root/SD-VBS/benchmarks/sift/src/matlab/siftormx.c
diff options
context:
space:
mode:
Diffstat (limited to 'SD-VBS/benchmarks/sift/src/matlab/siftormx.c')
-rw-r--r--SD-VBS/benchmarks/sift/src/matlab/siftormx.c251
1 files changed, 0 insertions, 251 deletions
diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftormx.c b/SD-VBS/benchmarks/sift/src/matlab/siftormx.c
deleted file mode 100644
index 40380a2..0000000
--- a/SD-VBS/benchmarks/sift/src/matlab/siftormx.c
+++ /dev/null
@@ -1,251 +0,0 @@
1/* file: siftormx.c
2** author: Andrea Vedaldi
3** description: Computes peaks of orientation histogram.
4**/
5
6/* AUTORIGHTS
7Copyright (c) 2006 The Regents of the University of California.
8All Rights Reserved.
9
10Created by Andrea Vedaldi
11UCLA Vision Lab - Department of Computer Science
12
13Permission to use, copy, modify, and distribute this software and its
14documentation for educational, research and non-profit purposes,
15without fee, and without a written agreement is hereby granted,
16provided that the above copyright notice, this paragraph and the
17following three paragraphs appear in all copies.
18
19This software program and documentation are copyrighted by The Regents
20of the University of California. The software program and
21documentation are supplied "as is", without any accompanying services
22from The Regents. The Regents does not warrant that the operation of
23the program will be uninterrupted or error-free. The end-user
24understands that the program was developed for research purposes and
25is advised not to rely exclusively on the program for any reason.
26
27This software embodies a method for which the following patent has
28been issued: "Method and apparatus for identifying scale invariant
29features in an image and use of same for locating an object in an
30image," David G. Lowe, US Patent 6,711,293 (March 23,
312004). Provisional application filed March 8, 1999. Asignee: The
32University of British Columbia.
33
34IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
35FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES,
36INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
37ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN
38ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF
39CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
40LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
41A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
42BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
43MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
44*/
45
46#include"mex.h"
47#include<stdlib.h>
48#include<string.h>
49#include<math.h>
50
51#include<mexutils.c>
52
53#define greater(a,b) a > b
54#define min(a,b) (((a)<(b))?(a):(b))
55#define max(a,b) (((a)>(b))?(a):(b))
56
57const double win_factor = 1.5 ;
58#define NBINS 36
59
60void
61mexFunction(int nout, mxArray *out[],
62 int nin, const mxArray *in[])
63{
64 int M,N,S,smin,K ;
65 const int* dimensions ;
66 const double* P_pt ;
67 const double* G_pt ;
68 double* TH_pt ;
69 double sigma0 ;
70 double H_pt [ NBINS ] ;
71
72 enum {IN_P=0,IN_G,IN_S,IN_SMIN,IN_SIGMA0} ;
73 enum {OUT_Q=0} ;
74
75 /* -----------------------------------------------------------------
76 ** Check the arguments
77 ** -------------------------------------------------------------- */
78 if (nin != 5) {
79 mexErrMsgTxt("Exactly five input arguments required.");
80 } else if (nout > 1) {
81 mexErrMsgTxt("Too many output arguments.");
82 }
83
84 if( !uIsRealScalar(in[IN_S]) ) {
85 mexErrMsgTxt("S should be a real scalar") ;
86 }
87
88 if( !uIsRealScalar(in[IN_SMIN]) ) {
89 mexErrMsgTxt("SMIN should be a real scalar") ;
90 }
91
92 if( !uIsRealScalar(in[IN_SIGMA0]) ) {
93 mexErrMsgTxt("SIGMA0 should be a real scalar") ;
94 }
95
96 if( !uIsRealMatrix(in[IN_P],3,-1)) {
97 mexErrMsgTxt("P should be a 3xK real matrix") ;
98 }
99
100 if(mxGetNumberOfDimensions(in[IN_G]) != 3) {
101 mexErrMsgTxt("SSO must be a three dimensional array") ;
102 }
103
104 dimensions = mxGetDimensions(in[IN_G]) ;
105 M = dimensions[0] ;
106 N = dimensions[1] ;
107 S = (int)(*mxGetPr(in[IN_S])) ;
108 smin = (int)(*mxGetPr(in[IN_SMIN])) ;
109 sigma0 = *mxGetPr(in[IN_SIGMA0]) ;
110
111 K = mxGetN(in[IN_P]) ;
112 P_pt = mxGetPr(in[IN_P]) ;
113 G_pt = mxGetPr(in[IN_G]) ;
114
115
116 /* If the input array is empty, then output an empty array as well. */
117 if(K == 0) {
118 out[OUT_Q] = mxCreateDoubleMatrix(4,0,mxREAL) ;
119 return ;
120 }
121
122 /* ------------------------------------------------------------------
123 * Do the job
124 * --------------------------------------------------------------- */
125 {
126 int p ;
127 const int yo = 1 ;
128 const int xo = M ;
129 const int so = M*N ;
130
131 int buffer_size = K*4 ;
132 double* buffer_start = (double*) mxMalloc( buffer_size *sizeof(double)) ;
133 double* buffer_iterator = buffer_start ;
134 double* buffer_end = buffer_start + buffer_size ;
135
136 for(p = 0 ; p < K ; ++p, TH_pt += 2) {
137 const double x = *P_pt++ ;
138 const double y = *P_pt++ ;
139 const double s = *P_pt++ ;
140 int xi = ((int) (x+0.5)) ; /* Round them off. */
141 int yi = ((int) (y+0.5)) ;
142 int si = ((int) (s+0.5)) - smin ;
143 int xs ;
144 int ys ;
145 double sigmaw = win_factor * sigma0 * pow(2, ((double)s) / S) ;
146 int W = (int) ceil(3.0 * sigmaw) ;
147 int bin ;
148 const double* pt ;
149
150 /* Make sure that the rounded off keypoint index is within bound.
151 */
152 if(xi < 0 ||
153 xi > N-1 ||
154 yi < 0 ||
155 yi > M-1 ||
156 si < 0 ||
157 si > dimensions[2]-1 ) {
158 mexPrintf("Dropping %d: W %d x %d y %d si [%d,%d,%d,%d]\n",p,W,xi,yi,si,M,N,dimensions[2]) ;
159 continue ;
160 }
161
162 /* Clear histogram buffer. */
163 {
164 int i ;
165 for(i = 0 ; i < NBINS ; ++i)
166 H_pt[i] = 0 ;
167 }
168
169 pt = G_pt + xi*xo + yi*yo + si*so ;
170
171#define at(dx,dy) (*(pt + (dx)*xo + (dy)*yo))
172
173 for(xs = max(-W, 1-xi) ; xs <= min(+W, N -2 -xi) ; ++xs) {
174 for(ys = max(-W, 1-yi) ; ys <= min(+W, M -2 -yi) ; ++ys) {
175 double Dx = 0.5 * ( at(xs+1,ys) - at(xs-1,ys) ) ;
176 double Dy = 0.5 * ( at(xs,ys+1) - at(xs,ys-1) ) ;
177 double dx = ((double)(xi+xs)) - x;
178 double dy = ((double)(yi+ys)) - y;
179
180 if(dx*dx + dy*dy > W*W+0.5) continue ;
181
182 {
183 double win = exp( - (dx*dx + dy*dy)/(2*sigmaw*sigmaw) ) ;
184 double mod = sqrt(Dx*Dx + Dy*Dy) ;
185 double theta = fmod(atan2(Dy, Dx) + 2*M_PI, 2*M_PI) ;
186 bin = (int) floor( NBINS * theta / (2*M_PI) ) ;
187 H_pt[bin] += mod*win ;
188 }
189 }
190 }
191
192 /* Smooth histogram */
193 {
194 int iter, i ;
195 for (iter = 0; iter < 6; iter++) {
196 double prev;
197 prev = H_pt[NBINS-1];
198 for (i = 0; i < NBINS; i++) {
199 float newh = (prev + H_pt[i] + H_pt[(i+1) % NBINS]) / 3.0;
200 prev = H_pt[i] ;
201 H_pt[i] = newh ;
202 }
203 }
204 }
205
206 /* Find most strong peaks. */
207 {
208 int i ;
209 double maxh = H_pt[0] ;
210 for(i = 1 ; i < NBINS ; ++i)
211 maxh = max(maxh, H_pt[i]) ;
212
213 for(i = 0 ; i < NBINS ; ++i) {
214 double h0 = H_pt[i] ;
215 double hm = H_pt[(i-1+NBINS) % NBINS] ;
216 double hp = H_pt[(i+1+NBINS) % NBINS] ;
217
218 if( h0 > 0.8*maxh && h0 > hm && h0 > hp ) {
219
220 double di = -0.5 * (hp-hm) / (hp+hm-2*h0) ; /*di=0;*/
221 double th = 2*M_PI*(i+di+0.5)/NBINS ;
222
223 if( buffer_iterator == buffer_end ) {
224 int offset = buffer_iterator - buffer_start ;
225 buffer_size += 4*max(1, K/16) ;
226 buffer_start = (double*) mxRealloc(buffer_start,
227 buffer_size*sizeof(double)) ;
228 buffer_end = buffer_start + buffer_size ;
229 buffer_iterator = buffer_start + offset ;
230 }
231
232 *buffer_iterator++ = x ;
233 *buffer_iterator++ = y ;
234 *buffer_iterator++ = s ;
235 *buffer_iterator++ = th ;
236 }
237 } /* Scan histogram */
238 } /* Find peaks */
239 }
240
241 /* Save back the result. */
242 {
243 double* result ;
244 int NL = (buffer_iterator - buffer_start)/4 ;
245 out[OUT_Q] = mxCreateDoubleMatrix(4, NL, mxREAL) ;
246 result = mxGetPr(out[OUT_Q]);
247 memcpy(result, buffer_start, sizeof(double) * 4 * NL) ;
248 }
249 mxFree(buffer_start) ;
250 }
251}
generated by cgit v1.2.2 (git 2.25.0) at 2025-04-30 17:37:46 -0400