/** @file scalespace.c ** @brief Scale Space - Definition ** @author Karel Lenc ** @author Andrea Vedaldi ** @author Michal Perdoch **/ /* Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson. All rights reserved. This file is part of the VLFeat library and is made available under the terms of the BSD license (see the COPYING file). */ /** @page scalespace Gaussian Scale Space (GSS) @author Karel Lenc @author Andrea Vedaldi @author Michal Perdoch @tableofcontents @ref scalespace.h implements a Gaussian scale space, a data structure representing an image at multiple resolutions @cite{witkin83scale-space} @cite{koenderink84the-structure} @cite{lindeberg94scale-space}. Scale spaces have many use, including the detection of co-variant local features @cite{lindeberg98principles} such as SIFT, Hessian-Affine, Harris-Affine, Harris-Laplace, etc. @ref scalespace-starting demonstreates how to use the C API to compute the scalespace of an image. For further details refer to: - @subpage scalespace-fundamentals @section scalespace-starting Getting started Given an input image `image`, the following example uses the ::VlScaleSpace object to compute its Gaussian scale space and return the image `level` at scale `(o,s)`, where `o` is the octave and `s` is the octave subdivision or sublevel: @code float* level ; VlScaleSpace ss = vl_scalespace_new(imageWidth, imageHeight) ; vl_scalespace_put_image(ss, image) ; level = vl_scalespace_get_level(ss, o, s) ; @endcode The image `level` is obtained by convolving `image` by a Gaussian filter of isotropic standard deviation given by @code double sigma = vl_scalespace_get_sigma(ss, o, s) ; @endcode The resolution of `level` is in general different from the resolution of `image` and is determined by the octave `o`. It can be obtained as follows: @code VlScaleSpaceOctaveGeometry ogeom = vl_scalespace_get_octave_geometry(ss, o) ; ogeom.width // width of level (in number of pixels) ogeom.height // height of level (in number of pixels) ogeom.step // spatial sampling step @endcode The parameter `ogeom.step` is the sampling step relatively to the sampling of the input image `image`. The ranges of valid octaves and scale sublevels can be obtained as @code VlScaleSpaceGeometry geom = vl_scalespace_get_geometry(ss) ; geom.firstOctave // Index of the fisrt octave geom.lastOctave // Index of the last octave geom.octaveResolution ; // Number of octave subdivisions geom.octaveFirstSubdivision // Index of the first octave subdivision geom.octaveLastSubdivision // Index of the last octave subdivision @endcode So for example `o` minimum value is `geom.firstOctave` and maximum value is `geom.lastOctave`. The subdivision index `s` naturally spans the range 0 to `geom.octaveResolution-1`. However, the scale space object is flexible in that it allows different ranges of subdivisions to be computed and `s` varies in the range `geom.octaveFirstSubdivision` to `geom.octaveLastSubdivision`. See @ref scalespace-fundamentals for further details. The geometry of the scale space can be customized upon creation, as follows: @code VlScaleSpaceGeometry geom = vl_scalespace_get_default_geometry(imageWidth, imageHeight) ; geom.firstOctave = -1 ; geom.octaveFirstSubdivision = -1 ; geom.octaveLastSubdivision = geom.octaveResolution ; VlScaleSpacae ss = vl_scalespace_new_with_geometry (geom) ; @endcode @page scalespace-fundamentals Gaussian scale space fundamentals @tableofcontents This page discusses the notion of *Gaussian scale space* and the relative data structure. For the C API see @ref scalespace.h and @ref scalespace-starting. A *scale space* is representation of an image at multiple resolution levels. An image is a function $\ell(x,y)$ of two coordinates $x$, $y$; the scale space $\ell(x,y,\sigma)$ adds a third coordinate $\sigma$ indexing the *scale*. Here the focus is the Gaussian scale space, where the image $\ell(x,y,\sigma)$ is obtained by smoothing $\ell(x,y)$ by a Gaussian kernel of isotropic standard deviation $\sigma$. @section scalespace-definition Scale space definition Formally, the *Gaussian scale space* of an image $\ell(x,y)$ is defined as \[ \ell(x,y,\sigma) = [g_{\sigma} * \ell](x,y,\sigma) \] where $g_\sigma$ denotes a 2D Gaussian kernel of isotropic standard deviation $\sigma$: \[ g_{\sigma}(x,y) = \frac{1}{2\pi\sigma^2} \exp\left( - \frac{x^2 + y^2}{2\sigma^2} \right). \] An important detail is that the algorithm computing the scale space assumes that the input image $\ell(x,y)$ is pre-smoothed, roughly capturing the effect of the finite pixel size in a CCD. This is modelled by assuming that the input is not $\ell(x,y)$, but $\ell(x,y,\sigma_n)$, where $\sigma_n$ is a *nominal smoothing*, usually taken to be 0.5 (half a pixel standard deviation). This also means that $\sigma = \sigma_n = 0.5$ is the *finest scale* that can actually be computed. The scale space structure stores samples of the function $\ell(x,y,\sigma)$. The density of the sampling of the spatial coordinates $x$ and $y$ is adjusted as a function of the scale $\sigma$, corresponding to the intuition that images at a coarse resolution can be sampled more coarsely without loss of information. Thus, the scale space has the structure of a *pyramid*: a collection of digital images sampled at progressively coarser spatial resolution and hence of progressively smaller size (in pixels). The following figure illustrates the scale space pyramid structure: @image html scalespace-basic.png "A scalespace structure with 2 octaves and S=3 subdivisions per octave" The pyramid is organised in a number of *octaves*, indexed by a parameter `o`. Each octave is further subdivided into *sublevels*, indexed by a parameter `s`. These are related to the scale $\sigma$ by the equation \[ \sigma(s,o) = \sigma_o 2^{\displaystyle o + \frac{s}{\mathtt{octaveResolution}}} \] where `octaveResolution` is the resolution of the octave subsampling $\sigma_0$ is the *base smoothing*. At each octave the spatial resolution is doubled, in the sense that samples are take with a step of \[ \mathtt{step} = 2^o. \] Hence, denoting as `level[i,j]` the corresponding samples, one has $\ell(x,y,\sigma) = \mathtt{level}[i,j]$, where \[ (x,y) = (i,j) \times \mathtt{step}, \quad \sigma = \sigma(o,s), \quad 0 \leq i < \mathtt{lwidth}, \quad 0 \leq j < \mathtt{lheight}, \] where \[ \mathtt{lwidth} = \lfloor \frac{\mathtt{width}}{2^\mathtt{o}}\rfloor, \quad \mathtt{lheight} = \lfloor \frac{\mathtt{height}}{2^\mathtt{o}}\rfloor. \] @section scalespace-geometry Scale space geometry In addition to the parameters discussed above, the geometry of the data stored in a scale space structure depends on the range of allowable octaves `o` and scale sublevels `s`. While `o` may range in any reasonable value given the size of the input image `image`, usually its minimum value is either 0 or -1. The latter corresponds to doubling the resolution of the image in the first octave of the scale space and it is often used in feature extraction. While there is no information added to the image by upsampling in this manner, fine scale filters, including derivative filters, are much easier to compute by upsalmpling first. The maximum practical value is dictated by the image resolution, as it should be $2^o\leq\min\{\mathtt{width},\mathtt{height}\}$. VLFeat has the flexibility of specifying the range of `o` using the `firstOctave` and `lastOctave` parameters of the ::VlScaleSpaceGeometry structure. The sublevel `s` varies naturally in the range $\{0,\dots,\mathtt{octaveResolution}-1\}$. However, it is often convenient to store a few extra levels per octave (e.g. to compute the local maxima of a function in scale or the Difference of Gaussian cornerness measure). Thus VLFeat scale space structure allows this parameter to vary in an arbitrary range, specified by the parameters `octaveFirstSubdivision` and `octaveLastSubdivision` of ::VlScaleSpaceGeometry. Overall the possible values of the indexes `o` and `s` are: \[ \mathtt{firstOctave} \leq o \leq \mathtt{lastOctave}, \qquad \mathtt{octaveFirstSubdivision} \leq s \leq \mathtt{octaveLastSubdivision}. \] Note that, depending on these ranges, there could be *redundant pairs* of indexes `o` and `s` that represent the *same* pyramid level at more than one sampling resolution. In practice, the ability to generate such redundant information is very useful in algorithms using scalespaces, as coding multiscale operations using a fixed sampling resolution is far easier. For example, the DoG feature detector computes the scalespace with three redundant levels per octave, as follows: @image html scalespace.png "A scalespace containing redundant representation of certain scale levels." @section scalespace-algorithm Algorithm and limitations Given $\ell(x,y,\sigma_n)$, any of a vast number digitial filtering techniques can be used to compute the scale levels. Presently, VLFeat uses a basic FIR implementation of the Gaussian filters. The FIR implementation is obtained by sampling the Gaussian function and re-normalizing it to have unit norm. This simple construction does not account properly for sampling effect, which may be a problem for very small Gausisan kernels. As a rule of thumb, such filters work sufficiently well for, say, standard deviation $\sigma$ at least 1.6 times the sampling step. A work around to apply this basic FIR implementation to very small Gaussian filters is to upsample the image first. The limitations on the FIR filters have relatively important for the pyramid construction, as the latter is obtained by *incremental smoothing*: each successive level is obtained from the previous one by adding the needed amount of smoothing. In this manner, the size of the FIR filters remains small, which makes them efficient; at the same time, for what discussed, excessively small filters are not represented properly. */ #include "scalespace.h" #include "mathop.h" #include #include #include #include #include /** @file scalespace.h ** @struct VlScaleSpace ** @brief Scale space class ** ** This is an opaque class used to compute the scale space of an ** image. **/ struct _VlScaleSpace { VlScaleSpaceGeometry geom ; /**< Geometry of the scale space */ float **octaves ; /**< Data */ } ; /* ---------------------------------------------------------------- */ /** @brief Get the default geometry for a given image size. ** @param width image width. ** @param height image height. ** @return the default scale space geometry. ** ** Both @a width and @a height must be at least one pixel wide. **/ VlScaleSpaceGeometry vl_scalespace_get_default_geometry (vl_size width, vl_size height) { VlScaleSpaceGeometry geom ; assert(width >= 1) ; assert(height >= 1) ; geom.width = width ; geom.height = height ; geom.firstOctave = 0 ; geom.lastOctave = VL_MAX(floor(vl_log2_d(VL_MIN(width, height))) - 3, 0) ; geom.octaveResolution= 3 ; geom.octaveFirstSubdivision = 0 ; geom.octaveLastSubdivision = geom.octaveResolution - 1 ; geom.baseScale = 1.6 * pow(2.0, 1.0 / geom.octaveResolution) ; geom.nominalScale = 0.5 ; return geom ; } #define is_valid_geometry(geom) (\ geom.firstOctave <= geom.lastOctave && \ geom.octaveResolution >= 1 && \ geom.octaveFirstSubdivision <= geom.octaveLastSubdivision && \ geom.baseScale >= 0.0 && \ geom.nominalScale >= 0.0) /** @brief Check scale space geometries for equality ** @param a first geometry. ** @param b second geometry. ** @return true if equal. **/ vl_bool vl_scalespacegeometry_is_equal (VlScaleSpaceGeometry a, VlScaleSpaceGeometry b) { return a.width == b.width && a.height == b.height && a.firstOctave == b.firstOctave && a.lastOctave == b.lastOctave && a.octaveResolution == b.octaveResolution && a.octaveFirstSubdivision == b.octaveLastSubdivision && a.baseScale == b.baseScale && a.nominalScale == b.nominalScale ; } /** @brief Get the geometry of the scale space. ** @param self object. ** @return the scale space geometry. **/ VlScaleSpaceGeometry vl_scalespace_get_geometry (VlScaleSpace const * self) { return self->geom ; } /** @brief Get the geometry of an octave of the scalespace. ** @param self object. ** @param o octave index. ** @return the geometry of octave @a o. **/ VlScaleSpaceOctaveGeometry vl_scalespace_get_octave_geometry (VlScaleSpace const * self, vl_index o) { VlScaleSpaceOctaveGeometry ogeom ; ogeom.width = VL_SHIFT_LEFT(self->geom.width, -o) ; ogeom.height = VL_SHIFT_LEFT(self->geom.height, -o) ; ogeom.step = pow(2.0, o) ; return ogeom ; } /** @brief Get the data of a scale space level ** @param self object. ** @param o octave index. ** @param s level index. ** @return pointer to the data for octave @a o, level @a s. ** ** The octave index @a o must be in the range @c firstOctave ** to @c lastOctave and the scale index @a s must be in the ** range @c octaveFirstSubdivision to @c octaveLastSubdivision. **/ float * vl_scalespace_get_level (VlScaleSpace *self, vl_index o, vl_index s) { VlScaleSpaceOctaveGeometry ogeom = vl_scalespace_get_octave_geometry(self,o) ; float * octave ; assert(self) ; assert(o >= self->geom.firstOctave) ; assert(o <= self->geom.lastOctave) ; assert(s >= self->geom.octaveFirstSubdivision) ; assert(s <= self->geom.octaveLastSubdivision) ; octave = self->octaves[o - self->geom.firstOctave] ; return octave + ogeom.width * ogeom.height * (s - self->geom.octaveFirstSubdivision) ; } /** @brief Get the data of a scale space level (const) ** @param self object. ** @param o octave index. ** @param s level index. ** @return pointer to the data for octave @a o, level @a s. ** ** This function is the same as ::vl_scalespace_get_level but reutrns ** a @c const pointer to the data. **/ float const * vl_scalespace_get_level_const (VlScaleSpace const * self, vl_index o, vl_index s) { return vl_scalespace_get_level((VlScaleSpace*)self, o, s) ; } /** ------------------------------------------------------------------ ** @brief Get the scale of a given octave and sublevel ** @param self object. ** @param o octave index. ** @param s sublevel index. ** ** The function returns the scale $\sigma(o,s)$ as a function of the ** octave index @a o and sublevel @a s. **/ double vl_scalespace_get_level_sigma (VlScaleSpace const *self, vl_index o, vl_index s) { return self->geom.baseScale * pow(2.0, o + (double) s / self->geom.octaveResolution) ; } /** ------------------------------------------------------------------ ** @internal @brief Upsample the rows and take the transpose ** @param destination output image. ** @param source input image. ** @param width input image width. ** @param height input image height. ** ** The output image has dimensions @a height by 2 @a width (so the ** destination buffer must be at least as big as two times the ** input buffer). ** ** Upsampling is performed by linear interpolation. **/ static void copy_and_upsample (float *destination, float const *source, vl_size width, vl_size height) { vl_index x, y, ox, oy ; float v00, v10, v01, v11 ; assert(destination) ; assert(source) ; for(y = 0 ; y < (signed)height ; ++y) { oy = (y < ((signed)height - 1)) * width ; v10 = source[0] ; v11 = source[oy] ; for(x = 0 ; x < (signed)width ; ++x) { ox = x < ((signed)width - 1) ; v00 = v10 ; v01 = v11 ; v10 = source[ox] ; v11 = source[ox + oy] ; destination[0] = v00 ; destination[1] = 0.5f * (v00 + v10) ; destination[2*width] = 0.5f * (v00 + v01) ; destination[2*width+1] = 0.25f * (v00 + v01 + v10 + v11) ; destination += 2 ; source ++; } destination += 2*width ; } } /** ------------------------------------------------------------------ ** @internal @brief Downsample ** @param destination output imgae buffer. ** @param source input image buffer. ** @param width input image width. ** @param height input image height. ** @param numOctaves octaves (non negative). ** ** The function downsamples the image @a d times, reducing it to @c ** 1/2^d of its original size. The parameters @a width and @a height ** are the size of the input image. The destination image @a dst is ** assumed to be floor(width/2^d) pixels wide and ** floor(height/2^d) pixels high. **/ static void copy_and_downsample (float *destination, float const *source, vl_size width, vl_size height, vl_size numOctaves) { vl_index x, y ; vl_size step = 1 << numOctaves ; /* step = 2^numOctaves */ assert(destination) ; assert(source) ; if (numOctaves == 0) { memcpy(destination, source, sizeof(float) * width * height) ; } else { for(y = 0 ; y < (signed)height ; y += step) { float const *p = source + y * width ; for(x = 0 ; x < (signed)width - ((signed)step - 1) ; x += step) { *destination++ = *p ; p += step ; } } } } /* ---------------------------------------------------------------- */ /** @brief Create a new scale space object ** @param width image width. ** @param height image height. ** @return new scale space object. ** ** This function is the same as ::vl_scalespace_new_with_geometry() ** but it uses ::vl_scalespace_get_default_geometry to initialise ** the geometry of the scale space from the image size. ** ** @sa ::vl_scalespace_new_with_geometry(), ::vl_scalespace_delete(). **/ VlScaleSpace * vl_scalespace_new (vl_size width, vl_size height) { VlScaleSpaceGeometry geom ; geom = vl_scalespace_get_default_geometry(width, height) ; return vl_scalespace_new_with_geometry(geom) ; } /** ------------------------------------------------------------------ ** @brief Create a new scale space with the specified geometry ** @param geom scale space geomerty. ** @return new scale space object. ** ** If the geometry is not valid (see ::VlScaleSpaceGeometry), the ** result is unpredictable. ** ** The function returns `NULL` if it was not possible to allocate the ** object because of an out-of-memory condition. ** ** @sa ::VlScaleSpaceGeometry, ::vl_scalespace_delete(). **/ VlScaleSpace * vl_scalespace_new_with_geometry (VlScaleSpaceGeometry geom) { vl_index o ; vl_size numSublevels = geom.octaveLastSubdivision - geom.octaveFirstSubdivision + 1 ; vl_size numOctaves = geom.lastOctave - geom.firstOctave + 1 ; VlScaleSpace *self ; assert(is_valid_geometry(geom)) ; numOctaves = geom.lastOctave - geom.firstOctave + 1 ; numSublevels = geom.octaveLastSubdivision - geom.octaveFirstSubdivision + 1 ; self = vl_calloc(1, sizeof(VlScaleSpace)) ; if (self == NULL) goto err_alloc_self ; self->geom = geom ; self->octaves = vl_calloc(numOctaves, sizeof(float*)) ; if (self->octaves == NULL) goto err_alloc_octave_list ; for (o = self->geom.firstOctave ; o <= self->geom.lastOctave ; ++o) { VlScaleSpaceOctaveGeometry ogeom = vl_scalespace_get_octave_geometry(self,o) ; vl_size octaveSize = ogeom.width * ogeom.height * numSublevels ; self->octaves[o - self->geom.firstOctave] = vl_malloc(octaveSize * sizeof(float)) ; if (self->octaves[o - self->geom.firstOctave] == NULL) goto err_alloc_octaves; } return self ; err_alloc_octaves: for (o = self->geom.firstOctave ; o <= self->geom.lastOctave ; ++o) { if (self->octaves[o - self->geom.firstOctave]) { vl_free(self->octaves[o - self->geom.firstOctave]) ; } } err_alloc_octave_list: vl_free(self) ; err_alloc_self: return NULL ; } /* ---------------------------------------------------------------- */ /** @brief Create a new copy of the object ** @param self object to copy from. ** ** The function returns `NULL` if the copy cannot be made due to an ** out-of-memory condition. **/ VlScaleSpace * vl_scalespace_new_copy (VlScaleSpace* self) { vl_index o ; VlScaleSpace * copy = vl_scalespace_new_shallow_copy(self) ; if (copy == NULL) return NULL ; for (o = self->geom.firstOctave ; o <= self->geom.lastOctave ; ++o) { VlScaleSpaceOctaveGeometry ogeom = vl_scalespace_get_octave_geometry(self,o) ; vl_size numSubevels = self->geom.octaveLastSubdivision - self->geom.octaveFirstSubdivision + 1; memcpy(copy->octaves[o - self->geom.firstOctave], self->octaves[o - self->geom.firstOctave], ogeom.width * ogeom.height * numSubevels * sizeof(float)) ; } return copy ; } /* ---------------------------------------------------------------- */ /** @brief Create a new shallow copy of the object ** @param self object to copy from. ** ** The function works like ::vl_scalespace_new_copy() but only allocates ** the scale space, without actually copying the data. **/ VlScaleSpace * vl_scalespace_new_shallow_copy (VlScaleSpace* self) { return vl_scalespace_new_with_geometry (self->geom) ; } /* ---------------------------------------------------------------- */ /** @brief Delete object ** @param self object to delete. ** @sa ::vl_scalespace_new() **/ void vl_scalespace_delete (VlScaleSpace * self) { if (self) { if (self->octaves) { vl_index o ; for (o = self->geom.firstOctave ; o <= self->geom.lastOctave ; ++o) { if (self->octaves[o - self->geom.firstOctave]) { vl_free(self->octaves[o - self->geom.firstOctave]) ; } } vl_free(self->octaves) ; } vl_free(self) ; } } /* ---------------------------------------------------------------- */ /** @internal @brief Fill octave starting from the first level ** @param self object instance. ** @param o octave to process. ** ** The function takes the first sublevel of octave @a o (the one at ** sublevel `octaveFirstLevel` and iteratively ** smoothes it to obtain the other octave levels. **/ void _vl_scalespace_fill_octave (VlScaleSpace *self, vl_index o) { vl_index s ; VlScaleSpaceOctaveGeometry ogeom = vl_scalespace_get_octave_geometry(self, o) ; for(s = self->geom.octaveFirstSubdivision + 1 ; s <= self->geom.octaveLastSubdivision ; ++s) { double sigma = vl_scalespace_get_level_sigma(self, o, s) ; double previousSigma = vl_scalespace_get_level_sigma(self, o, s - 1) ; double deltaSigma = sqrtf(sigma*sigma - previousSigma*previousSigma) ; float* level = vl_scalespace_get_level (self, o, s) ; float* previous = vl_scalespace_get_level (self, o, s-1) ; vl_imsmooth_f (level, ogeom.width, previous, ogeom.width, ogeom.height, ogeom.width, deltaSigma / ogeom.step, deltaSigma / ogeom.step) ; } } /** ------------------------------------------------------------------ ** @internal @brief Initialize the first level of an octave from an image ** @param self ::VlScaleSpace object instance. ** @param image image data. ** @param o octave to start. ** ** The function initializes the first level of octave @a o from ** image @a image. The dimensions of the image are the ones set ** during the creation of the ::VlScaleSpace object instance. **/ static void _vl_scalespace_start_octave_from_image (VlScaleSpace *self, float const *image, vl_index o) { float *level ; double sigma, imageSigma ; vl_index op ; assert(self) ; assert(image) ; assert(o >= self->geom.firstOctave) ; assert(o <= self->geom.lastOctave) ; /* * Copy the image to self->geom.octaveFirstSubdivision of octave o, upscaling or * downscaling as needed. */ level = vl_scalespace_get_level(self, VL_MAX(0, o), self->geom.octaveFirstSubdivision) ; copy_and_downsample(level, image, self->geom.width, self->geom.height, VL_MAX(0, o)) ; for (op = -1 ; op >= o ; --op) { VlScaleSpaceOctaveGeometry ogeom = vl_scalespace_get_octave_geometry(self, op + 1) ; float *succLevel = vl_scalespace_get_level(self, op + 1, self->geom.octaveFirstSubdivision) ; level = vl_scalespace_get_level(self, op, self->geom.octaveFirstSubdivision) ; copy_and_upsample(level, succLevel, ogeom.width, ogeom.height) ; } /* * Adjust the smoothing of the first level just initialised, accounting * for the fact that the input image is assumed to be a nominal scale * level. */ sigma = vl_scalespace_get_level_sigma(self, o, self->geom.octaveFirstSubdivision) ; imageSigma = self->geom.nominalScale ; if (sigma > imageSigma) { VlScaleSpaceOctaveGeometry ogeom = vl_scalespace_get_octave_geometry(self, o) ; double deltaSigma = sqrt (sigma*sigma - imageSigma*imageSigma) ; level = vl_scalespace_get_level (self, o, self->geom.octaveFirstSubdivision) ; vl_imsmooth_f (level, ogeom.width, level, ogeom.width, ogeom.height, ogeom.width, deltaSigma / ogeom.step, deltaSigma / ogeom.step) ; } } /** @internal @brief Initialize the first level of an octave from the previous octave ** @param self object. ** @param o octave to initialize. ** ** The function initializes the first level of octave @a o from the ** content of octave o - 1. **/ static void _vl_scalespace_start_octave_from_previous_octave (VlScaleSpace *self, vl_index o) { double sigma, prevSigma ; float *level, *prevLevel ; vl_index prevLevelIndex ; VlScaleSpaceOctaveGeometry ogeom ; assert(self) ; assert(o > self->geom.firstOctave) ; /* must not be the first octave */ assert(o <= self->geom.lastOctave) ; /* * From the previous octave pick the level which is closer to * self->geom.octaveFirstSubdivision in this octave. * The is self->geom.octaveFirstSubdivision + self->numLevels since there are * self->geom.octaveResolution levels in an octave, provided that * this value does not exceed self->geom.octaveLastSubdivision. */ prevLevelIndex = VL_MIN(self->geom.octaveFirstSubdivision + (signed)self->geom.octaveResolution, self->geom.octaveLastSubdivision) ; prevLevel = vl_scalespace_get_level (self, o - 1, prevLevelIndex) ; level = vl_scalespace_get_level (self, o, self->geom.octaveFirstSubdivision) ; ogeom = vl_scalespace_get_octave_geometry(self, o - 1) ; copy_and_downsample (level, prevLevel, ogeom.width, ogeom.height, 1) ; /* * Add remaining smoothing, if any. */ sigma = vl_scalespace_get_level_sigma(self, o, self->geom.octaveFirstSubdivision) ; prevSigma = vl_scalespace_get_level_sigma(self, o - 1, prevLevelIndex) ; if (sigma > prevSigma) { VlScaleSpaceOctaveGeometry ogeom = vl_scalespace_get_octave_geometry(self, o) ; double deltaSigma = sqrt (sigma*sigma - prevSigma*prevSigma) ; level = vl_scalespace_get_level (self, o, self->geom.octaveFirstSubdivision) ; /* todo: this may fail due to an out-of-memory condition */ vl_imsmooth_f (level, ogeom.width, level, ogeom.width, ogeom.height, ogeom.width, deltaSigma / ogeom.step, deltaSigma / ogeom.step) ; } } /** @brief Initialise Scale space with new image ** @param self ::VlScaleSpace object instance. ** @param image image to process. ** ** Compute the data of all the defined octaves and scales of the scale ** space @a self. **/ void vl_scalespace_put_image (VlScaleSpace *self, float const *image) { vl_index o ; _vl_scalespace_start_octave_from_image(self, image, self->geom.firstOctave) ; _vl_scalespace_fill_octave(self, self->geom.firstOctave) ; for (o = self->geom.firstOctave + 1 ; o <= self->geom.lastOctave ; ++o) { _vl_scalespace_start_octave_from_previous_octave(self, o) ; _vl_scalespace_fill_octave(self, o) ; } }