liop.c

Go to the documentation of this file.

/*
Copyright (C) 2013 Hana Sarbortova and Andrea Vedaldi.
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 "liop.h"
#include "mathop.h"
#include "imopv.h"
#include <string.h>

#define DEFAULT_INTENSITY_THRESHOLD -(5.0/255)
#define DEFAULT_RADIUS 6.0
#define DEFAULT_NUM_SPATIAL_BINS 6
#define DEFAULT_NUM_NEIGHBOURS 4

/* ---------------------------------------------------------------- */
/*                                                 Helper functions */
/* ---------------------------------------------------------------- */

static
vl_int factorial(vl_int num)
{
  vl_int result = 1;
  while(num > 1){
    result = num*result;
    num--;
  }
  return result ;
}

VL_INLINE vl_index get_permutation_index(vl_uindex *permutation, vl_size size){
  vl_index index = 0 ;
  vl_index i ;
  vl_index j ;

  for (i = 0 ; i < (signed)size ; ++i) {
    index = index * ((signed)size - i) + permutation[i] ;
    for (j = i + 1 ; j < (signed)size ; ++j) {
      if (permutation[j] > permutation[i]) { permutation[j] -- ; }
    }
  }
  return index ;
}

/* instantiate two quick sort algorithms */
VL_INLINE float patch_cmp (VlLiopDesc * liop, vl_index i, vl_index j)
{
  vl_index ii = liop->patchPermutation[i] ;
  vl_index jj = liop->patchPermutation[j] ;
  return liop->patchIntensities[ii] - liop->patchIntensities[jj] ;
}

VL_INLINE void patch_swap (VlLiopDesc * liop, vl_index i, vl_index j)
{
  vl_index tmp = liop->patchPermutation[i] ;
  liop->patchPermutation[i] = liop->patchPermutation[j] ;
  liop->patchPermutation[j] = tmp ;
}

#define VL_QSORT_prefix patch
#define VL_QSORT_array VlLiopDesc*
#define VL_QSORT_cmp patch_cmp
#define VL_QSORT_swap patch_swap
#include "qsort-def.h"

VL_INLINE float neigh_cmp (VlLiopDesc * liop, vl_index i, vl_index j)
{
  vl_index ii = liop->neighPermutation[i] ;
  vl_index jj = liop->neighPermutation[j] ;
  return liop->neighIntensities[ii] - liop->neighIntensities[jj] ;
}

VL_INLINE void neigh_swap (VlLiopDesc * liop, vl_index i, vl_index j)
{
  vl_index tmp = liop->neighPermutation[i] ;
  liop->neighPermutation[i] = liop->neighPermutation[j] ;
  liop->neighPermutation[j] = tmp ;
}

#define VL_QSORT_prefix neigh
#define VL_QSORT_array VlLiopDesc*
#define VL_QSORT_cmp neigh_cmp
#define VL_QSORT_swap neigh_swap
#include "qsort-def.h"

/* ---------------------------------------------------------------- */
/*                                            Construct and destroy */
/* ---------------------------------------------------------------- */

