dsift.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 "dsift.h"
#include "pgm.h"
#include "mathop.h"
#include "imopv.h"
#include <math.h>
#include <string.h>

float *
_vl_dsift_new_kernel (int binSize, int numBins, int binIndex, double windowSize)
{
  int filtLen = 2 * binSize - 1 ;
  float * ker = vl_malloc (sizeof(float) * filtLen) ;
  float * kerIter = ker ;
  float delta = binSize * (binIndex - 0.5F * (numBins - 1)) ;
  /*
  float sigma = 0.5F * ((numBins - 1) * binSize + 1) ;
  float sigma = 0.5F * ((numBins) * binSize) ;
  */
  float sigma = (float) binSize * (float) windowSize ;
  int x ;

  for (x = - binSize + 1 ; x <= + binSize - 1 ; ++ x) {
    float z = (x - delta) / sigma ;
    *kerIter++ = (1.0F - fabsf(x) / binSize) *
      ((binIndex >= 0) ? expf(- 0.5F * z*z) : 1.0F) ;
  }
  return ker ;
}

static float
_vl_dsift_get_bin_window_mean
(int binSize, int numBins, int binIndex, double windowSize)
{
  float delta = binSize * (binIndex - 0.5F * (numBins - 1)) ;
  /*float sigma = 0.5F * ((numBins - 1) * binSize + 1) ;*/
  float sigma = (float) binSize * (float) windowSize ;
  int x ;

  float acc = 0.0 ;
  for (x = - binSize + 1 ; x <= + binSize - 1 ; ++ x) {
    float z = (x - delta) / sigma ;
    acc += ((binIndex >= 0) ? expf(- 0.5F * z*z) : 1.0F) ;
  }
  return acc /= (2 * binSize - 1) ;
}

VL_INLINE float
_vl_dsift_normalize_histogram (float * begin, float * end)
{
  float * iter ;
  float  norm = 0.0F ;

  for (iter = begin ; iter < end ; ++ iter) {
    norm += (*iter) * (*iter) ;
  }
  norm = vl_fast_sqrt_f (norm) + VL_EPSILON_F ;

  for (iter = begin; iter < end ; ++ iter) {
    *iter /= norm ;
  }
  return norm ;
}

static void
_vl_dsift_free_buffers (VlDsiftFilter* self)
{
  if (self->frames) {
    vl_free(self->frames) ;
    self->frames = NULL ;
  }
  if (self->descrs) {
    vl_free(self->descrs) ;
    self->descrs = NULL ;
  }
  if (self->grads) {
    int t ;
    for (t = 0 ; t < self->numGradAlloc ; ++t)
      if (self->grads[t]) vl_free(self->grads[t]) ;
    vl_free(self->grads) ;
    self->grads = NULL ;
  }
  self->numFrameAlloc = 0 ;
  self->numBinAlloc = 0 ;
  self->numGradAlloc = 0 ;
}

VL_EXPORT void
_vl_dsift_update_buffers (VlDsiftFilter * self)
{
  int x1 = self->boundMinX ;
  int x2 = self->boundMaxX ;
  int y1 = self->boundMinY ;
  int y2 = self->boundMaxY ;

  int rangeX = x2 - x1 - (self->geom.numBinX - 1) * self->geom.binSizeX ;
  int rangeY = y2 - y1 - (self->geom.numBinY - 1) * self->geom.binSizeY ;

  int numFramesX = (rangeX >= 0) ? rangeX / self->stepX + 1 : 0 ;
  int numFramesY = (rangeY >= 0) ? rangeY / self->stepY + 1 : 0 ;

  self->numFrames = numFramesX * numFramesY ;
  self->descrSize = self->geom.numBinT *
                    self->geom.numBinX *
                    self->geom.numBinY ;
}

