homkermap.c

Go to the documentation of this file.

/*
Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
Copyright (C) 2013 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).
*/

/* ---------------------------------------------------------------- */
#ifndef VL_HOMKERMAP_INSTANTIATING
/* ---------------------------------------------------------------- */

#include "homkermap.h"
#include "mathop.h"
#include <math.h>

struct _VlHomogeneousKernelMap
{
  VlHomogeneousKernelType kernelType ;
  double gamma ;
  VlHomogeneousKernelMapWindowType windowType ;
  vl_size order ;
  double period ;
  vl_size numSubdivisions ;
  double subdivision  ;
  vl_index minExponent ;
  vl_index maxExponent ;
  double * table ;
} ;

VL_INLINE double
vl_homogeneouskernelmap_get_spectrum (VlHomogeneousKernelMap const * self, double omega)
{
  assert (self) ;
  switch (self->kernelType) {
    case VlHomogeneousKernelIntersection:
      return (2.0 / VL_PI) / (1 + 4 * omega*omega) ;
    case VlHomogeneousKernelChi2:
      return 2.0 / (exp(VL_PI * omega) + exp(-VL_PI * omega)) ;
    case VlHomogeneousKernelJS:
      return (2.0 / log(4.0)) *
      2.0 / (exp(VL_PI * omega) + exp(-VL_PI * omega)) /
      (1 + 4 * omega*omega) ;
    default:
      abort() ;
  }
}

/* helper */
VL_INLINE double sinc(double x)
{
  if (x == 0.0) return 1.0 ;
  return sin(x) / x ;
}

VL_INLINE double
vl_homogeneouskernelmap_get_smooth_spectrum (VlHomogeneousKernelMap const * self, double omega)
{
  double kappa_hat = 0 ;
  double omegap ;
  double epsilon = 1e-2 ;
  double const omegaRange = 2.0 / (self->period * epsilon) ;
  double const domega = 2 * omegaRange / (2 * 1024.0 + 1) ;
  assert (self) ;
  switch (self->windowType) {
    case VlHomogeneousKernelMapWindowUniform:
      kappa_hat = vl_homogeneouskernelmap_get_spectrum(self, omega) ;
      break ;
    case VlHomogeneousKernelMapWindowRectangular:
      for (omegap = - omegaRange ; omegap <= omegaRange ; omegap += domega) {
        double win = sinc((self->period/2.0) * omegap) ;
        win *= (self->period/(2.0*VL_PI)) ;
        kappa_hat += win * vl_homogeneouskernelmap_get_spectrum(self, omegap + omega) ;
      }
      kappa_hat *= domega ;
      /* project on the postivie orthant (see PAMI) */
      kappa_hat = VL_MAX(kappa_hat, 0.0) ;
      break ;
    default:
      abort() ;
  }
  return kappa_hat ;
}

/* ---------------------------------------------------------------- */
/*                                     Constructors and destructors */
/* ---------------------------------------------------------------- */

VlHomogeneousKernelMap *
vl_homogeneouskernelmap_new (VlHomogeneousKernelType kernelType,
                             double gamma,
                             vl_size order,
                             double period,
                             VlHomogeneousKernelMapWindowType windowType)
{
  int tableWidth, tableHeight ;
  VlHomogeneousKernelMap * self = vl_malloc(sizeof(VlHomogeneousKernelMap)) ;
  if (! self) return NULL ;

  assert(gamma > 0) ;

  assert(kernelType == VlHomogeneousKernelIntersection ||
         kernelType == VlHomogeneousKernelChi2 ||
         kernelType == VlHomogeneousKernelJS) ;

  assert(windowType == VlHomogeneousKernelMapWindowUniform ||
         windowType == VlHomogeneousKernelMapWindowRectangular) ;

  if (period < 0) {
    switch (windowType) {
    case VlHomogeneousKernelMapWindowUniform:
      switch (kernelType) {
      case VlHomogeneousKernelChi2:         period = 5.86 * sqrt(order + 0)  + 3.65 ; break ;
      case VlHomogeneousKernelJS:           period = 6.64 * sqrt(order + 0)  + 7.24 ; break ;
      case VlHomogeneousKernelIntersection: period = 2.38 * log(order + 0.8) + 5.6 ; break ;
      }
      break ;
    case VlHomogeneousKernelMapWindowRectangular:
      switch (kernelType) {
      case VlHomogeneousKernelChi2:         period = 8.80 * sqrt(order + 4.44) - 12.6 ; break ;
      case VlHomogeneousKernelJS:           period = 9.63 * sqrt(order + 1.00) - 2.93;  break ;
      case VlHomogeneousKernelIntersection: period = 2.00 * log(order + 0.99)  + 3.52 ; break ;
      }
      break ;
    }
    period = VL_MAX(period, 1.0) ;
  }

  self->kernelType = kernelType ;
  self->windowType = windowType ;
  self->gamma = gamma ;
  self->order = order ;
  self->period = period ;
  self->numSubdivisions = 8 + 8*order ;
  self->subdivision = 1.0 / self->numSubdivisions ;
  self->minExponent = -20 ;
  self->maxExponent = 8 ;

  tableHeight = (int) (2*self->order + 1) ;
  tableWidth = (int) (self->numSubdivisions * (self->maxExponent - self->minExponent + 1)) ;
  self->table = vl_malloc (sizeof(double) *
                           (tableHeight * tableWidth + 2*(1+self->order))) ;
  if (! self->table) {
    vl_free(self) ;
    return NULL ;
  }

  {
    vl_index exponent ;
    vl_uindex i, j ;
    double * tablep = self->table ;
    double * kappa = self->table + tableHeight * tableWidth ;
    double * freq = kappa + (1+self->order) ;
    double L = 2.0 * VL_PI / self->period ;

    /* precompute the sampled periodicized spectrum */
    j = 0 ;
    i = 0 ;
    while (i <= self->order) {
      freq[i] = j ;
      kappa[i] = vl_homogeneouskernelmap_get_smooth_spectrum(self, j * L) ;
      ++ j ;
      if (kappa[i] > 0 || j >= 3*i) ++ i ;
    }

    /* fill table */
    for (exponent  = self->minExponent ;
         exponent <= self->maxExponent ; ++ exponent) {

      double x, Lxgamma, Llogx, xgamma ;
      double sqrt2kappaLxgamma ;
      double mantissa = 1.0 ;

      for (i = 0 ; i < self->numSubdivisions ;
           ++i, mantissa += self->subdivision) {
        x = ldexp(mantissa, (int)exponent) ;
        xgamma = pow(x, self->gamma) ;
        Lxgamma = L * xgamma ;
        Llogx = L * log(x) ;

        *tablep++ = sqrt(Lxgamma * kappa[0]) ;
        for (j = 1 ; j <= self->order ; ++j) {
          sqrt2kappaLxgamma = sqrt(2.0 * Lxgamma * kappa[j]) ;
          *tablep++ = sqrt2kappaLxgamma * cos(freq[j] * Llogx) ;
          *tablep++ = sqrt2kappaLxgamma * sin(freq[j] * Llogx) ;
        }
      } /* next mantissa */
    } /* next exponent */
  }
  return self ;
}

