1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
|
#include "sift.h"
/**
SIFTLOCALMAX Find local maximizers
Returns the indexes of the local maximizers of
the 3-dimensional array F.
Say, we have a Q-dimensional array F.
A local maximizer is an element whose value is greater than the
value of all its neighbors. The neighbors of an element i1...iQ
are the subscripts j1...jQ such that iq-1 <= jq <= iq (excluding
i1...iQ itself). For example, if Q=1 the neighbors of an element
are its predecessor and successor in the linear order; if Q=2, its
neighbors are the elements immediately to its north, south, west,
est, north-west, north-est, south-west and south-est
(8-neighborhood).
Points on the boundary of F are ignored (and never selected as
local maximizers).
SEL=SIFTLOCALMAX(F,THRESH) accepts an element as a mazimizer only
if it is at least THRES greater than all its neighbors.
SEL=SIFTLOCALMAX(F,THRESH,P) look for neighbors only in the first
P dimensions of the Q-dimensional array F. This is useful to
process F in ``slices''.
REMARK. Matrices (2-array) with a singleton dimension are
interpreted as vectors (1-array). So for example SIFTLOCALMAX([0 1
0]) and SIFTLOCALMAX([0 1 0]') both return 2 as an aswer. However,
if [0 1 0] is to be interpreted as a 1x2 matrix, then the correct
answer is the empty set, as all elements are on the boundary.
Unfortunately MATLAB does not distinguish between vectors and
2-matrices with a singleton dimension. To forece the
interpretation of all matrices as 2-arrays, use
SIFTLOCALMAX(F,TRESH,2) (but note that in this case the result is
always empty!).
**/
#define NDIMS 3
F2D* siftlocalmax(F2D* in, float thresh, int intervals, int M, int N)
{
int ndims, ptoffset=0, maxIter = 0;
int offsets[NDIMS];
int midx[NDIMS];
int dims[NDIMS];
I2D* neighbors ;
int nneighbors ;
F2D* out;
/**
We pass dogss[o], which has >1 intervals
If >1 intervals, then the dimensions of in[F] will be 1 x intervals
If =1 interval, then the dimensions of in[F] will be the size of the dogss[o] image
We will check for this condition.
**/
ndims = NDIMS; /* Since we have intervals number of images of size M x N */
dims[0] = M;
dims[1] = N;
dims[2] = intervals-1;
/*
If there are only two dimensions and if one is singleton, then
assume that a vector has been provided as input (and treat this
as a COLUMN matrix with p=1). We do this because Matlab does not
distinguish between vectors and 1xN or Mx1 matrices and because
the cases 1xN and Mx1 are trivial (the result is alway empty).
*/
/* ------------------------------------------------------------------
* Do the job
* --------------------------------------------------------------- */
int maxima_size = M*N ;
I2D* maxima_start = iMallocHandle(1, maxima_size);
int* maxima_iterator = maxima_start->data ;
int* maxima_end = maxima_start->data + maxima_size ;
int i,h,o,j;
F2D* pt;
pt = in;
/* Compute the offsets between dimensions. */
offsets[0] = 1 ;
for(i = 1 ; i < ndims ; ++i)
{
offsets[i] = offsets[i-1]*dims[i-1] ;
}
/* Multi-index. */
for(i = 0 ; i < ndims ; ++i)
midx[i] = 1 ;
/* Neighbors. */
nneighbors = 1 ;
o=0 ;
for(i = 0 ; i < ndims ; ++i)
{
nneighbors *= 3 ;
midx[i] = -1 ;
o -= offsets[i] ;
}
nneighbors -= 1 ;
neighbors = iMallocHandle(1,nneighbors);
/* Precompute offsets from offset(-1,...,-1,0,...0) to
* offset(+1,...,+1,0,...,0). */
i = 0 ;
while(1)
{
if(o != 0)
{
asubsref(neighbors, i) = o ;
i++;
}
h = 0 ;
while( o += offsets[h], (++midx[h]) > 1 )
{
o -= 3*offsets[h] ;
midx[h] = -1 ;
if(++h >= ndims)
goto stop ;
}
}
stop: ;
/* Starts at the corner (1,1,...,1,0,0,...0) */
for(h = 0 ; h < ndims ; ++h)
{
midx[h] = 1 ;
ptoffset += offsets[h] ;
}
/* ---------------------------------------------------------------
* Loop
* ------------------------------------------------------------ */
/*
If any dimension in the first P is less than 3 elements wide
then just return the empty matrix (if we proceed without doing
anything we break the carry reporting algorithm below).
*/
for(h=0 ; h < ndims ; ++h)
if(dims[h] < 3)
goto end ;
while(1)
{
/* Propagate carry along multi index midx */
h = 0 ;
while(midx[h] >= dims[h] - 1)
{
midx[h] = 1 ;
if(++h >= ndims)
goto next_layer ;
++midx[h] ;
}
/* Scan neighbors */
{
float v = asubsref(pt, ptoffset); //*pt ;
int is_greater = (v>=thresh) ? 1 : 0;
assert(ptoffset < pt->width*pt->height);
i = 0 ;
while(is_greater && i < nneighbors)
{
if(v > asubsref(pt, ptoffset + asubsref(neighbors,i)))
is_greater *= 1;
else
is_greater *= 0;
i = i+1;
}
/* Add the local maximum */
if(is_greater)
{
/* Need more space? */
if( &(maxima_iterator[maxIter]) == maxima_end)
{
int *temp, i, j;
maxima_size += M*N ;
free(maxima_start);
maxima_start = iMallocHandle(1, maxima_size);
maxima_end = maxima_start->data + maxima_size ;
maxima_iterator = maxima_end - M*N ;
maxIter = 0;
}
maxima_iterator[maxIter] = ptoffset + 1;
maxIter++;
}
/* Go to next element */
ptoffset += 1 ;
++midx[0] ;
continue ;
next_layer: ;
if( h >= ndims )
goto end ;
while((++midx[h]) >= dims[h])
{
midx[h] = 0 ;
if(++h >= ndims)
goto end ;
}
}
}
end:;
/* Return. */
{
int i=0;
out = fMallocHandle(1, maxIter);
for(i=0; i<maxIter; i++)
asubsref(out,i) = maxima_iterator[i] ;
}
/* Release space. */
free(neighbors) ;
free(maxima_start);
return out;
}
|