static void
_vl_dsift_alloc_buffers (VlDsiftFilter* self)
{
  _vl_dsift_update_buffers (self) ;
  {
    int numFrameAlloc = vl_dsift_get_keypoint_num (self) ;
    int numBinAlloc   = vl_dsift_get_descriptor_size (self) ;
    int numGradAlloc  = self->geom.numBinT ;

    /* see if we need to update the buffers */
    if (numBinAlloc != self->numBinAlloc ||
        numGradAlloc != self->numGradAlloc ||
        numFrameAlloc != self->numFrameAlloc) {

      int t ;

      _vl_dsift_free_buffers(self) ;

      self->frames = vl_malloc(sizeof(VlDsiftKeypoint) * numFrameAlloc) ;
      self->descrs = vl_malloc(sizeof(float) * numBinAlloc * numFrameAlloc) ;
      self->grads  = vl_malloc(sizeof(float*) * numGradAlloc) ;
      for (t = 0 ; t < numGradAlloc ; ++t) {
        self->grads[t] =
          vl_malloc(sizeof(float) * self->imWidth * self->imHeight) ;
      }
      self->numBinAlloc = numBinAlloc ;
      self->numGradAlloc = numGradAlloc ;
      self->numFrameAlloc = numFrameAlloc ;
    }
  }
}

VL_EXPORT VlDsiftFilter *
vl_dsift_new (int imWidth, int imHeight)
{
  VlDsiftFilter * self = vl_malloc (sizeof(VlDsiftFilter)) ;
  self->imWidth  = imWidth ;
  self->imHeight = imHeight ;

  self->stepX = 5 ;
  self->stepY = 5 ;

  self->boundMinX = 0 ;
  self->boundMinY = 0 ;
  self->boundMaxX = imWidth - 1 ;
  self->boundMaxY = imHeight - 1 ;

  self->geom.numBinX = 4 ;
  self->geom.numBinY = 4 ;
  self->geom.numBinT = 8 ;
  self->geom.binSizeX = 5 ;
  self->geom.binSizeY = 5 ;

  self->useFlatWindow = VL_FALSE ;
  self->windowSize = 2.0 ;

  self->convTmp1 = vl_malloc(sizeof(float) * self->imWidth * self->imHeight) ;
  self->convTmp2 = vl_malloc(sizeof(float) * self->imWidth * self->imHeight) ;

  self->numBinAlloc = 0 ;
  self->numFrameAlloc = 0 ;
  self->numGradAlloc = 0 ;

  self->descrSize = 0 ;
  self->numFrames = 0 ;
  self->grads = NULL ;
  self->frames = NULL ;
  self->descrs = NULL ;

  _vl_dsift_update_buffers(self) ;
  return self ;
}

VL_EXPORT VlDsiftFilter *
vl_dsift_new_basic (int imWidth, int imHeight, int step, int binSize)
{
  VlDsiftFilter* self = vl_dsift_new(imWidth, imHeight) ;
  VlDsiftDescriptorGeometry geom = *vl_dsift_get_geometry(self) ;
  geom.binSizeX = binSize ;
  geom.binSizeY = binSize ;
  vl_dsift_set_geometry(self, &geom) ;
  vl_dsift_set_steps(self, step, step) ;
  return self ;
}

VL_EXPORT void
vl_dsift_delete (VlDsiftFilter * self)
{
  _vl_dsift_free_buffers (self) ;
  if (self->convTmp2) vl_free (self->convTmp2) ;
  if (self->convTmp1) vl_free (self->convTmp1) ;
  vl_free (self) ;
}