VlLiopDesc *
vl_liopdesc_new (vl_int numNeighbours, vl_int numSpatialBins,
                 float radius, vl_size sideLength)
{
  vl_index i, t ;
  VlLiopDesc * self = vl_calloc(sizeof(VlLiopDesc), 1);

  assert(radius <= sideLength/2) ;

  self->numNeighbours = numNeighbours ;
  self->numSpatialBins = numSpatialBins ;
  self->neighRadius = radius ;
  self->intensityThreshold = DEFAULT_INTENSITY_THRESHOLD ;

  self->dimension = factorial(numNeighbours) * numSpatialBins ;

  /*
   Precompute a list of pixels within a circular patch inside
   the square image. Leave a suitable marging for sampling around
   these pixels.
   */

  self->patchSize = 0 ;
  self->patchPixels = vl_malloc(sizeof(vl_uindex)*sideLength*sideLength) ;
  self->patchSideLength = sideLength ;

  {
    vl_index x, y ;
    vl_index center = (sideLength - 1) / 2 ;
    double t = center - radius + 0.6 ;
    vl_index t2 = (vl_index) (t * t) ;
    for (y = 0 ; y < (signed)sideLength ; ++y) {
      for (x = 0 ; x < (signed)sideLength ; ++x) {
        vl_index dx = x - center ;
        vl_index dy = y - center ;
        if (x == 0 && y == 0) continue ;
        if (dx*dx + dy*dy <= t2) {
          self->patchPixels[self->patchSize++] = x + y * sideLength ;
        }
      }
    }
  }

  self->patchIntensities = vl_malloc(sizeof(vl_uindex)*self->patchSize) ;
  self->patchPermutation = vl_malloc(sizeof(vl_uindex)*self->patchSize) ;

  /*
   Precompute the samples in the circular neighbourhood of each
   measurement point.
   */

  self->neighPermutation = vl_malloc(sizeof(vl_uindex) * self->numNeighbours) ;
  self->neighIntensities = vl_malloc(sizeof(float) * self->numNeighbours) ;
  self->neighSamplesX = vl_calloc(sizeof(double), self->numNeighbours * self->patchSize) ;
  self->neighSamplesY = vl_calloc(sizeof(double), self->numNeighbours * self->patchSize) ;

  for (i = 0 ; i < (signed)self->patchSize ; ++i) {
    vl_index pixel ;
    double x, y ;
    double dangle = 2*VL_PI / (double)self->numNeighbours ;
    double angle0 ;
    vl_index center = (sideLength - 1) / 2 ;

    pixel = self->patchPixels[i] ;
    x = (pixel % (signed)self->patchSideLength) - center ;
    y = (pixel / (signed)self->patchSideLength) - center ;

    angle0 = atan2(y,x) ;

    for (t = 0 ; t < (signed)self->numNeighbours ; ++t) {
      double x1 = x + radius * cos(angle0 + dangle * t) + center ;
      double y1 = y + radius * sin(angle0 + dangle * t) + center ;
      self->neighSamplesX[t + (signed)self->numNeighbours * i] = x1 ;
      self->neighSamplesY[t + (signed)self->numNeighbours * i] = y1 ;
    }
  }
  return self ;
}

VlLiopDesc * vl_liopdesc_new_basic (vl_size sideLength)
{
  return vl_liopdesc_new(DEFAULT_NUM_NEIGHBOURS,
                         DEFAULT_NUM_SPATIAL_BINS,
                         DEFAULT_RADIUS,
                         sideLength) ;
}

void
vl_liopdesc_delete (VlLiopDesc * self)
{
  vl_free (self->patchPixels) ;
  vl_free (self->patchIntensities) ;
  vl_free (self->patchPermutation) ;
  vl_free (self->neighPermutation) ;
  vl_free (self->neighIntensities) ;
  vl_free (self->neighSamplesX) ;
  vl_free (self->neighSamplesY) ;
  vl_free (self) ;
}

/* ---------------------------------------------------------------- */
/*                                          Compute LIOP descriptor */
/* ---------------------------------------------------------------- */

