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