Bug Summary

File:ccv_sift.c
Warning:line 258, column 15
Dereference of null pointer

Annotated Source Code

Press '?' to see keyboard shortcuts

clang -cc1 -cc1 -triple x86_64-unknown-linux-gnu -analyze -disable-free -clear-ast-before-backend -disable-llvm-verifier -discard-value-names -main-file-name ccv_sift.c -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -pic-is-pie -mframe-pointer=none -menable-no-infs -menable-no-nans -fapprox-func -funsafe-math-optimizations -fno-signed-zeros -mreassociate -freciprocal-math -ffp-contract=fast -fno-rounding-math -ffast-math -ffinite-math-only -complex-range=basic -mconstructor-aliases -funwind-tables=2 -target-cpu x86-64 -target-feature +sse2 -tune-cpu generic -debugger-tuning=gdb -fdebug-compilation-dir=/home/liu/actions-runner/_work/ccv/ccv/lib -fcoverage-compilation-dir=/home/liu/actions-runner/_work/ccv/ccv/lib -resource-dir /usr/local/lib/clang/19 -I . -I /usr/local/cuda/include -D HAVE_CBLAS -D HAVE_LIBPNG -D HAVE_LIBJPEG -D HAVE_FFTW3 -D HAVE_PTHREAD -D HAVE_LIBLINEAR -D HAVE_TESSERACT -D HAVE_AVCODEC -D HAVE_AVFORMAT -D HAVE_AVUTIL -D HAVE_SWSCALE -D HAVE_SSE2 -D HAVE_GSL -D HAVE_CUDA -D HAVE_CUDNN -D HAVE_NCCL -D USE_SYSTEM_CUB -I /usr/local/include -internal-isystem /usr/local/lib/clang/19/include -internal-isystem /usr/local/include -internal-isystem /usr/lib/gcc/x86_64-linux-gnu/12/../../../../x86_64-linux-gnu/include -internal-externc-isystem /usr/include/x86_64-linux-gnu -internal-externc-isystem /include -internal-externc-isystem /usr/include -O3 -ferror-limit 19 -fgnuc-version=4.2.1 -fskip-odr-check-in-gmf -vectorize-loops -vectorize-slp -analyzer-output=html -faddrsig -D__GCC_HAVE_DWARF2_CFI_ASM=1 -o /home/liu/actions-runner/_work/ccv/ccv/_analyze/2024-10-18-174139-54102-1 -x c ccv_sift.c
1/* The code is adopted from VLFeat with heavily rewrite.
2 * The original code is licenced under 2-clause BSD license,
3 * should be compatible with New BSD Licence used by ccv.
4 * The original Copyright:
5 *
6 * Copyright (C) 2007-12, Andrea Vedaldi and Brian Fulkerson
7 * All rights reserved.
8 *
9 * Redistribution and use in source and binary forms, with or without
10 * modification, are permitted provided that the following conditions are
11 * met:
12 * 1. Redistributions of source code must retain the above copyright
13 * notice, this list of conditions and the following disclaimer.
14 * 2. Redistributions in binary form must reproduce the above copyright
15 * notice, this list of conditions and the following disclaimer in the
16 * documentation and/or other materials provided with the
17 * distribution.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
22 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
23 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
24 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
25 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
26 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
27 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
28 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
29 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30 */
31
32#include "ccv.h"
33#include "ccv_internal.h"
34
35const ccv_sift_param_t ccv_sift_default_params = {
36 .noctaves = 3,
37 .nlevels = 6,
38 .up2x = 1,
39 .edge_threshold = 10,
40 .norm_threshold = 0,
41 .peak_threshold = 0,
42};
43
44inline static double _ccv_keypoint_interpolate(float N9[3][9], int ix, int iy, int is, ccv_keypoint_t* kp)
45{
46 double Dxx = N9[1][3] - 2 * N9[1][4] + N9[1][5];
47 double Dyy = N9[1][1] - 2 * N9[1][4] + N9[1][7];
48 double Dxy = (N9[1][8] - N9[1][6] - N9[1][2] + N9[1][0]) * 0.25;
49 double score = (Dxx + Dyy) * (Dxx + Dyy) / (Dxx * Dyy - Dxy * Dxy);
50 double Dx = (N9[1][5] - N9[1][3]) * 0.5;
51 double Dy = (N9[1][7] - N9[1][1]) * 0.5;
52 double Ds = (N9[2][4] - N9[0][4]) * 0.5;
53 double Dxs = (N9[2][5] + N9[0][3] - N9[2][3] - N9[0][5]) * 0.25;
54 double Dys = (N9[2][7] + N9[0][1] - N9[2][1] - N9[0][7]) * 0.25;
55 double Dss = N9[0][4] - 2 * N9[1][4] + N9[2][4];
56 double A[3][3] = { { Dxx, Dxy, Dxs },
57 { Dxy, Dyy, Dys },
58 { Dxs, Dys, Dss } };
59 double b[3] = { -Dx, -Dy, -Ds };
60 /* Gauss elimination */
61 int i, j, ii, jj;
62 for(j = 0; j < 3; j++)
63 {
64 double maxa = 0;
65 double maxabsa = 0;
66 int maxi = j;
67 double tmp;
68
69 /* look for the maximally stable pivot */
70 for (i = j; i < 3; i++)
71 {
72 double a = A[i][j];
73 double absa = fabs(a);
74 if (absa > maxabsa)
75 {
76 maxa = a;
77 maxabsa = absa;
78 maxi = i;
79 }
80 }
81
82 /* if singular give up */
83 if (maxabsa < 1e-10f)
84 {
85 b[0] = b[1] = b[2] = 0;
86 break;
87 }
88
89 i = maxi;
90
91 /* swap j-th row with i-th row and normalize j-th row */
92 for(jj = j; jj < 3; jj++)
93 {
94 tmp = A[i][jj];
95 A[i][jj] = A[j][jj];
96 A[j][jj] = tmp;
97 A[j][jj] /= maxa;
98 }
99 tmp = b[j];
100 b[j] = b[i];
101 b[i] = tmp;
102 b[j] /= maxa;
103
104 /* elimination */
105 for (ii = j + 1; ii < 3; ii++)
106 {
107 double x = A[ii][j];
108 for (jj = j; jj < 3; jj++)
109 A[ii][jj] -= x * A[j][jj];
110 b[ii] -= x * b[j];
111 }
112 }
113
114 /* backward substitution */
115 for (i = 2; i > 0; i--)
116 {
117 double x = b[i];
118 for (ii = i - 1; ii >= 0; ii--)
119 b[ii] -= x * A[ii][i];
120 }
121 kp->x = ix + ccv_min(ccv_max(b[0], -1), 1)({ typeof (({ typeof (b[0]) _a = (b[0]); typeof (-1) _b = (-1
); (_a > _b) ? _a : _b; })) _a = (({ typeof (b[0]) _a = (b
[0]); typeof (-1) _b = (-1); (_a > _b) ? _a : _b; })); typeof
(1) _b = (1); (_a < _b) ? _a : _b; });
122 kp->y = iy + ccv_min(ccv_max(b[1], -1), 1)({ typeof (({ typeof (b[1]) _a = (b[1]); typeof (-1) _b = (-1
); (_a > _b) ? _a : _b; })) _a = (({ typeof (b[1]) _a = (b
[1]); typeof (-1) _b = (-1); (_a > _b) ? _a : _b; })); typeof
(1) _b = (1); (_a < _b) ? _a : _b; });
123 kp->regular.scale = is + b[2];
124 return score;
125}
126
127static float _ccv_mod_2pi(float x)
128{
129 while (x > 2 * CCV_PI(3.141592653589793))
130 x -= 2 * CCV_PI(3.141592653589793);
131 while (x < 0)
132 x += 2 * CCV_PI(3.141592653589793);
133 return x;
134}
135
136static int _ccv_floor(float x)
137{
138 int xi = (int)x;
139 if (x >= 0 || (float)xi == x)
140 return xi;
141 return xi - 1;
142}
143
144#define EXPN_SZ256 256 /* fast_expn table size */
145#define EXPN_MAX25.0 25.0 /* fast_expn table max */
146static double _ccv_expn_tab[EXPN_SZ256 + 1]; /* fast_expn table */
147static int _ccv_expn_init = 0;
148
149static inline double _ccv_expn(double x)
150{
151 double a, b, r;
152 int i;
153 assert(0 <= x && x <= EXPN_MAX)((void) sizeof ((0 <= x && x <= 25.0) ? 1 : 0),
__extension__ ({ if (0 <= x && x <= 25.0) ; else
__assert_fail ("0 <= x && x <= EXPN_MAX", "ccv_sift.c"
, 153, __extension__ __PRETTY_FUNCTION__); }));
154 if (x > EXPN_MAX25.0)
155 return 0.0;
156 x *= EXPN_SZ256 / EXPN_MAX25.0;
157 i = (int)x;
158 r = x - i;
159 a = _ccv_expn_tab[i];
160 b = _ccv_expn_tab[i + 1];
161 return a + r * (b - a);
162}
163
164static void _ccv_precomputed_expn()
165{
166 int i;
167 for(i = 0; i < EXPN_SZ256 + 1; i++)
168 _ccv_expn_tab[i] = exp(-(double)i * (EXPN_MAX25.0 / EXPN_SZ256));
169 _ccv_expn_init = 1;
170}
171
172void ccv_sift(ccv_dense_matrix_t* a, ccv_array_t** _keypoints, ccv_dense_matrix_t** _desc, int type, ccv_sift_param_t params)
173{
174 assert(CCV_GET_CHANNEL(a->type) == CCV_C1)((void) sizeof ((((a->type) & 0xFFF) == CCV_C1) ? 1 : 0
), __extension__ ({ if (((a->type) & 0xFFF) == CCV_C1)
; else __assert_fail ("CCV_GET_CHANNEL(a->type) == CCV_C1"
, "ccv_sift.c", 174, __extension__ __PRETTY_FUNCTION__); }));
1
Assuming the condition is true
→
2
←
Taking true branch
→
175 ccv_dense_matrix_t** g = (ccv_dense_matrix_t**)alloca(sizeof(ccv_dense_matrix_t*) * (params.nlevels + 1) * (params.up2x ? params.noctaves + 1 : params.noctaves))__builtin_alloca (sizeof(ccv_dense_matrix_t*) * (params.nlevels
+ 1) * (params.up2x ? params.noctaves + 1 : params.noctaves)
);
3
←
Assuming field 'up2x' is 0
→
4
←
'?' condition is false
→
176 memset(g, 0, sizeof(ccv_dense_matrix_t*) * (params.nlevels + 1) * (params.up2x
4.1
Field 'up2x' is 0
 ? params.noctaves + 1 : params.noctaves));
5
←
'?' condition is false
→
177 ccv_dense_matrix_t** dog = (ccv_dense_matrix_t**)alloca(sizeof(ccv_dense_matrix_t*) * (params.nlevels - 1) * (params.up2x ? params.noctaves + 1 : params.noctaves))__builtin_alloca (sizeof(ccv_dense_matrix_t*) * (params.nlevels
- 1) * (params.up2x ? params.noctaves + 1 : params.noctaves)
);
6
←
'?' condition is false
→
178 memset(dog, 0, sizeof(ccv_dense_matrix_t*) * (params.nlevels - 1) * (params.up2x
6.1
Field 'up2x' is 0
 ? params.noctaves + 1 : params.noctaves));
7
←
'?' condition is false
→
8
←
Storing null pointer value
→
179 ccv_dense_matrix_t** th = (ccv_dense_matrix_t**)alloca(sizeof(ccv_dense_matrix_t*) * (params.nlevels - 3) * (params.up2x ? params.noctaves + 1 : params.noctaves))__builtin_alloca (sizeof(ccv_dense_matrix_t*) * (params.nlevels
- 3) * (params.up2x ? params.noctaves + 1 : params.noctaves)
);
9
←
'?' condition is false
→
180 memset(th, 0, sizeof(ccv_dense_matrix_t*) * (params.nlevels - 3) * (params.up2x
9.1
Field 'up2x' is 0
 ? params.noctaves + 1 : params.noctaves));
10
←
'?' condition is false
→
181 ccv_dense_matrix_t** md = (ccv_dense_matrix_t**)alloca(sizeof(ccv_dense_matrix_t*) * (params.nlevels - 3) * (params.up2x ? params.noctaves + 1 : params.noctaves))__builtin_alloca (sizeof(ccv_dense_matrix_t*) * (params.nlevels
- 3) * (params.up2x ? params.noctaves + 1 : params.noctaves)
);
11
←
'?' condition is false
→
182 memset(md, 0, sizeof(ccv_dense_matrix_t*) * (params.nlevels - 3) * (params.up2x
11.1
Field 'up2x' is 0
 ? params.noctaves + 1 : params.noctaves));
12
←
'?' condition is false
→
183 if (params.up2x
12.1
Field 'up2x' is 0
)
13
←
Taking false branch
→
184 {
185 g += params.nlevels + 1;
186 dog += params.nlevels - 1;
187 th += params.nlevels - 3;
188 md += params.nlevels - 3;
189 }
190 ccv_array_t* keypoints = *_keypoints;
191 int custom_keypoints = 0;
192 if (keypoints == 0)
14
←
Assuming 'keypoints' is equal to null
→
15
←
Taking true branch
→
193 keypoints = *_keypoints = ccv_array_new(sizeof(ccv_keypoint_t), 10, 0);
194 else
195 custom_keypoints = 1;
196 int i, j, k, x, y;
197 double sigma0 = 1.6;
198 double sigmak = pow(2.0, 1.0 / (params.nlevels - 3));
199 double dsigma0 = sigma0 * sigmak * sqrt(1.0 - 1.0 / (sigmak * sigmak));
200 if (params.up2x
15.1
Field 'up2x' is 0
)
16
←
Taking false branch
→
201 {
202 ccv_sample_up(a, &g[-(params.nlevels + 1)], 0, 0, 0);
203 /* since there is a gaussian filter in sample_up function already,
204 * the default sigma for upsampled image is sqrt(2) */
205 double sd = sqrt(sigma0 * sigma0 - 2.0);
206 ccv_blur(g[-(params.nlevels + 1)], &g[-(params.nlevels + 1) + 1], CCV_32F | CCV_C1, sd);
207 ccv_matrix_free(g[-(params.nlevels + 1)]);
208 for (j = 1; j < params.nlevels; j++)
209 {
210 sd = dsigma0 * pow(sigmak, j - 1);
211 ccv_blur(g[-(params.nlevels + 1) + j], &g[-(params.nlevels + 1) + j + 1], 0, sd);
212 ccv_subtract(g[-(params.nlevels + 1) + j + 1], g[-(params.nlevels + 1) + j], (ccv_matrix_t**)&dog[-(params.nlevels - 1) + j - 1], 0);
213 if (j > 1 && j < params.nlevels - 1)
214 ccv_gradient(g[-(params.nlevels + 1) + j], &th[-(params.nlevels - 3) + j - 2], 0, &md[-(params.nlevels - 3) + j - 2], 0, 1, 1);
215 ccv_matrix_free(g[-(params.nlevels + 1) + j]);
216 }
217 ccv_matrix_free(g[-1]);
218 }
219 double sd = sqrt(sigma0 * sigma0 - 0.25);
220 g[0] = a;
221 /* generate gaussian pyramid (g, dog) & gradient pyramid (th, md) */
222 ccv_blur(g[0], &g[1], CCV_32F | CCV_C1, sd);
223 for (j = 1; j < params.nlevels; j++)
17
←
Assuming 'j' is >= field 'nlevels'
→
18
←
Loop condition is false. Execution continues on line 232
→
224 {
225 sd = dsigma0 * pow(sigmak, j - 1);
226 ccv_blur(g[j], &g[j + 1], 0, sd);
227 ccv_subtract(g[j + 1], g[j], (ccv_matrix_t**)&dog[j - 1], 0);
228 if (j > 1 && j < params.nlevels - 1)
229 ccv_gradient(g[j], &th[j - 2], 0, &md[j - 2], 0, 1, 1);
230 ccv_matrix_free(g[j]);
231 }
232 ccv_matrix_free(g[params.nlevels]);
233 for (i = 1; i < params.noctaves; i++)
19
←
Assuming 'i' is >= field 'noctaves'
→
20
←
Loop condition is false. Execution continues on line 251
→
234 {
235 ccv_sample_down(g[(i - 1) * (params.nlevels + 1)], &g[i * (params.nlevels + 1)], 0, 0, 0);
236 if (i - 1 > 0)
237 ccv_matrix_free(g[(i - 1) * (params.nlevels + 1)]);
238 sd = sqrt(sigma0 * sigma0 - 0.25);
239 ccv_blur(g[i * (params.nlevels + 1)], &g[i * (params.nlevels + 1) + 1], CCV_32F | CCV_C1, sd);
240 for (j = 1; j < params.nlevels; j++)
241 {
242 sd = dsigma0 * pow(sigmak, j - 1);
243 ccv_blur(g[i * (params.nlevels + 1) + j], &g[i * (params.nlevels + 1) + j + 1], 0, sd);
244 ccv_subtract(g[i * (params.nlevels + 1) + j + 1], g[i * (params.nlevels + 1) + j], (ccv_matrix_t**)&dog[i * (params.nlevels - 1) + j - 1], 0);
245 if (j > 1 && j < params.nlevels - 1)
246 ccv_gradient(g[i * (params.nlevels + 1) + j], &th[i * (params.nlevels - 3) + j - 2], 0, &md[i * (params.nlevels - 3) + j - 2], 0, 1, 1);
247 ccv_matrix_free(g[i * (params.nlevels + 1) + j]);
248 }
249 ccv_matrix_free(g[i * (params.nlevels + 1) + params.nlevels]);
250 }
251 ccv_matrix_free(g[(params.noctaves - 1) * (params.nlevels + 1)]);
252 if (!custom_keypoints
20.1
'custom_keypoints' is 0
)
253 {
254 /* detect keypoint */
255 for (i = (params.up2x
21.1
Field 'up2x' is 0
 ? -1 : 0); i < params.noctaves; i++)
21
←
Taking true branch
→
22
←
'?' condition is false
→
23
←
The value 0 is assigned to 'i'
→
24
←
Assuming 'i' is < field 'noctaves'
→
25
←
Loop condition is true. Entering loop body
→
256 {
257 double s = pow(2.0, i);
258 int rows = dog[i * (params.nlevels - 1)]->rows;
26
←
Dereference of null pointer
259 int cols = dog[i * (params.nlevels - 1)]->cols;
260 for (j = 1; j < params.nlevels - 2; j++)
261 {
262 float* bf = dog[i * (params.nlevels - 1) + j - 1]->data.f32 + cols;
263 float* cf = dog[i * (params.nlevels - 1) + j]->data.f32 + cols;
264 float* uf = dog[i * (params.nlevels - 1) + j + 1]->data.f32 + cols;
265 for (y = 1; y < rows - 1; y++)
266 {
267 for (x = 1; x < cols - 1; x++)
268 {
269 float v = cf[x];
270#define locality_if(CMP, SGN) \
271 (v CMP ## = SGN params.peak_threshold && v CMP cf[x - 1] && v CMP cf[x + 1] && \
272 v CMP cf[x - cols - 1] && v CMP cf[x - cols] && v CMP cf[x - cols + 1] && \
273 v CMP cf[x + cols - 1] && v CMP cf[x + cols] && v CMP cf[x + cols + 1] && \
274 v CMP bf[x - 1] && v CMP bf[x] && v CMP bf[x + 1] && \
275 v CMP bf[x - cols - 1] && v CMP bf[x - cols] && v CMP bf[x - cols + 1] && \
276 v CMP bf[x + cols - 1] && v CMP bf[x + cols] && v CMP bf[x + cols + 1] && \
277 v CMP uf[x - 1] && v CMP uf[x] && v CMP uf[x + 1] && \
278 v CMP uf[x - cols - 1] && v CMP uf[x - cols] && v CMP uf[x - cols + 1] && \
279 v CMP uf[x + cols - 1] && v CMP uf[x + cols] && v CMP uf[x + cols + 1])
280 if (locality_if(<, -) || locality_if(>, +))
281 {
282 ccv_keypoint_t kp;
283 int ix = x, iy = y;
284 double score = -1;
285 int cvg = 0;
286 int offset = ix + (iy - y) * cols;
287 /* iteratively converge to meet subpixel accuracy */
288 for (k = 0; k < 5; k++)
289 {
290 offset = ix + (iy - y) * cols;
291 float N9[3][9] = { { bf[offset - cols - 1], bf[offset - cols], bf[offset - cols + 1],
292 bf[offset - 1], bf[offset], bf[offset + 1],
293 bf[offset + cols - 1], bf[offset + cols], bf[offset + cols + 1] },
294 { cf[offset - cols - 1], cf[offset - cols], cf[offset - cols + 1],
295 cf[offset - 1], cf[offset], cf[offset + 1],
296 cf[offset + cols - 1], cf[offset + cols], cf[offset + cols + 1] },
297 { uf[offset - cols - 1], uf[offset - cols], uf[offset - cols + 1],
298 uf[offset - 1], uf[offset], uf[offset + 1],
299 uf[offset + cols - 1], uf[offset + cols], uf[offset + cols + 1] } };
300 score = _ccv_keypoint_interpolate(N9, ix, iy, j, &kp);
301 if (kp.x >= 1 && kp.x <= cols - 2 && kp.y >= 1 && kp.y <= rows - 2)
302 {
303 int nx = (int)(kp.x + 0.5);
304 int ny = (int)(kp.y + 0.5);
305 if (ix == nx && iy == ny)
306 break;
307 ix = nx;
308 iy = ny;
309 } else {
310 cvg = -1;
311 break;
312 }
313 }
314 if (cvg == 0 && fabs(cf[offset]) > params.peak_threshold && score >= 0 && score < (params.edge_threshold + 1) * (params.edge_threshold + 1) / params.edge_threshold && kp.regular.scale > 0 && kp.regular.scale < params.nlevels - 1)
315 {
316 kp.x *= s;
317 kp.y *= s;
318 kp.octave = i;
319 kp.level = j;
320 kp.regular.scale = sigma0 * sigmak * pow(2.0, kp.regular.scale / (double)(params.nlevels - 3));
321 ccv_array_push(keypoints, &kp);
322 }
323 }
324#undef locality_if
325 }
326 bf += cols;
327 cf += cols;
328 uf += cols;
329 }
330 }
331 }
332 }
333 /* repeatable orientation/angle (p.s. it will push more keypoints (with different angles) to array) */
334 float const winf = 1.5;
335 double bins[36];
336 int kpnum = keypoints->rnum;
337 if (!_ccv_expn_init)
338 _ccv_precomputed_expn();
339 for (i = 0; i < kpnum; i++)
340 {
341 ccv_keypoint_t* kp = (ccv_keypoint_t*)ccv_array_get(keypoints, i)((void*)(((char*)((keypoints)->data)) + (size_t)(keypoints
)->rsize * (size_t)(i)));
342 float ds = pow(2.0, kp->octave);
343 float dx = kp->x / ds;
344 float dy = kp->y / ds;
345 int ix = (int)(dx + 0.5);
346 int iy = (int)(dy + 0.5);
347 float const sigmaw = winf * kp->regular.scale;
348 int wz = ccv_max((int)(3.0 * sigmaw + 0.5), 1)({ typeof ((int)(3.0 * sigmaw + 0.5)) _a = ((int)(3.0 * sigmaw
+ 0.5)); typeof (1) _b = (1); (_a > _b) ? _a : _b; });
349 ccv_dense_matrix_t* tho = th[kp->octave * (params.nlevels - 3) + kp->level - 1];
350 ccv_dense_matrix_t* mdo = md[kp->octave * (params.nlevels - 3) + kp->level - 1];
351 assert(tho->rows == mdo->rows && tho->cols == mdo->cols)((void) sizeof ((tho->rows == mdo->rows && tho->
cols == mdo->cols) ? 1 : 0), __extension__ ({ if (tho->
rows == mdo->rows && tho->cols == mdo->cols)
; else __assert_fail ("tho->rows == mdo->rows && tho->cols == mdo->cols"
, "ccv_sift.c", 351, __extension__ __PRETTY_FUNCTION__); }));
352 if (ix >= 0 && ix < tho->cols && iy >=0 && iy < tho->rows)
353 {
354 float* theta = tho->data.f32 + ccv_max(iy - wz, 0)({ typeof (iy - wz) _a = (iy - wz); typeof (0) _b = (0); (_a >
_b) ? _a : _b; }) * tho->cols;
355 float* magnitude = mdo->data.f32 + ccv_max(iy - wz, 0)({ typeof (iy - wz) _a = (iy - wz); typeof (0) _b = (0); (_a >
_b) ? _a : _b; }) * mdo->cols;
356 memset(bins, 0, 36 * sizeof(double));
357 /* oriented histogram with bilinear interpolation */
358 for (y = ccv_max(iy - wz, 0)({ typeof (iy - wz) _a = (iy - wz); typeof (0) _b = (0); (_a >
_b) ? _a : _b; }); y <= ccv_min(iy + wz, tho->rows - 1)({ typeof (iy + wz) _a = (iy + wz); typeof (tho->rows - 1)
_b = (tho->rows - 1); (_a < _b) ? _a : _b; }); y++)
359 {
360 for (x = ccv_max(ix - wz, 0)({ typeof (ix - wz) _a = (ix - wz); typeof (0) _b = (0); (_a >
_b) ? _a : _b; }); x <= ccv_min(ix + wz, tho->cols - 1)({ typeof (ix + wz) _a = (ix + wz); typeof (tho->cols - 1)
_b = (tho->cols - 1); (_a < _b) ? _a : _b; }); x++)
361 {
362 float r2 = (x - dx) * (x - dx) + (y - dy) * (y - dy);
363 if (r2 > wz * wz + 0.6)
364 continue;
365 float weight = _ccv_expn(r2 / (2.0 * sigmaw * sigmaw));
366 float fbin = theta[x] * 0.1;
367 int ibin = _ccv_floor(fbin - 0.5);
368 float rbin = fbin - ibin - 0.5;
369 /* bilinear interpolation */
370 bins[(ibin + 36) % 36] += (1 - rbin) * magnitude[x] * weight;
371 bins[(ibin + 1) % 36] += rbin * magnitude[x] * weight;
372 }
373 theta += tho->cols;
374 magnitude += mdo->cols;
375 }
376 /* smoothing histogram */
377 for (j = 0; j < 6; j++)
378 {
379 double first = bins[0];
380 double prev = bins[35];
381 for (k = 0; k < 35; k++)
382 {
383 double nb = (prev + bins[k] + bins[k + 1]) / 3.0;
384 prev = bins[k];
385 bins[k] = nb;
386 }
387 bins[35] = (prev + bins[35] + first) / 3.0;
388 }
389 int maxib = 0;
390 for (j = 1; j < 36; j++)
391 if (bins[j] > bins[maxib])
392 maxib = j;
393 double maxb = bins[maxib];
394 double bm = bins[(maxib + 35) % 36];
395 double bp = bins[(maxib + 1) % 36];
396 double di = -0.5 * (bp - bm) / (bp + bm - 2 * maxb);
397 kp->regular.angle = 2 * CCV_PI(3.141592653589793) * (maxib + di + 0.5) / 36.0;
398 maxb *= 0.8;
399 for (j = 0; j < 36; j++)
400 if (j != maxib)
401 {
402 bm = bins[(j + 35) % 36];
403 bp = bins[(j + 1) % 36];
404 if (bins[j] > maxb && bins[j] > bm && bins[j] > bp)
405 {
406 di = -0.5 * (bp - bm) / (bp + bm - 2 * bins[j]);
407 ccv_keypoint_t nkp = *kp;
408 nkp.regular.angle = 2 * CCV_PI(3.141592653589793) * (j + di + 0.5) / 36.0;
409 ccv_array_push(keypoints, &nkp);
410 }
411 }
412 }
413 }
414 /* calculate descriptor */
415 if (_desc != 0)
416 {
417 ccv_dense_matrix_t* desc = *_desc = ccv_dense_matrix_new(keypoints->rnum, 128, CCV_32F | CCV_C1, 0, 0);
418 float* fdesc = desc->data.f32;
419 memset(fdesc, 0, sizeof(float) * keypoints->rnum * 128);
420 for (i = 0; i < keypoints->rnum; i++)
421 {
422 ccv_keypoint_t* kp = (ccv_keypoint_t*)ccv_array_get(keypoints, i)((void*)(((char*)((keypoints)->data)) + (size_t)(keypoints
)->rsize * (size_t)(i)));
423 float ds = pow(2.0, kp->octave);
424 float dx = kp->x / ds;
425 float dy = kp->y / ds;
426 int ix = (int)(dx + 0.5);
427 int iy = (int)(dy + 0.5);
428 double SBP = 3.0 * kp->regular.scale;
429 int wz = ccv_max((int)(SBP * sqrt(2.0) * 2.5 + 0.5), 1)({ typeof ((int)(SBP * sqrt(2.0) * 2.5 + 0.5)) _a = ((int)(SBP
* sqrt(2.0) * 2.5 + 0.5)); typeof (1) _b = (1); (_a > _b)
? _a : _b; });
430 ccv_dense_matrix_t* tho = th[kp->octave * (params.nlevels - 3) + kp->level - 1];
431 ccv_dense_matrix_t* mdo = md[kp->octave * (params.nlevels - 3) + kp->level - 1];
432 assert(tho->rows == mdo->rows && tho->cols == mdo->cols)((void) sizeof ((tho->rows == mdo->rows && tho->
cols == mdo->cols) ? 1 : 0), __extension__ ({ if (tho->
rows == mdo->rows && tho->cols == mdo->cols)
; else __assert_fail ("tho->rows == mdo->rows && tho->cols == mdo->cols"
, "ccv_sift.c", 432, __extension__ __PRETTY_FUNCTION__); }));
433 assert(ix >= 0 && ix < tho->cols && iy >=0 && iy < tho->rows)((void) sizeof ((ix >= 0 && ix < tho->cols &&
iy >=0 && iy < tho->rows) ? 1 : 0), __extension__
({ if (ix >= 0 && ix < tho->cols &&
iy >=0 && iy < tho->rows) ; else __assert_fail
("ix >= 0 && ix < tho->cols && iy >=0 && iy < tho->rows"
, "ccv_sift.c", 433, __extension__ __PRETTY_FUNCTION__); }));
434 float* theta = tho->data.f32 + ccv_max(iy - wz, 0)({ typeof (iy - wz) _a = (iy - wz); typeof (0) _b = (0); (_a >
_b) ? _a : _b; }) * tho->cols;
435 float* magnitude = mdo->data.f32 + ccv_max(iy - wz, 0)({ typeof (iy - wz) _a = (iy - wz); typeof (0) _b = (0); (_a >
_b) ? _a : _b; }) * mdo->cols;
436 float ca = cos(kp->regular.angle);
437 float sa = sin(kp->regular.angle);
438 float sigmaw = 2.0;
439 /* sidenote: NBP = 4, NBO = 8 */
440 for (y = ccv_max(iy - wz, 0)({ typeof (iy - wz) _a = (iy - wz); typeof (0) _b = (0); (_a >
_b) ? _a : _b; }); y <= ccv_min(iy + wz, tho->rows - 1)({ typeof (iy + wz) _a = (iy + wz); typeof (tho->rows - 1)
_b = (tho->rows - 1); (_a < _b) ? _a : _b; }); y++)
441 {
442 for (x = ccv_max(ix - wz, 0)({ typeof (ix - wz) _a = (ix - wz); typeof (0) _b = (0); (_a >
_b) ? _a : _b; }); x <= ccv_min(ix + wz, tho->cols - 1)({ typeof (ix + wz) _a = (ix + wz); typeof (tho->cols - 1)
_b = (tho->cols - 1); (_a < _b) ? _a : _b; }); x++)
443 {
444 float nx = (ca * (x - dx) + sa * (y - dy)) / SBP;
445 float ny = (-sa * (x - dx) + ca * (y - dy)) / SBP;
446 float nt = 8.0 * _ccv_mod_2pi(theta[x] * CCV_PI(3.141592653589793) / 180.0 - kp->regular.angle) / (2.0 * CCV_PI(3.141592653589793));
447 float weight = _ccv_expn((nx * nx + ny * ny) / (2.0 * sigmaw * sigmaw));
448 int binx = _ccv_floor(nx - 0.5);
449 int biny = _ccv_floor(ny - 0.5);
450 int bint = _ccv_floor(nt);
451 float rbinx = nx - (binx + 0.5);
452 float rbiny = ny - (biny + 0.5);
453 float rbint = nt - bint;
454 int dbinx, dbiny, dbint;
455 /* Distribute the current sample into the 8 adjacent bins*/
456 for(dbinx = 0; dbinx < 2; dbinx++)
457 for(dbiny = 0; dbiny < 2; dbiny++)
458 for(dbint = 0; dbint < 2; dbint++)
459 if (binx + dbinx >= -2 && binx + dbinx < 2 && biny + dbiny >= -2 && biny + dbiny < 2)
460 fdesc[(2 + biny + dbiny) * 32 + (2 + binx + dbinx) * 8 + (bint + dbint) % 8] += weight * magnitude[x] * fabs(1 - dbinx - rbinx) * fabs(1 - dbiny - rbiny) * fabs(1 - dbint - rbint);
461 }
462 theta += tho->cols;
463 magnitude += mdo->cols;
464 }
465 ccv_dense_matrix_t tm = ccv_dense_matrix(1, 128, CCV_32F | CCV_C1, fdesc, 0);
466 ccv_dense_matrix_t* tmp = &tm;
467 double norm = ccv_normalize(&tm, (ccv_matrix_t**)&tmp, 0, CCV_L2_NORM);
468 int num = (ccv_min(iy + wz, tho->rows - 1)({ typeof (iy + wz) _a = (iy + wz); typeof (tho->rows - 1)
_b = (tho->rows - 1); (_a < _b) ? _a : _b; }) - ccv_max(iy - wz, 0)({ typeof (iy - wz) _a = (iy - wz); typeof (0) _b = (0); (_a >
_b) ? _a : _b; }) + 1) * (ccv_min(ix + wz, tho->cols - 1)({ typeof (ix + wz) _a = (ix + wz); typeof (tho->cols - 1)
_b = (tho->cols - 1); (_a < _b) ? _a : _b; }) - ccv_max(ix - wz, 0)({ typeof (ix - wz) _a = (ix - wz); typeof (0) _b = (0); (_a >
_b) ? _a : _b; }) + 1);
469 if (params.norm_threshold && norm < params.norm_threshold * num)
470 {
471 for (j = 0; j < 128; j++)
472 fdesc[j] = 0;
473 } else {
474 for (j = 0; j < 128; j++)
475 if (fdesc[j] > 0.2)
476 fdesc[j] = 0.2;
477 ccv_normalize(&tm, (ccv_matrix_t**)&tmp, 0, CCV_L2_NORM);
478 }
479 fdesc += 128;
480 }
481 }
482 for (i = (params.up2x ? -(params.nlevels - 1) : 0); i < (params.nlevels - 1) * params.noctaves; i++)
483 ccv_matrix_free(dog[i]);
484 for (i = (params.up2x ? -(params.nlevels - 3) : 0); i < (params.nlevels - 3) * params.noctaves; i++)
485 {
486 ccv_matrix_free(th[i]);
487 ccv_matrix_free(md[i]);
488 }
489}