1 // This file is part of OpenMVG, an Open Multiple View Geometry C++ library.
2 
3 // Copyright (c) 2014 Pierre MOULON.
4 
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #ifndef OPENMVG_FEATURES_DIPOLE_DIPOLE_DESCRIPTOR_HPP
10 #define OPENMVG_FEATURES_DIPOLE_DIPOLE_DESCRIPTOR_HPP
11 
12 
13 #include "openMVG/features/feature.hpp"
14 #include "openMVG/features/descriptor.hpp"
15 #include "openMVG/image/image_container.hpp"
16 #include "openMVG/image/sample.hpp"
17 
18 //------------------
19 //-- Bibliography --
20 //------------------
21 //- [1] "New local descriptor based on dissociated dipoles"
22 //- Authors: Alexis Joly.
23 //- Date: December 2007.
24 //- Conference: CIVR.
25 
26 namespace openMVG
27 {
28 namespace features
29 {
30   // Create a DIPOLE describer [1].
31   //
32   // Note :
33   // - Angle is in radians.
34   // - data the output array (must be allocated to 20 values).
35   template<typename Real>
PickNaiveDipole(const image::Image<Real> & image,float x,float y,float scale,float angle,float * data)36   void PickNaiveDipole
37   (
38     const image::Image<Real> & image,
39     float x,
40     float y,
41     float scale,
42     float angle,
43     float * data
44   )
45   {
46     // Use bilinear sampling
47     const image::Sampler2d<image::SamplerLinear> sampler;
48     // Setup the rotation center.
49     const float & cx = x, & cy = y;
50 
51     const float lambda1 = scale;
52     const float lambda2 = lambda1 / 2.0f;
53     static const float angleSubdiv = 2.0f * M_PI / 12.0f;
54     const float memoizeCos[12] =
55     {
56       std::cos(angle + 0.f * angleSubdiv),
57       std::cos(angle + 1.f * angleSubdiv),
58       std::cos(angle + 2.f * angleSubdiv),
59       std::cos(angle + 3.f * angleSubdiv),
60       std::cos(angle + 4.f * angleSubdiv),
61       std::cos(angle + 5.f * angleSubdiv),
62       std::cos(angle + 6.f * angleSubdiv),
63       std::cos(angle + 7.f * angleSubdiv),
64       std::cos(angle + 8.f * angleSubdiv),
65       std::cos(angle + 9.f * angleSubdiv),
66       std::cos(angle + 10.f * angleSubdiv),
67       std::cos(angle + 11.f * angleSubdiv)
68     };
69     const float memoizeSin[12] =
70     {
71       std::sin(angle + 0.f * angleSubdiv),
72       std::sin(angle + 1.f * angleSubdiv),
73       std::sin(angle + 2.f * angleSubdiv),
74       std::sin(angle + 3.f * angleSubdiv),
75       std::sin(angle + 4.f * angleSubdiv),
76       std::sin(angle + 5.f * angleSubdiv),
77       std::sin(angle + 6.f * angleSubdiv),
78       std::sin(angle + 7.f * angleSubdiv),
79       std::sin(angle + 8.f * angleSubdiv),
80       std::sin(angle + 9.f * angleSubdiv),
81       std::sin(angle + 10.f * angleSubdiv),
82       std::sin(angle + 11.f * angleSubdiv)
83     };
84 
85     Vecf dipoleF1(12);
86     for (int i = 0; i < 12; ++i)
87     {
88       const float xi = cx + lambda1 * memoizeCos[i];
89       const float yi = cy + lambda1 * memoizeSin[i];
90       dipoleF1(i) = sampler(image, yi, xi);
91     }
92     Matf A(8,12);
93     A <<  0, 0, 0, 1, 0, 0, 0, 0, 0,-1, 0, 0,
94           0,-1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
95           0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 1,
96           0, 0, 0, 0, 1, 0, 0,-1, 0, 0, 0, 0,
97           0, 0, 0, 0, 0, 0, 1, 0, 0,-1, 0, 0,
98           0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,-1,
99           0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
100           1, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0;
101 
102     // Add the second order F2 dipole
103     Vecf dipoleF2(12);
104     for (int i = 0; i < 12; ++i)
105     {
106       const float xi = cx + (lambda1 + lambda2) * memoizeCos[i];
107       const float yi = cy + (lambda1 + lambda2) * memoizeSin[i];
108 
109       const float xii = cx + (lambda1 - lambda2) * memoizeCos[i];
110       const float yii = cy + (lambda1 - lambda2) * memoizeSin[i];
111 
112       // Bilinear interpolation
113       dipoleF2(i) = sampler(image, yi, xi) - sampler(image, yii, xii);
114     }
115     // Normalize to be affine luminance invariant (a*I(x,y)+b).
116     Map<Vecf> dataMap( data, 20);
117     dataMap.block<8,1>(0,0) = (A * dipoleF1).normalized();
118     dataMap.block<12,1>(8,0) = dipoleF2.normalized();
119   }
120 
121   // Pick an angular smoothed dipole
122   template<typename Real>
PickASDipole(const image::Image<Real> & image,float x,float y,float scale,float angle,float * data)123   void PickASDipole
124   (
125     const image::Image<Real> & image,
126     float x,
127     float y,
128     float scale,
129     float angle,
130     float * data)
131   {
132     const image::Sampler2d<image::SamplerLinear> sampler;
133     // Setup the rotation center.
134     const float & cx = x, & cy = y;
135 
136     const float lambda1 = scale;
137     const float lambda2 = lambda1 / 2.0f;
138     const float angleSubdiv = 2.0f * M_PI / 12.0f;
139 
140     //-- First order dipole:
141     Vecf dipoleF1(12);
142     for (int i = 0; i < 12; ++i)
143     {
144       const float xi = cx + lambda1 * std::cos(angle + i * angleSubdiv);
145       const float yi = cy + lambda1 * std::sin(angle + i * angleSubdiv);
146       const float xi0 = cx + lambda1 * std::cos(angle + i * angleSubdiv - angleSubdiv/2.0);
147       const float yi0 = cy + lambda1 * std::sin(angle + i * angleSubdiv - angleSubdiv/2.0);
148       const float xi3 = cx + lambda1 * std::cos(angle + i * angleSubdiv + angleSubdiv/2.0);
149       const float yi3 = cy + lambda1 * std::sin(angle + i * angleSubdiv + angleSubdiv/2.0);
150       // Bilinear interpolation
151       dipoleF1(i) =
152         (sampler(image, yi, xi) +
153          sampler(image, yi0, xi0) +
154          sampler(image, yi3, xi3))/3.0f;
155     }
156     Matf A(8,12);
157     A <<  0, 0, 0, 1, 0, 0, 0, 0, 0,-1, 0, 0,
158           0,-1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
159           0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 1,
160           0, 0, 0, 0, 1, 0, 0,-1, 0, 0, 0, 0,
161           0, 0, 0, 0, 0, 0, 1, 0, 0,-1, 0, 0,
162           0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,-1,
163           0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
164           1, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0;
165 
166     //-- Second order dipole:
167     Vecf dipoleF2(12);
168     for (int i = 0; i < 12; ++i)
169     {
170       const double angleSample = i * angleSubdiv;
171       const float xi = cx + (lambda1 + lambda2) * std::cos(angle + angleSample);
172       const float yi = cy + (lambda1 + lambda2) * std::sin(angle + angleSample);
173       const float xi0 = cx + (lambda1 + lambda2) * std::cos(angle + angleSample - angleSubdiv/2.0);
174       const float yi0 = cy + (lambda1 + lambda2) * std::sin(angle + angleSample - angleSubdiv/2.0);
175       const float xi3 = cx + (lambda1 + lambda2) * std::cos(angle + angleSample + angleSubdiv/2.0);
176       const float yi3 = cy + (lambda1 + lambda2) * std::sin(angle + angleSample + angleSubdiv/2.0);
177 
178       const float xii = cx + (lambda1 - lambda2) * std::cos(angle + angleSample);
179       const float yii = cy + (lambda1 - lambda2) * std::sin(angle + angleSample);
180       const float xii0 = cx + (lambda1 - lambda2) * std::cos(angle + angleSample - angleSubdiv/2.0);
181       const float yii0 = cy + (lambda1 - lambda2) * std::sin(angle + angleSample - angleSubdiv/2.0);
182       const float xii3 = cx + (lambda1 - lambda2) * std::cos(angle + angleSample + angleSubdiv/2.0);
183       const float yii3 = cy + (lambda1 - lambda2) * std::sin(angle + angleSample + angleSubdiv/2.0);
184 
185       // Bilinear interpolation
186       dipoleF2(i) =
187         (sampler(image, yi, xi)   +
188          sampler(image, yi0, xi0) +
189          sampler(image, yi3, xi3)) /3.0f
190         -
191        (sampler(image, yii, xii)   +
192         sampler(image, yii0, xii0) +
193         sampler(image, yii3, xii3)) /3.0f;
194     }
195     // Normalize to be affine luminance invariant (a*I(x,y)+b).
196     Map<Vecf> dataMap( data, 20);
197     dataMap.block<8,1>(0,0) = (A * dipoleF1).normalized();
198     dataMap.block<12,1>(8,0) = dipoleF2.normalized();
199   }
200 
201    /**
202     ** @brief Compute DIPOLE descriptor for a given interest point
203     ** @param Li Input image
204     ** @param ipt Input interest point
205     ** @param desc output descriptor (floating point descriptor)
206     ** @param bAngularSmoothedDipole Tell if we must extract an upright or an
207     **  angular smoothed dipole
208     ** @param magnif_factor Scaling factor used to rescale the dipole sampling
209     **/
210   template<typename Real>
ComputeDipoleDescriptor(const image::Image<Real> & Li,const SIOPointFeature & ipt,Descriptor<float,20> & desc,bool bAngularSmoothedDipole=true,const float magnif_factor=3.5f)211   void ComputeDipoleDescriptor
212   (
213     const image::Image<Real> & Li,
214     const SIOPointFeature & ipt,
215     Descriptor<float, 20> & desc,
216     bool bAngularSmoothedDipole = true,
217     const float magnif_factor = 3.5f
218   )
219   {
220     if (bAngularSmoothedDipole)
221       PickASDipole(Li, ipt.x(), ipt.y(), ipt.scale() * magnif_factor, ipt.orientation(), &desc[0]);
222     else
223       PickNaiveDipole(Li, ipt.x(), ipt.y(), ipt.scale() * magnif_factor, ipt.orientation(), &desc[0]);
224   }
225 
226 } // namespace features
227 } // namespace openMVG
228 
229 #endif // OPENMVG_FEATURES_DIPOLE_DIPOLE_DESCRIPTOR_HPP
230