VL_INLINE void
_vl_dsift_with_gaussian_window (VlDsiftFilter * self)
{
  int binx, biny, bint ;
  int framex, framey ;
  float *xker, *yker ;

  int Wx = self->geom.binSizeX - 1 ;
  int Wy = self->geom.binSizeY - 1 ;

  for (biny = 0 ; biny < self->geom.numBinY ; ++biny) {

    yker = _vl_dsift_new_kernel (self->geom.binSizeY,
                                 self->geom.numBinY,
                                 biny,
                                 self->windowSize) ;

    for (binx = 0 ; binx < self->geom.numBinX ; ++binx) {

      xker = _vl_dsift_new_kernel(self->geom.binSizeX,
                                  self->geom.numBinX,
                                  binx,
                                  self->windowSize) ;

      for (bint = 0 ; bint < self->geom.numBinT ; ++bint) {

        vl_imconvcol_vf (self->convTmp1, self->imHeight,
                         self->grads[bint], self->imWidth, self->imHeight,
                         self->imWidth,
                         yker, -Wy, +Wy, 1,
                         VL_PAD_BY_CONTINUITY|VL_TRANSPOSE) ;

        vl_imconvcol_vf (self->convTmp2, self->imWidth,
                         self->convTmp1, self->imHeight, self->imWidth,
                         self->imHeight,
                         xker, -Wx, +Wx, 1,
                         VL_PAD_BY_CONTINUITY|VL_TRANSPOSE) ;

        {
          float *dst = self->descrs
            + bint
            + binx * self->geom.numBinT
            + biny * (self->geom.numBinX * self->geom.numBinT)  ;

          float *src = self->convTmp2 ;

          int frameSizeX = self->geom.binSizeX * (self->geom.numBinX - 1) + 1 ;
          int frameSizeY = self->geom.binSizeY * (self->geom.numBinY - 1) + 1 ;
          int descrSize = vl_dsift_get_descriptor_size (self) ;

          for (framey  = self->boundMinY ;
               framey <= self->boundMaxY - frameSizeY + 1 ;
               framey += self->stepY) {
            for (framex  = self->boundMinX ;
                 framex <= self->boundMaxX - frameSizeX + 1 ;
                 framex += self->stepX) {
              *dst = src [(framex + binx * self->geom.binSizeX) * 1 +
                          (framey + biny * self->geom.binSizeY) * self->imWidth]  ;
              dst += descrSize ;
            } /* framex */
          } /* framey */
        }

      } /* for bint */
      vl_free (xker) ;
    } /* for binx */
    vl_free (yker) ;
  } /* for biny */
}

VL_INLINE void
_vl_dsift_with_flat_window (VlDsiftFilter* self)
{
  int binx, biny, bint ;
  int framex, framey ;

  /* for each orientation bin */
  for (bint = 0 ; bint < self->geom.numBinT ; ++bint) {

    vl_imconvcoltri_f (self->convTmp1, self->imHeight,
                       self->grads [bint], self->imWidth, self->imHeight,
                       self->imWidth,
                       self->geom.binSizeY, /* filt size */
                       1, /* subsampling step */
                       VL_PAD_BY_CONTINUITY|VL_TRANSPOSE) ;

    vl_imconvcoltri_f (self->convTmp2, self->imWidth,
                       self->convTmp1, self->imHeight, self->imWidth,
                       self->imHeight,
                       self->geom.binSizeX,
                       1,
                       VL_PAD_BY_CONTINUITY|VL_TRANSPOSE) ;

    for (biny = 0 ; biny < self->geom.numBinY ; ++biny) {

      /*
      This fast version of DSIFT does not use a proper Gaussian
      weighting scheme for the gradiens that are accumulated on the
      spatial bins. Instead each spatial bins is accumulated based on
      the triangular kernel only, equivalent to bilinear interpolation
      plus a flat, rather than Gaussian, window. Eventually, however,
      the magnitude of the spatial bins in the SIFT descriptor is
      reweighted by the average of the Gaussian window on each bin.
      */

      float wy = _vl_dsift_get_bin_window_mean
        (self->geom.binSizeY, self->geom.numBinY, biny,
         self->windowSize) ;

      /* The convolution functions vl_imconvcoltri_* convolve by a
       * triangular kernel with unit integral. Instead for SIFT the
       * triangular kernel should have unit height. This is
       * compensated for by multiplying by the bin size:
       */

      wy *= self->geom.binSizeY ;

      for (binx = 0 ; binx < self->geom.numBinX ; ++binx) {
        float w ;
        float wx = _vl_dsift_get_bin_window_mean (self->geom.binSizeX,
                                                  self->geom.numBinX,
                                                  binx,
                                                  self->windowSize) ;

        float *dst = self->descrs
          + bint
          + binx * self->geom.numBinT
          + biny * (self->geom.numBinX * self->geom.numBinT)  ;

        float *src = self->convTmp2 ;

        int frameSizeX = self->geom.binSizeX * (self->geom.numBinX - 1) + 1 ;
        int frameSizeY = self->geom.binSizeY * (self->geom.numBinY - 1) + 1 ;
        int descrSize = vl_dsift_get_descriptor_size (self) ;

        wx *= self->geom.binSizeX ;
        w = wx * wy ;

        for (framey  = self->boundMinY ;
             framey <= self->boundMaxY - frameSizeY + 1 ;
             framey += self->stepY) {
          for (framex  = self->boundMinX ;
               framex <= self->boundMaxX - frameSizeX + 1 ;
               framex += self->stepX) {
            *dst = w * src [(framex + binx * self->geom.binSizeX) * 1 +
                            (framey + biny * self->geom.binSizeY) * self->imWidth]  ;
            dst += descrSize ;
          } /* framex */
        } /* framey */
      } /* binx */
    } /* biny */
  } /* bint */
}

