1 // This file is part of OpenMVG, an Open Multiple View Geometry C++ library.
2 
3 // Copyright (c) 2015 Pierre MOULON, Romuald PERROT.
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_SIFT_HIERARCHICAL_GAUSSIAN_SCALE_SPACE_HPP
10 #define OPENMVG_FEATURES_SIFT_HIERARCHICAL_GAUSSIAN_SCALE_SPACE_HPP
11 
12 #include <algorithm>
13 #include <vector>
14 
15 #include "openMVG/features/sift/octaver.hpp"
16 #include "openMVG/image/image_filtering.hpp"
17 #include "openMVG/image/image_resampling.hpp"
18 #include "openMVG/numeric/numeric.h"
19 #include "openMVG/system/logger.hpp"
20 
21 namespace openMVG{
22 namespace features{
23 
24 struct Octave{
25   int octave_level;                   // the octave level
26   float delta;                        // sampling rate in this octave
27   std::vector<float> sigmas;          // sigma values
28   std::vector<image::Image<float>> slices;  // octave slice (from fine to coarse)
29 
30 };
31 
32 struct GaussianScaleSpaceParams
33 {
34   float sigma_min; // minimal level of blur featured in the scale-space
35   float delta_min; // minimal inter-pixel distance featured in the scale-space
36   float sigma_in;  // assumed level of blur in the input image
37   int supplementary_levels; //necessary for post processing (i.e DOG)
38 
39   /*
40 
41   Explanation for the supplementary_levels parameters:
42    For some computation you need supplementary levels,
43     to ensure manipulation accross the gaussian pyramid.
44    i.e. SIFT: when maxima/minima are searched over three levels
45 
46 
47                        Slice 0 => ex. 1 slice per octave was asked
48                      /
49                DOG 0
50              /       \
51            /           Slice 1
52          /           /
53  MAX 0   ----  DOG 1
54          \           \
55            \           Slice 2
56              \       /
57                DOG 2
58                      \
59                        Slice 3 => 3 supplementary level are required to produce 3 DoGs
60 
61     to resume : Max computation required n+2 Dog so it implies n+3 Slice
62   */
63 
GaussianScaleSpaceParamsopenMVG::features::GaussianScaleSpaceParams64   GaussianScaleSpaceParams(
65     const float sigma_min_ = 1.6f,
66     const float delta_min_ = 1.0f,
67     const float sigma_in_  = 0.5f,
68     const int supplementary_levels_ = 0
69     ):
70     sigma_min(sigma_min_),
71     delta_min(delta_min_),
72     sigma_in(sigma_in_),
73     supplementary_levels(supplementary_levels_)
74   {  }
75 };
76 
77 struct HierarchicalGaussianScaleSpace: public Octaver<Octave>
78 {
79   /**
80   * @brief HierarchicalGaussianScaleSpace constructor
81   * @param nb_octave Number of Octave
82   * @param nb_slice Number of slice per octave
83   * @param params Parameters of the Gaussian scale space
84   */
HierarchicalGaussianScaleSpaceopenMVG::features::HierarchicalGaussianScaleSpace85   HierarchicalGaussianScaleSpace(
86     const int nb_octave = 6,
87     const int nb_slice = 3,
88     const GaussianScaleSpaceParams & params =
89       std::move(GaussianScaleSpaceParams())
90   ) :Octaver<Octave>(nb_octave, nb_slice),
91     m_params(params),
92     m_cur_octave_id(0)
93   {
94   }
95 
96   /**
97   * @brief Set Initial image and update nb_octave if necessary
98   * @param img Input image
99   */
SetImageopenMVG::features::HierarchicalGaussianScaleSpace100   virtual void SetImage(const image::Image<float> & img)
101   {
102     const double sigma_extra =
103       sqrt(Square(m_params.sigma_min) - Square(m_params.sigma_in)) / m_params.delta_min;
104     if (m_params.delta_min == 1.0f)
105     {
106       image::ImageGaussianFilter(img, sigma_extra, m_cur_base_octave_image);
107     }
108     else  // delta_min == 1
109     {
110       if (m_params.delta_min == 0.5f)
111       {
112         image::Image<float> tmp;
113         ImageUpsample(img, tmp);
114         image::ImageGaussianFilter(tmp, sigma_extra, m_cur_base_octave_image);
115       }
116       else
117       {
118         OPENMVG_LOG_ERROR
119           << "Upsampling or downsampling with delta equal to: "
120           << m_params.delta_min << " is not yet implemented";
121       }
122     }
123     //-- Limit the size of the last octave to be at least 32x32 pixels
124     const int nbOctaveMax = std::ceil(std::log2( std::min(m_cur_base_octave_image.Width(), m_cur_base_octave_image.Height())/32));
125     m_nb_octave = std::min(m_nb_octave, nbOctaveMax);
126   }
127 
128   /**
129   * @brief Compute a full octave
130   * @param[out] oct Computed octave
131   * @retval true If an octave have been computed
132   * @retval false if no octave have been computed (process ended)
133   */
NextOctaveopenMVG::features::HierarchicalGaussianScaleSpace134   virtual bool NextOctave(Octave & octave)
135   {
136     if (m_cur_octave_id >= m_nb_octave)
137     {
138       return false;
139     }
140     else
141     {
142       octave.octave_level = m_cur_octave_id;
143       if (m_cur_octave_id == 0)
144       {
145         octave.delta = m_params.delta_min;
146       }
147       else
148       {
149         octave.delta *= 2.0f;
150       }
151 
152       // init the "blur"/sigma scale spaces values
153       octave.slices.resize(m_nb_slice + m_params.supplementary_levels);
154       octave.sigmas.resize(m_nb_slice + m_params.supplementary_levels);
155       for (int s = 0; s < m_nb_slice  + m_params.supplementary_levels; ++s)
156       {
157         octave.sigmas[s] =
158           octave.delta / m_params.delta_min * m_params.sigma_min * pow(2.0,(float)s/(float)m_nb_slice);
159       }
160 
161       // Build the octave iteratively
162       octave.slices[0] = m_cur_base_octave_image;
163       for (int s = 1; s < octave.sigmas.size(); ++s)
164       {
165         // Iterative blurring the previous image
166         const image::Image<float> & im_prev = octave.slices[s-1];
167         image::Image<float> & im_next = octave.slices[s];
168         const double sig_prev = octave.sigmas[s-1];
169         const double sig_next = octave.sigmas[s];
170         const double sigma_extra = sqrt(Square(sig_next) - Square(sig_prev)) / octave.delta;
171 
172         image::ImageGaussianFilter(im_prev, sigma_extra, im_next);
173       }
174       /*
175       // Debug: Export DoG scale space on disk
176       for (int s = 0; s < octave.sigmas.size(); ++s)
177       {
178         std::stringstream os;
179         os << "DoG_out_00" << m_cur_octave_id << "_s" << "00" << s << ".png";
180         image::WriteImage(os.str().c_str(), octave.slices[s]);
181       }
182       */
183 
184       // Prepare for next octave computation -> Decimation
185       ++m_cur_octave_id;
186       if (m_cur_octave_id < m_nb_octave)
187       {
188         // Decimate => sigma * 2 for the next iteration
189         const int index = (m_params.supplementary_levels == 0) ? 1 : m_params.supplementary_levels;
190         ImageDecimate(octave.slices[octave.sigmas.size()-index], m_cur_base_octave_image);
191       }
192       return true;
193     }
194   }
195 
196 protected:
197   GaussianScaleSpaceParams m_params;  // The Gaussian scale space parameters
198   image::Image<float> m_cur_base_octave_image; // The image that will be used to generate the next octave
199   int m_cur_octave_id; // The current Octave id [0 -> Octaver::m_nb_octave]
200 };
201 
202 } // namespace features
203 } // namespace openMVG
204 
205 #endif // OPENMVG_FEATURES_SIFT_HIERARCHICAL_GAUSSIAN_SCALE_SPACE_HPP
206