void
vl_liopdesc_process (VlLiopDesc * self, float * desc, float const * patch)
{
  vl_index i,t ;
  vl_index offset,numPermutations ;
  vl_index spatialBinArea, spatialBinEnd, spatialBinIndex ;
  float threshold ;

  memset(desc, 0, sizeof(float) * self->dimension) ;

  /*
   * Sort pixels in the patch by increasing intensity.
   */

  for (i = 0 ; i < (signed)self->patchSize ; ++i) {
    vl_index pixel = self->patchPixels[i] ;
    self->patchIntensities[i] = patch[pixel] ;
    self->patchPermutation[i] = i ;
  }
  patch_sort(self, self->patchSize) ;

  /*
   * Tune the threshold if needed.
   */

  if (self->intensityThreshold < 0) {
    i = self->patchPermutation[0] ;
    t = self->patchPermutation[self->patchSize-1] ;
    threshold = - self->intensityThreshold
    * (self->patchIntensities[t] - self->patchIntensities[i]);
  } else {
    threshold = self->intensityThreshold ;
  }

  /*
   * Process pixels in order of increasing intenisity, dividing them into
   * spatial bins on the fly.
   */

  numPermutations = factorial(self->numNeighbours) ;
  spatialBinArea = self->patchSize / self->numSpatialBins ;
  spatialBinEnd = spatialBinArea ;
  spatialBinIndex = 0 ;
  offset = 0 ;

  for (i = 0 ; i < (signed)self->patchSize ; ++i) {
    vl_index permIndex ;
    double *sx, *sy ;

    /* advance to the next spatial bin if needed */
    if (i >= (signed)spatialBinEnd && spatialBinIndex < (signed)self->numSpatialBins - 1) {
      spatialBinEnd += spatialBinArea ;
      spatialBinIndex ++ ;
      offset += numPermutations ;
    }

    /* get intensities of neighbours of the current patch element and sort them */
    sx = self->neighSamplesX + self->numNeighbours * self->patchPermutation[i] ;
    sy = self->neighSamplesY + self->numNeighbours * self->patchPermutation[i] ;
    for (t = 0 ; t < self->numNeighbours ; ++t) {
      double x = *sx++ ;
      double y = *sy++ ;

      /* bilinear interpolation */
      vl_index ix = vl_floor_d(x) ;
      vl_index iy = vl_floor_d(y) ;

      double wx = x - ix ;
      double wy = y - iy ;

      double a = 0, b = 0, c = 0, d = 0 ;

      int L = (int) self->patchSideLength ;

      if (ix >= 0   && iy >= 0  ) { a = patch[ix   + iy * L] ; }
      if (ix <  L-1 && iy >= 0  ) { b = patch[ix+1 + iy * L] ; }
      if (ix >= 0   && iy <  L-1) { c = patch[ix   + (iy+1) * L] ; }
      if (ix <  L-1 && iy <  L-1) { d = patch[ix+1 + (iy+1) * L] ; }

      self->neighPermutation[t] = t;
      self->neighIntensities[t] = (1 - wy) * (a + (b - a) * wx) + wy * (c + (d - c) * wx) ;
    }
    neigh_sort (self, self->numNeighbours) ;

    /* get permutation index */
    permIndex = get_permutation_index(self->neighPermutation, self->numNeighbours);

    /*
     * Compute weight according to difference in intensity values and
     * accumulate.
     */
    {
      int k, t ;
      float weight = 0 ;
      for(k = 0; k < self->numNeighbours ; ++k) {
        for(t = k + 1; t < self->numNeighbours; ++t){
          double a = self->neighIntensities[k] ;
          double b = self->neighIntensities[t] ;
          weight += (a > b + threshold || b > a + threshold) ;
        }
      }
      desc[permIndex + offset] += weight ;
    }
  }

  /* normalization */
  {
    float norm = 0;
    for(i = 0; i < (signed)self->dimension; i++) {
      norm += desc[i]*desc[i];
    }
    norm = VL_MAX(sqrt(norm), 1e-12) ;
    for(i = 0; i < (signed)self->dimension; i++){
      desc[i] /= norm ;
    }
  }
}


/* ---------------------------------------------------------------- */
/*                                              Getters and setters */
/* ---------------------------------------------------------------- */

vl_size
vl_liopdesc_get_dimension (VlLiopDesc const * self)
{
  return self->dimension ;
}


vl_size
vl_liopdesc_get_num_neighbours (VlLiopDesc const * self)
{
  assert(self) ;
  return self->numNeighbours ;
}

float
vl_liopdesc_get_intensity_threshold (VlLiopDesc const * self)
{
  assert(self) ;
  return self->intensityThreshold ;
}

void
vl_liopdesc_set_intensity_threshold (VlLiopDesc * self, float x)
{
  assert(self) ;
  self->intensityThreshold = x ;
}

double
vl_liopdesc_get_neighbourhood_radius (VlLiopDesc const * self)
{
  assert(self) ;
  return self->neighRadius ;
}

vl_size
vl_liopdesc_get_num_spatial_bins (VlLiopDesc const * self)
{
  assert(self) ;
  return self->numSpatialBins ;
}