scalespace.c

Go to the documentation of this file.

/*
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).
*/

#include "scalespace.h"
#include "mathop.h"

#include <assert.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdio.h>

struct _VlScaleSpace
{
  VlScaleSpaceGeometry geom ; 
  float **octaves ; 
} ;

/* ---------------------------------------------------------------- */
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)

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 ;
}

VlScaleSpaceGeometry
vl_scalespace_get_geometry (VlScaleSpace const * self)
{
  return self->geom ;
}

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 ;
}

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) ;
}

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) ;
}

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) ;
}

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 ;
  }
}

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 ;
      }
    }
  }
}

/* ---------------------------------------------------------------- */
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) ;
}

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 ;
}

/* ---------------------------------------------------------------- */
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 ;
}

/* ---------------------------------------------------------------- */
VlScaleSpace *
vl_scalespace_new_shallow_copy (VlScaleSpace* self)
{
  return vl_scalespace_new_with_geometry (self->geom) ;
}

/* ---------------------------------------------------------------- */
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) ;
  }
}

/* ---------------------------------------------------------------- */

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) ;
  }
}

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) ;
  }
}

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) ;
  }
}

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) ;
  }
}