void
vl_homogeneouskernelmap_delete (VlHomogeneousKernelMap * self)
{
  vl_free(self->table) ;
  self->table = NULL ;
  vl_free(self) ;
}

/* ---------------------------------------------------------------- */
/*                                     Retrieve data and parameters */
/* ---------------------------------------------------------------- */

vl_size
vl_homogeneouskernelmap_get_order (VlHomogeneousKernelMap const * self)
{
  assert(self) ;
  return self->order ;
}

vl_size
vl_homogeneouskernelmap_get_dimension (VlHomogeneousKernelMap const * self)
{
  assert(self) ;
  return 2 * self->order + 1 ;
}

VlHomogeneousKernelType
vl_homogeneouskernelmap_get_kernel_type (VlHomogeneousKernelMap const * self)
{
  assert(self) ;
  return self->kernelType ;
}

VlHomogeneousKernelMapWindowType
vl_homogeneouskernelmap_get_window_type (VlHomogeneousKernelMap const * self)
{
  assert(self) ;
  return self->windowType ;
}

/* ---------------------------------------------------------------- */
/*                                                     Process data */
/* ---------------------------------------------------------------- */

#define FLT VL_TYPE_FLOAT
#define VL_HOMKERMAP_INSTANTIATING
#include "homkermap.c"

#define FLT VL_TYPE_DOUBLE
#define VL_HOMKERMAP_INSTANTIATING
#include "homkermap.c"

/* VL_HOMKERMAP_INSTANTIATING */
#endif

/* ---------------------------------------------------------------- */
#ifdef VL_HOMKERMAP_INSTANTIATING
/* ---------------------------------------------------------------- */

#include "float.th"

void
VL_XCAT(vl_homogeneouskernelmap_evaluate_,SFX)
(VlHomogeneousKernelMap const * self,
 T * destination,
 vl_size stride,
 double x)
{
  /* break value into exponent and mantissa */
  int exponent ;
  int unsigned j ;
  double mantissa = frexp(x, &exponent) ;
  double sign = (mantissa >= 0.0) ? +1.0 : -1.0 ;
  mantissa *= 2*sign ;
  exponent -- ;

  if (mantissa == 0 ||
      exponent <= self->minExponent ||
      exponent >= self->maxExponent) {
    for (j = 0 ; j < 2*self->order+1 ; ++j) {
      *destination = (T) 0.0 ;
      destination += stride ;
    }
    return  ;
  }
  {
    vl_size featureDimension = 2*self->order + 1 ;
    double const * v1 = self->table +
    (exponent - self->minExponent) * self->numSubdivisions * featureDimension ;
    double const * v2 ;
    double f1, f2 ;

    mantissa -= 1.0 ;
    while (mantissa >= self->subdivision) {
      mantissa -= self->subdivision ;
      v1 += featureDimension ;
    }
    v2 = v1 + featureDimension ;
    for (j = 0 ; j < featureDimension ; ++j) {
      f1 = *v1++ ;
      f2 = *v2++ ;
      *destination = (T) sign * ((f2 - f1) * (self->numSubdivisions * mantissa) + f1) ;
      destination += stride ;
    }
  }
}

#undef FLT
#undef VL_HOMKERMAP_INSTANTIATING
/* VL_HOMKERMAP_INSTANTIATING */
#endif