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
|
/********************************
Author: Sravanthi Kota Venkata
********************************/
#include "sift.h"
F2D* resizeArray(F2D* array, int omin)
{
F2D* prev = NULL;
F2D* current = array;
int o;
if(omin<0)
{
for(o=1; o>=-omin; o--)
{
prev = current;
current = doubleSize(current);
fFreeHandle(prev);
}
}
if(omin>0)
{
for(o=1; o<= omin; o++)
{
prev = current;
current = halveSize(current);
fFreeHandle(prev);
}
}
return current;
}
/**
Returns the Gaussian scale space of image I. Image I is assumed to be
pre-smoothed at level SIGMAN. O,S,OMIN,SMIN,SMAX,SIGMA0 are the
parameters of the scale space.
**/
F2D** gaussianss(F2D* array, float sigman, int O, int S, int omin, int smin, int smax, float sigma0)
{
/* We compute the following items in the function
1. Smooth input image per octave
2. Smooth each octave for different intervals
3. Subtract each "interval-1" smooth image from "interval" image per octave. So, per octave, we have "interval" * DOG images.
4. So, octave * intervals * image
5. Note: At each octave, the image is scaled down in both x and y directions
*/
float k, dsigma0, dsigma;
int s, i, j, o, so, M, N, sbest;
int intervals = smax-smin+1;
float temp, target_sigma, prev_sigma;
F2D *TMP, **gss;
F2D* I = array;
// Scale multiplicative step
k = pow(2, (1.0/S));
dsigma0 = sigma0 * sqrt(1-(1.0/pow(k,2))); // Scale step factor
// If omin < 0, multiply the size of the image.
I = resizeArray(I, omin);
M = I->height;
N = I->width;
so = -smin+1; // Index offset
gss = malloc(O*intervals*sizeof(F2D*));
if(gss == NULL)
{
printf("Could not allocate memory\n");
return NULL;
}
/**
First octave
--------------------------------------------------------------------
The first level of the first octave has scale index (o,s) =
(omin,smin) and scale coordinate
sigma(omin,smin) = sigma0 2^omin k^smin
The input image I is at nominal scale sigman. Thus in order to get
the first level of the pyramid we need to apply a smoothing of
sqrt( (sigma0 2^omin k^smin)^2 - sigman^2 ).
As we have pre-scaled the image omin octaves (up or down,
depending on the sign of omin), we need to correct this value
by dividing by 2^omin, getting
sqrt( (sigma0 k^smin)^2 - (sigman/2^omin)^2 )
**/
temp = sqrt(pow((sigma0*pow(k,smin)),2) - pow((sigman/pow(2,omin)),2));
{
gss[0] = fSetArray(I->height, I->width, 0);
imsmooth(I, temp, gss[0] );
}
for(s=smin; s<smax; s++)
{
/**
Here we go from (omin,s-1) to (omin,s). The extra smoothing
standard deviation is
(sigma0 2^omin 2^(s/S) )^2 - (simga0 2^omin 2^(s/S-1/S) )^2
After dividing by 2^omin (to take into account the fact
that the image has been pre-scaled omin octaves), the
standard deviation of the smoothing kernel is
dsigma = sigma0 k^s sqrt(1-1/k^2)
**/
dsigma = pow(k,s+1) * dsigma0;
gss[s+so] = fSetArray(gss[s+so-1]->height, gss[s+so-1]->width, 0);
imsmooth( gss[(s+so-1)] , dsigma, gss[(s+so)] );
}
/** Other octaves **/
for(o=1; o<O; o++)
{
/**
We need to initialize the first level of octave (o,smin) from
the closest possible level of the previous octave. A level (o,s)
in this octave corrsponds to the level (o-1,s+S) in the previous
octave. In particular, the level (o,smin) correspnds to
(o-1,smin+S). However (o-1,smin+S) might not be among the levels
(o-1,smin), ..., (o-1,smax) that we have previously computed.
The closest pick is
/ smin+S if smin+S <= smax
(o-1,sbest) , sbest = |
\ smax if smin+S > smax
The amount of extra smoothing we need to apply is then given by
( sigma0 2^o 2^(smin/S) )^2 - ( sigma0 2^o 2^(sbest/S - 1) )^2
As usual, we divide by 2^o to cancel out the effect of the
downsampling and we get
( sigma 0 k^smin )^2 - ( sigma0 2^o k^(sbest - S) )^2
**/
sbest = MIN(smin+S-1, smax-1);
TMP = halveSize( gss[(o-1)*intervals+sbest+so]);
target_sigma = sigma0 * pow(k,smin);
prev_sigma = sigma0 * pow(k, (sbest+1)-S);
temp = sqrt(pow(target_sigma,2) - pow(prev_sigma, 2));
if(target_sigma > prev_sigma)
{
gss[o*intervals] = fSetArray(TMP->height, TMP->width, 0);
imsmooth(TMP, temp, gss[o*intervals] );
}
else
{
int i;
gss[o*intervals] = fSetArray(TMP->height, TMP->width, 0);
for(i=0; i<(TMP->height*TMP->width); i++)
asubsref(gss[o*intervals],i) = asubsref(TMP,i);
}
M = TMP->height;
N = TMP->width;
fFreeHandle(TMP);
for(s=smin; s<smax; s++)
{
// The other levels are determined as above for the first octave.
dsigma = pow(k,s+1) * dsigma0;
gss[o*intervals+s+so] = fSetArray(gss[o*intervals+s-1+so]->height, gss[o*intervals+s-1+so]->width, 0);
imsmooth( gss[o*intervals+s-1+so] , dsigma, gss[o*intervals+s+so]);
}
}
fFreeHandle(I);
return gss;
}
|