void vl_dsift_process (VlDsiftFilter* self, float const* im)
{
  int t, x, y ;

  /* update buffers */
  _vl_dsift_alloc_buffers (self) ;

  /* clear integral images */
  for (t = 0 ; t < self->geom.numBinT ; ++t)
    memset (self->grads[t], 0,
            sizeof(float) * self->imWidth * self->imHeight) ;

#undef at
#define at(x,y) (im[(y)*self->imWidth+(x)])

  /* Compute gradients, their norm, and their angle */

  for (y = 0 ; y < self->imHeight ; ++ y) {
    for (x = 0 ; x < self->imWidth ; ++ x) {
      float gx, gy ;
      float angle, mod, nt, rbint ;
      int bint ;

      /* y derivative */
      if (y == 0) {
        gy = at(x,y+1) - at(x,y) ;
      } else if (y == self->imHeight - 1) {
        gy = at(x,y) - at(x,y-1) ;
      } else {
        gy = 0.5F * (at(x,y+1) - at(x,y-1)) ;
      }

      /* x derivative */
      if (x == 0) {
        gx = at(x+1,y) - at(x,y) ;
      } else if (x == self->imWidth - 1) {
        gx = at(x,y) - at(x-1,y) ;
      } else {
        gx = 0.5F * (at(x+1,y) - at(x-1,y)) ;
      }

      /* angle and modulus */
      angle = vl_fast_atan2_f (gy,gx) ;
      mod = vl_fast_sqrt_f (gx*gx + gy*gy) ;

      /* quantize angle */
      nt = vl_mod_2pi_f (angle) * (self->geom.numBinT / (2*VL_PI)) ;
      bint = (int) vl_floor_f (nt) ;
      rbint = nt - bint ;

      /* write it back */
      self->grads [(bint    ) % self->geom.numBinT][x + y * self->imWidth] = (1 - rbint) * mod ;
      self->grads [(bint + 1) % self->geom.numBinT][x + y * self->imWidth] = (    rbint) * mod ;
    }
  }

  if (self->useFlatWindow) {
    _vl_dsift_with_flat_window(self) ;
  } else {
    _vl_dsift_with_gaussian_window(self) ;
  }

  {
    VlDsiftKeypoint* frameIter = self->frames ;
    float * descrIter = self->descrs ;
    int framex, framey, bint ;

    int frameSizeX = self->geom.binSizeX * (self->geom.numBinX - 1) + 1 ;
    int frameSizeY = self->geom.binSizeY * (self->geom.numBinY - 1) + 1 ;
    int descrSize = vl_dsift_get_descriptor_size (self) ;

    float deltaCenterX = 0.5F * self->geom.binSizeX * (self->geom.numBinX - 1) ;
    float deltaCenterY = 0.5F * self->geom.binSizeY * (self->geom.numBinY - 1) ;

    float normConstant = frameSizeX * frameSizeY ;

    for (framey  = self->boundMinY ;
         framey <= self->boundMaxY - frameSizeY + 1 ;
         framey += self->stepY) {

      for (framex  = self->boundMinX ;
           framex <= self->boundMaxX - frameSizeX + 1 ;
           framex += self->stepX) {

        frameIter->x    = framex + deltaCenterX ;
        frameIter->y    = framey + deltaCenterY ;

        /* mass */
        {
          float mass = 0 ;
          for (bint = 0 ; bint < descrSize ; ++ bint)
            mass += descrIter[bint] ;
          mass /= normConstant ;
          frameIter->norm = mass ;
        }

        /* L2 normalize */
        _vl_dsift_normalize_histogram (descrIter, descrIter + descrSize) ;

        /* clamp */
        for(bint = 0 ; bint < descrSize ; ++ bint)
          if (descrIter[bint] > 0.2F) descrIter[bint] = 0.2F ;

        /* L2 normalize */
        _vl_dsift_normalize_histogram (descrIter, descrIter + descrSize) ;

        frameIter ++ ;
        descrIter += descrSize ;
      } /* for framex */
    } /* for framey */
  }
}