covdet.c

Go to the documentation of this file.

/*
Copyright (C) 2013-14 Andrea Vedaldi.
Copyright (C) 2012 Karel Lenc, Andrea Vedaldi and Michal Perdoch.
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 "covdet.h"
#include <string.h>

static int
_vl_resize_buffer (void ** buffer, vl_size * bufferSize, vl_size targetSize) {
  void * newBuffer ;
  if (*buffer == NULL) {
    *buffer = vl_malloc(targetSize) ;
    if (*buffer) {
      *bufferSize = targetSize ;
      return VL_ERR_OK ;
    } else {
      *bufferSize = 0 ;
      return VL_ERR_ALLOC ;
    }
  }
  newBuffer = vl_realloc(*buffer, targetSize) ;
  if (newBuffer) {
    *buffer = newBuffer ;
    *bufferSize = targetSize ;
    return VL_ERR_OK ;
  } else {
    return VL_ERR_ALLOC ;
  }
}

static int
_vl_enlarge_buffer (void ** buffer, vl_size * bufferSize, vl_size targetSize) {
  if (*bufferSize >= targetSize) return VL_ERR_OK ;
  return _vl_resize_buffer(buffer,bufferSize,targetSize) ;
}

/* ---------------------------------------------------------------- */
/*                                            Finding local extrema */
/* ---------------------------------------------------------------- */

/* Todo: make this generally available in the library */

typedef struct _VlCovDetExtremum2
{
  vl_index xi ;
  vl_index yi ;
  float x ;
  float y ;
  float peakScore ;
  float edgeScore ;
} VlCovDetExtremum2 ;

typedef struct _VlCovDetExtremum3
{
  vl_index xi ;
  vl_index yi ;
  vl_index zi ;
  float x ;
  float y ;
  float z ;
  float peakScore ;
  float edgeScore ;
} VlCovDetExtremum3 ;

VL_EXPORT vl_size
vl_find_local_extrema_3 (vl_index ** extrema, vl_size * bufferSize,
                         float const * map,
                         vl_size width, vl_size height, vl_size depth,
                         double threshold) ;

VL_EXPORT vl_size
vl_find_local_extrema_2 (vl_index ** extrema, vl_size * bufferSize,
                         float const * map,
                         vl_size width, vl_size height,
                         double threshold) ;

VL_EXPORT vl_bool
vl_refine_local_extreum_3 (VlCovDetExtremum3 * refined,
                           float const * map,
                           vl_size width, vl_size height, vl_size depth,
                           vl_index x, vl_index y, vl_index z) ;

VL_EXPORT vl_bool
vl_refine_local_extreum_2 (VlCovDetExtremum2 * refined,
                           float const * map,
                           vl_size width, vl_size height,
                           vl_index x, vl_index y) ;

vl_size
vl_find_local_extrema_3 (vl_index ** extrema, vl_size * bufferSize,
                         float const * map,
                         vl_size width, vl_size height, vl_size depth,
                         double threshold)
{
  vl_index x, y, z ;
  vl_size const xo = 1 ;
  vl_size const yo = width ;
  vl_size const zo = width * height ;
  float const *pt = map + xo + yo + zo ;

  vl_size numExtrema = 0 ;
  vl_size requiredSize = 0 ;

#define CHECK_NEIGHBORS_3(v,CMP,SGN)     (\
v CMP ## = SGN threshold &&               \
v CMP *(pt + xo) &&                       \
v CMP *(pt - xo) &&                       \
v CMP *(pt + zo) &&                       \
v CMP *(pt - zo) &&                       \
v CMP *(pt + yo) &&                       \
v CMP *(pt - yo) &&                       \
\
v CMP *(pt + yo + xo) &&                  \
v CMP *(pt + yo - xo) &&                  \
v CMP *(pt - yo + xo) &&                  \
v CMP *(pt - yo - xo) &&                  \
\
v CMP *(pt + xo      + zo) &&             \
v CMP *(pt - xo      + zo) &&             \
v CMP *(pt + yo      + zo) &&             \
v CMP *(pt - yo      + zo) &&             \
v CMP *(pt + yo + xo + zo) &&             \
v CMP *(pt + yo - xo + zo) &&             \
v CMP *(pt - yo + xo + zo) &&             \
v CMP *(pt - yo - xo + zo) &&             \
\
v CMP *(pt + xo      - zo) &&             \
v CMP *(pt - xo      - zo) &&             \
v CMP *(pt + yo      - zo) &&             \
v CMP *(pt - yo      - zo) &&             \
v CMP *(pt + yo + xo - zo) &&             \
v CMP *(pt + yo - xo - zo) &&             \
v CMP *(pt - yo + xo - zo) &&             \
v CMP *(pt - yo - xo - zo) )

  for (z = 1 ; z < (signed)depth - 1 ; ++z) {
    for (y = 1 ; y < (signed)height - 1 ; ++y) {
      for (x = 1 ; x < (signed)width - 1 ; ++x) {
        float value = *pt ;
        if (CHECK_NEIGHBORS_3(value,>,+) || CHECK_NEIGHBORS_3(value,<,-)) {
          numExtrema ++ ;
          requiredSize += sizeof(vl_index) * 3 ;
          if (*bufferSize < requiredSize) {
            int err = _vl_resize_buffer((void**)extrema, bufferSize,
                                        requiredSize + 2000 * 3 * sizeof(vl_index)) ;
            if (err != VL_ERR_OK) abort() ;
          }
          (*extrema) [3 * (numExtrema - 1) + 0] = x ;
          (*extrema) [3 * (numExtrema - 1) + 1] = y ;
          (*extrema) [3 * (numExtrema - 1) + 2] = z ;
        }
        pt += xo ;
      }
      pt += 2*xo ;
    }
    pt += 2*yo ;
  }
  return numExtrema ;
}

vl_size
vl_find_local_extrema_2 (vl_index ** extrema, vl_size * bufferSize,
                         float const* map,
                         vl_size width, vl_size height,
                         double threshold)
{
  vl_index x, y ;
  vl_size const xo = 1 ;
  vl_size const yo = width ;
  float const *pt = map + xo + yo ;

  vl_size numExtrema = 0 ;
  vl_size requiredSize = 0 ;
#define CHECK_NEIGHBORS_2(v,CMP,SGN)     (\
v CMP ## = SGN threshold &&               \
v CMP *(pt + xo) &&                       \
v CMP *(pt - xo) &&                       \
v CMP *(pt + yo) &&                       \
v CMP *(pt - yo) &&                       \
\
v CMP *(pt + yo + xo) &&                  \
v CMP *(pt + yo - xo) &&                  \
v CMP *(pt - yo + xo) &&                  \
v CMP *(pt - yo - xo) )

  for (y = 1 ; y < (signed)height - 1 ; ++y) {
    for (x = 1 ; x < (signed)width - 1 ; ++x) {
      float value = *pt ;
      if (CHECK_NEIGHBORS_2(value,>,+) || CHECK_NEIGHBORS_2(value,<,-)) {
        numExtrema ++ ;
        requiredSize += sizeof(vl_index) * 2 ;
        if (*bufferSize < requiredSize) {
          int err = _vl_resize_buffer((void**)extrema, bufferSize,
                                      requiredSize + 2000 * 2 * sizeof(vl_index)) ;
          if (err != VL_ERR_OK) abort() ;
        }
        (*extrema) [2 * (numExtrema - 1) + 0] = x ;
        (*extrema) [2 * (numExtrema - 1) + 1] = y ;
      }
      pt += xo ;
    }
    pt += 2*xo ;
  }
  return numExtrema ;
}

VL_EXPORT vl_bool
vl_refine_local_extreum_3 (VlCovDetExtremum3 * refined,
                           float const * map,
                           vl_size width, vl_size height, vl_size depth,
                           vl_index x, vl_index y, vl_index z)
{
  vl_size const xo = 1 ;
  vl_size const yo = width ;
  vl_size const zo = width * height ;

  double Dx=0,Dy=0,Dz=0,Dxx=0,Dyy=0,Dzz=0,Dxy=0,Dxz=0,Dyz=0 ;
  double A [3*3], b [3] ;

#define at(dx,dy,dz) (*(pt + (dx)*xo + (dy)*yo + (dz)*zo))
#define Aat(i,j) (A[(i)+(j)*3])

  float const * pt ;
  vl_index dx = 0 ;
  vl_index dy = 0 ;
  /*vl_index dz = 0 ;*/
  vl_index iter ;
  int err ;

  assert (map) ;
  assert (1 <= x && x <= (signed)width - 2) ;
  assert (1 <= y && y <= (signed)height - 2) ;
  assert (1 <= z && z <= (signed)depth - 2) ;

  for (iter = 0 ; iter < 5 ; ++iter) {
    x += dx ;
    y += dy ;
    pt = map + x*xo + y*yo + z*zo ;

    /* compute the gradient */
    Dx = 0.5 * (at(+1,0,0) - at(-1,0,0)) ;
    Dy = 0.5 * (at(0,+1,0) - at(0,-1,0));
    Dz = 0.5 * (at(0,0,+1) - at(0,0,-1)) ;

    /* compute the Hessian */
    Dxx = (at(+1,0,0) + at(-1,0,0) - 2.0 * at(0,0,0)) ;
    Dyy = (at(0,+1,0) + at(0,-1,0) - 2.0 * at(0,0,0)) ;
    Dzz = (at(0,0,+1) + at(0,0,-1) - 2.0 * at(0,0,0)) ;

    Dxy = 0.25 * (at(+1,+1,0) + at(-1,-1,0) - at(-1,+1,0) - at(+1,-1,0)) ;
    Dxz = 0.25 * (at(+1,0,+1) + at(-1,0,-1) - at(-1,0,+1) - at(+1,0,-1)) ;
    Dyz = 0.25 * (at(0,+1,+1) + at(0,-1,-1) - at(0,-1,+1) - at(0,+1,-1)) ;

    /* solve linear system */
    Aat(0,0) = Dxx ;
    Aat(1,1) = Dyy ;
    Aat(2,2) = Dzz ;
    Aat(0,1) = Aat(1,0) = Dxy ;
    Aat(0,2) = Aat(2,0) = Dxz ;
    Aat(1,2) = Aat(2,1) = Dyz ;

    b[0] = - Dx ;
    b[1] = - Dy ;
    b[2] = - Dz ;

    err = vl_solve_linear_system_3(b, A, b) ;

    if (err != VL_ERR_OK) {
      b[0] = 0 ;
      b[1] = 0 ;
      b[2] = 0 ;
      break ;
    }

    /* Keep going if there is sufficient translation */

    dx = (b[0] > 0.6 && x < (signed)width - 2 ?  1 : 0)
    + (b[0] < -0.6 && x > 1 ? -1 : 0) ;

    dy = (b[1] > 0.6 && y < (signed)height - 2 ?  1 : 0)
    + (b[1] < -0.6 && y > 1 ? -1 : 0) ;

    if (dx == 0 && dy == 0) break ;
  }

  /* check threshold and other conditions */
  {
    double peakScore = at(0,0,0)
    + 0.5 * (Dx * b[0] + Dy * b[1] + Dz * b[2]) ;
    double alpha = (Dxx+Dyy)*(Dxx+Dyy) / (Dxx*Dyy - Dxy*Dxy) ;
    double edgeScore ;

    if (alpha < 0) {
      /* not an extremum */
      edgeScore = VL_INFINITY_D ;
    } else {
      edgeScore = (0.5*alpha - 1) + sqrt(VL_MAX(0.25*alpha - 1,0)*alpha) ;
    }

    refined->xi = x ;
    refined->yi = y ;
    refined->zi = z ;
    refined->x = x + b[0] ;
    refined->y = y + b[1] ;
    refined->z = z + b[2] ;
    refined->peakScore = peakScore ;
    refined->edgeScore = edgeScore ;

    return
    err == VL_ERR_OK &&
    vl_abs_d(b[0]) < 1.5 &&
    vl_abs_d(b[1]) < 1.5 &&
    vl_abs_d(b[2]) < 1.5 &&
    0 <= refined->x && refined->x <= (signed)width - 1 &&
    0 <= refined->y && refined->y <= (signed)height - 1 &&
    0 <= refined->z && refined->z <= (signed)depth - 1 ;
  }
#undef Aat
#undef at
}

VL_EXPORT vl_bool
vl_refine_local_extreum_2 (VlCovDetExtremum2 * refined,
                           float const * map,
                           vl_size width, vl_size height,
                           vl_index x, vl_index y)
{
  vl_size const xo = 1 ;
  vl_size const yo = width ;

  double Dx=0,Dy=0,Dxx=0,Dyy=0,Dxy=0;
  double A [2*2], b [2] ;

#define at(dx,dy) (*(pt + (dx)*xo + (dy)*yo ))
#define Aat(i,j) (A[(i)+(j)*2])

  float const * pt ;
  vl_index dx = 0 ;
  vl_index dy = 0 ;
  vl_index iter ;
  int err ;

  assert (map) ;
  assert (1 <= x && x <= (signed)width - 2) ;
  assert (1 <= y && y <= (signed)height - 2) ;

  for (iter = 0 ; iter < 5 ; ++iter) {
    x += dx ;
    y += dy ;
    pt = map + x*xo + y*yo  ;

    /* compute the gradient */
    Dx = 0.5 * (at(+1,0) - at(-1,0)) ;
    Dy = 0.5 * (at(0,+1) - at(0,-1));

    /* compute the Hessian */
    Dxx = (at(+1,0) + at(-1,0) - 2.0 * at(0,0)) ;
    Dyy = (at(0,+1) + at(0,-1) - 2.0 * at(0,0)) ;
    Dxy = 0.25 * (at(+1,+1) + at(-1,-1) - at(-1,+1) - at(+1,-1)) ;

    /* solve linear system */
    Aat(0,0) = Dxx ;
    Aat(1,1) = Dyy ;
    Aat(0,1) = Aat(1,0) = Dxy ;

    b[0] = - Dx ;
    b[1] = - Dy ;

    err = vl_solve_linear_system_2(b, A, b) ;

    if (err != VL_ERR_OK) {
      b[0] = 0 ;
      b[1] = 0 ;
      break ;
    }

    /* Keep going if there is sufficient translation */

    dx = (b[0] > 0.6 && x < (signed)width - 2 ?  1 : 0)
    + (b[0] < -0.6 && x > 1 ? -1 : 0) ;

    dy = (b[1] > 0.6 && y < (signed)height - 2 ?  1 : 0)
    + (b[1] < -0.6 && y > 1 ? -1 : 0) ;

    if (dx == 0 && dy == 0) break ;
  }

  /* check threshold and other conditions */
  {
    double peakScore = at(0,0) + 0.5 * (Dx * b[0] + Dy * b[1]) ;
    double alpha = (Dxx+Dyy)*(Dxx+Dyy) / (Dxx*Dyy - Dxy*Dxy) ;
    double edgeScore ;

    if (alpha < 0) {
      /* not an extremum */
      edgeScore = VL_INFINITY_D ;
    } else {
      edgeScore = (0.5*alpha - 1) + sqrt(VL_MAX(0.25*alpha - 1,0)*alpha) ;
    }

    refined->xi = x ;
    refined->yi = y ;
    refined->x = x + b[0] ;
    refined->y = y + b[1] ;
    refined->peakScore = peakScore ;
    refined->edgeScore = edgeScore ;

    return
    err == VL_ERR_OK &&
    vl_abs_d(b[0]) < 1.5 &&
    vl_abs_d(b[1]) < 1.5 &&
    0 <= refined->x && refined->x <= (signed)width - 1 &&
    0 <= refined->y && refined->y <= (signed)height - 1 ;
  }
#undef Aat
#undef at
}

/* ---------------------------------------------------------------- */
/*                                                Covarant detector */
/* ---------------------------------------------------------------- */

#define VL_COVDET_MAX_NUM_ORIENTATIONS 4
#define VL_COVDET_MAX_NUM_LAPLACIAN_SCALES 4
#define VL_COVDET_AA_PATCH_RESOLUTION 20
#define VL_COVDET_AA_MAX_NUM_ITERATIONS 15
#define VL_COVDET_OR_NUM_ORIENTATION_HISTOGAM_BINS 36
#define VL_COVDET_AA_RELATIVE_INTEGRATION_SIGMA 3
#define VL_COVDET_AA_RELATIVE_DERIVATIVE_SIGMA 1
#define VL_COVDET_AA_MAX_ANISOTROPY 5
#define VL_COVDET_AA_CONVERGENCE_THRESHOLD 1.001
#define VL_COVDET_AA_ACCURATE_SMOOTHING VL_FALSE
#define VL_COVDET_AA_PATCH_EXTENT (3*VL_COVDET_AA_RELATIVE_INTEGRATION_SIGMA)
#define VL_COVDET_OR_ADDITIONAL_PEAKS_RELATIVE_SIZE 0.8
#define VL_COVDET_LAP_NUM_LEVELS 10
#define VL_COVDET_LAP_PATCH_RESOLUTION 16
#define VL_COVDET_LAP_DEF_PEAK_THRESHOLD 0.01
#define VL_COVDET_DOG_DEF_PEAK_THRESHOLD VL_COVDET_LAP_DEF_PEAK_THRESHOLD
#define VL_COVDET_DOG_DEF_EDGE_THRESHOLD 10.0
#define VL_COVDET_HARRIS_DEF_PEAK_THRESHOLD 0.000002
#define VL_COVDET_HARRIS_DEF_EDGE_THRESHOLD 10.0
#define VL_COVDET_HESSIAN_DEF_PEAK_THRESHOLD 0.003
#define VL_COVDET_HESSIAN_DEF_EDGE_THRESHOLD 10.0

struct _VlCovDet
{
  VlScaleSpace *gss ;        
  VlScaleSpace *css ;        
  VlCovDetMethod method ;    
  double peakThreshold ;     
  double edgeThreshold ;     
  double lapPeakThreshold;   
  vl_size octaveResolution ; 
  vl_index firstOctave ;     
  double nonExtremaSuppression ;
  vl_size numNonExtremaSuppressed ;

  VlCovDetFeature *features ;
  vl_size numFeatures ;
  vl_size numFeatureBufferSize ;

  float * patch ;
  vl_size patchBufferSize ;

  vl_bool transposed ;
  VlCovDetFeatureOrientation orientations [VL_COVDET_MAX_NUM_ORIENTATIONS] ;
  VlCovDetFeatureLaplacianScale scales [VL_COVDET_MAX_NUM_LAPLACIAN_SCALES] ;

  vl_bool aaAccurateSmoothing ;
  float aaPatch [(2*VL_COVDET_AA_PATCH_RESOLUTION+1)*(2*VL_COVDET_AA_PATCH_RESOLUTION+1)] ;
  float aaPatchX [(2*VL_COVDET_AA_PATCH_RESOLUTION+1)*(2*VL_COVDET_AA_PATCH_RESOLUTION+1)] ;
  float aaPatchY [(2*VL_COVDET_AA_PATCH_RESOLUTION+1)*(2*VL_COVDET_AA_PATCH_RESOLUTION+1)] ;
  float aaMask [(2*VL_COVDET_AA_PATCH_RESOLUTION+1)*(2*VL_COVDET_AA_PATCH_RESOLUTION+1)] ;

  float lapPatch [(2*VL_COVDET_LAP_PATCH_RESOLUTION+1)*(2*VL_COVDET_LAP_PATCH_RESOLUTION+1)] ;
  float laplacians [(2*VL_COVDET_LAP_PATCH_RESOLUTION+1)*(2*VL_COVDET_LAP_PATCH_RESOLUTION+1)*VL_COVDET_LAP_NUM_LEVELS] ;
  vl_size numFeaturesWithNumScales [VL_COVDET_MAX_NUM_LAPLACIAN_SCALES + 1] ;
}  ;

VlEnumerator vlCovdetMethods [VL_COVDET_METHOD_NUM] = {
  {"DoG" ,              (vl_index)VL_COVDET_METHOD_DOG               },
  {"Hessian",           (vl_index)VL_COVDET_METHOD_HESSIAN           },
  {"HessianLaplace",    (vl_index)VL_COVDET_METHOD_HESSIAN_LAPLACE   },
  {"HarrisLaplace",     (vl_index)VL_COVDET_METHOD_HARRIS_LAPLACE    },
  {"MultiscaleHessian", (vl_index)VL_COVDET_METHOD_MULTISCALE_HESSIAN},
  {"MultiscaleHarris",  (vl_index)VL_COVDET_METHOD_MULTISCALE_HARRIS },
  {0,                   0                                            }
} ;

VlCovDet *
vl_covdet_new (VlCovDetMethod method)
{
  VlCovDet * self = vl_calloc(sizeof(VlCovDet),1) ;
  self->method = method ;
  self->octaveResolution = 3 ;
  self->firstOctave = -1 ;
  switch (self->method) {
    case VL_COVDET_METHOD_DOG :
      self->peakThreshold = VL_COVDET_DOG_DEF_PEAK_THRESHOLD ;
      self->edgeThreshold = VL_COVDET_DOG_DEF_EDGE_THRESHOLD ;
      self->lapPeakThreshold = 0  ; /* not used */
      break ;
    case VL_COVDET_METHOD_HARRIS_LAPLACE:
    case VL_COVDET_METHOD_MULTISCALE_HARRIS:
      self->peakThreshold = VL_COVDET_HARRIS_DEF_PEAK_THRESHOLD ;
      self->edgeThreshold = VL_COVDET_HARRIS_DEF_EDGE_THRESHOLD ;
      self->lapPeakThreshold = VL_COVDET_LAP_DEF_PEAK_THRESHOLD ;
      break ;
    case VL_COVDET_METHOD_HESSIAN :
    case VL_COVDET_METHOD_HESSIAN_LAPLACE:
    case VL_COVDET_METHOD_MULTISCALE_HESSIAN:
      self->peakThreshold = VL_COVDET_HESSIAN_DEF_PEAK_THRESHOLD ;
      self->edgeThreshold = VL_COVDET_HESSIAN_DEF_EDGE_THRESHOLD ;
      self->lapPeakThreshold = VL_COVDET_LAP_DEF_PEAK_THRESHOLD ;
      break;
    default:
      assert(0) ;
  }

  self->nonExtremaSuppression = 0.5 ;
  self->features = NULL ;
  self->numFeatures = 0 ;
  self->numFeatureBufferSize = 0 ;
  self->patch = NULL ;
  self->patchBufferSize = 0 ;
  self->transposed = VL_FALSE ;
  self->aaAccurateSmoothing = VL_COVDET_AA_ACCURATE_SMOOTHING ;

  {
    vl_index const w = VL_COVDET_AA_PATCH_RESOLUTION ;
    vl_index i,j ;
    double step = (2.0 * VL_COVDET_AA_PATCH_EXTENT) / (2*w+1) ;
    double sigma = VL_COVDET_AA_RELATIVE_INTEGRATION_SIGMA ;
    for (j = -w ; j <= w ; ++j) {
      for (i = -w ; i <= w ; ++i) {
        double dx = i*step/sigma ;
        double dy = j*step/sigma ;
        self->aaMask[(i+w) + (2*w+1)*(j+w)] = exp(-0.5*(dx*dx+dy*dy)) ;
      }
    }
  }

  {
    /*
     Covers one octave of Laplacian filters, from sigma=1 to sigma=2.
     The spatial sampling step is 0.5.
     */
    vl_index s ;
    for (s = 0 ; s < VL_COVDET_LAP_NUM_LEVELS ; ++s) {
      double sigmaLap = pow(2.0, -0.5 +
                            (double)s / (VL_COVDET_LAP_NUM_LEVELS - 1)) ;
      double const sigmaImage = 1.0 / sqrt(2.0) ;
      double const step = 0.5 * sigmaImage ;
      double const sigmaDelta = sqrt(sigmaLap*sigmaLap - sigmaImage*sigmaImage) ;
      vl_size const w = VL_COVDET_LAP_PATCH_RESOLUTION ;
      vl_size const num = 2 * w + 1  ;
      float * pt = self->laplacians + s * (num * num) ;

      memset(pt, 0, num * num * sizeof(float)) ;

#define at(x,y) pt[(x+w)+(y+w)*(2*w+1)]
      at(0,0) = - 4.0 ;
      at(-1,0) = 1.0 ;
      at(+1,0) = 1.0 ;
      at(0,1) = 1.0 ;
      at(0,-1) = 1.0 ;
#undef at

      vl_imsmooth_f(pt, num,
                    pt, num, num, num,
                    sigmaDelta / step, sigmaDelta / step) ;

#if 0
      {
        char name [200] ;
        snprintf(name, 200, "/tmp/%f-lap.pgm", sigmaDelta) ;
        vl_pgm_write_f(name, pt, num, num) ;
      }
#endif

    }
  }
  return self ;
}

void
vl_covdet_reset (VlCovDet * self)
{
  if (self->features) {
    vl_free(self->features) ;
    self->features = NULL ;
  }
  if (self->css) {
    vl_scalespace_delete(self->css) ;
    self->css = NULL ;
  }
  if (self->gss) {
    vl_scalespace_delete(self->gss) ;
    self->gss = NULL ;
  }
}

void
vl_covdet_delete (VlCovDet * self)
{
  vl_covdet_reset(self) ;
  if (self->patch) vl_free (self->patch) ;
  vl_free(self) ;
}

int
vl_covdet_append_feature (VlCovDet * self, VlCovDetFeature const * feature)
{
  vl_size requiredSize ;
  assert(self) ;
  assert(feature) ;
  self->numFeatures ++ ;
  requiredSize = self->numFeatures * sizeof(VlCovDetFeature) ;
  if (requiredSize > self->numFeatureBufferSize) {
    int err = _vl_resize_buffer((void**)&self->features, &self->numFeatureBufferSize,
                                (self->numFeatures + 1000) * sizeof(VlCovDetFeature)) ;
    if (err) {
      self->numFeatures -- ;
      return err ;
    }
  }
  self->features[self->numFeatures - 1] = *feature ;
  return VL_ERR_OK ;
}

/* ---------------------------------------------------------------- */
/*                                              Process a new image */
/* ---------------------------------------------------------------- */

int
vl_covdet_put_image (VlCovDet * self,
                     float const * image,
                     vl_size width, vl_size height)
{
  vl_size const minOctaveSize = 16 ;
  vl_index lastOctave ;
  vl_index octaveFirstSubdivision ;
  vl_index octaveLastSubdivision ;
  VlScaleSpaceGeometry geom = vl_scalespace_get_default_geometry(width,height) ;

  assert (self) ;
  assert (image) ;
  assert (width >= 1) ;
  assert (height >= 1) ;

  /* (minOctaveSize - 1) 2^lastOctave <= min(width,height) - 1 */
  lastOctave = vl_floor_d(vl_log2_d(VL_MIN((double)width-1,(double)height-1) / (minOctaveSize - 1))) ;

  if (self->method == VL_COVDET_METHOD_DOG) {
    octaveFirstSubdivision = -1 ;
    octaveLastSubdivision = self->octaveResolution + 1 ;
  } else if (self->method == VL_COVDET_METHOD_HESSIAN) {
    octaveFirstSubdivision = -1 ;
    octaveLastSubdivision = self->octaveResolution ;
  } else {
    octaveFirstSubdivision = 0 ;
    octaveLastSubdivision = self->octaveResolution - 1 ;
  }

  geom.width = width ;
  geom.height = height ;
  geom.firstOctave = self->firstOctave ;
  geom.lastOctave = lastOctave ;
  geom.octaveResolution = self->octaveResolution ;
  geom.octaveFirstSubdivision = octaveFirstSubdivision ;
  geom.octaveLastSubdivision = octaveLastSubdivision ;

  if (self->gss == NULL ||
      ! vl_scalespacegeometry_is_equal (geom,
                                        vl_scalespace_get_geometry(self->gss)))
  {
    if (self->gss) vl_scalespace_delete(self->gss) ;
    self->gss = vl_scalespace_new_with_geometry(geom) ;
    if (self->gss == NULL) return VL_ERR_ALLOC ;
  }
  vl_scalespace_put_image(self->gss, image) ;
  return VL_ERR_OK ;
}

/* ---------------------------------------------------------------- */
/*                                              Cornerness measures */
/* ---------------------------------------------------------------- */

static void
_vl_det_hessian_response (float * hessian,
                          float const * image,
                          vl_size width, vl_size height,
                          double step, double sigma)
{
  float factor = (float) pow(sigma/step, 4.0) ;
  vl_index const xo = 1 ; /* x-stride */
  vl_index const yo = width;  /* y-stride */
  vl_size r, c;

  float p11, p12, p13, p21, p22, p23, p31, p32, p33;

  /* setup input pointer to be centered at 0,1 */
  float const *in = image + yo ;

  /* setup output pointer to be centered at 1,1 */
  float *out = hessian + xo + yo;

  /* move 3x3 window and convolve */
  for (r = 1; r < height - 1; ++r)
  {
    /* fill in shift registers at the beginning of the row */
    p11 = in[-yo]; p12 = in[xo - yo];
    p21 = in[  0]; p22 = in[xo     ];
    p31 = in[+yo]; p32 = in[xo + yo];
    /* setup input pointer to (2,1) of the 3x3 square */
    in += 2;
    for (c = 1; c < width - 1; ++c)
    {
      float Lxx, Lyy, Lxy;
      /* fetch remaining values (last column) */
      p13 = in[-yo]; p23 = *in; p33 = in[+yo];

      /* Compute 3x3 Hessian values from pixel differences. */
      Lxx = (-p21 + 2*p22 - p23);
      Lyy = (-p12 + 2*p22 - p32);
      Lxy = ((p11 - p31 - p13 + p33)/4.0f);

      /* normalize and write out */
      *out = (Lxx * Lyy - Lxy * Lxy) * factor ;

      /* move window */
      p11=p12; p12=p13;
      p21=p22; p22=p23;
      p31=p32; p32=p33;

      /* move input/output pointers */
      in++; out++;
    }
    out += 2;
  }

  /* Copy the computed values to borders */
  in = hessian + yo + xo ;
  out = hessian + xo ;

  /* Top row without corners */
  memcpy(out, in, (width - 2)*sizeof(float));
  out--;
  in -= yo;

  /* Left border columns without last row */
  for (r = 0; r < height - 1; r++){
    *out = *in;
    *(out + yo - 1) = *(in + yo - 3);
    in += yo;
    out += yo;
  }

  /* Bottom corners */
  in -= yo;
  *out = *in;
  *(out + yo - 1) = *(in + yo - 3);

  /* Bottom row without corners */
  out++;
  memcpy(out, in, (width - 2)*sizeof(float));
}

static void
_vl_harris_response (float * harris,
                     float const * image,
                     vl_size width, vl_size height,
                     double step, double sigma,
                     double sigmaI, double alpha)
{
  float factor = (float) pow(sigma/step, 4.0) ;
  vl_index k ;

  float * LxLx ;
  float * LyLy ;
  float * LxLy ;

  LxLx = vl_malloc(sizeof(float) * width * height) ;
  LyLy = vl_malloc(sizeof(float) * width * height) ;
  LxLy = vl_malloc(sizeof(float) * width * height) ;

  vl_imgradient_f (LxLx, LyLy, 1, width, image, width, height, width) ;

  for (k = 0 ; k < (signed)(width * height) ; ++k) {
    float dx = LxLx[k] ;
    float dy = LyLy[k] ;
    LxLx[k] = dx*dx ;
    LyLy[k] = dy*dy ;
    LxLy[k] = dx*dy ;
  }

  vl_imsmooth_f(LxLx, width, LxLx, width, height, width,
                sigmaI / step, sigmaI / step) ;

  vl_imsmooth_f(LyLy, width, LyLy, width, height, width,
                sigmaI / step, sigmaI / step) ;

  vl_imsmooth_f(LxLy, width, LxLy, width, height, width,
                sigmaI / step, sigmaI / step) ;

  for (k = 0 ; k < (signed)(width * height) ; ++k) {
    float a = LxLx[k] ;
    float b = LyLy[k] ;
    float c = LxLy[k] ;

    float determinant = a * b - c * c ;
    float trace = a + b ;

    harris[k] = factor * (determinant - alpha * (trace * trace)) ;
  }

  vl_free(LxLy) ;
  vl_free(LyLy) ;
  vl_free(LxLx) ;
}

static void
_vl_dog_response (float * dog,
                  float const * level1,
                  float const * level2,
                  vl_size width, vl_size height)
{
  vl_index k ;
  for (k = 0 ; k < (signed)(width*height) ; ++k) {
    dog[k] = level2[k] - level1[k] ;
  }
}

/* ---------------------------------------------------------------- */
/*                                                  Detect features */
/* ---------------------------------------------------------------- */

void
vl_covdet_detect (VlCovDet * self)
{
  VlScaleSpaceGeometry geom = vl_scalespace_get_geometry(self->gss) ;
  VlScaleSpaceGeometry cgeom ;
  float * levelxx = NULL ;
  float * levelyy = NULL ;
  float * levelxy = NULL ;
  vl_index o, s ;

  assert (self) ;
  assert (self->gss) ;

  /* clear previous detections if any */
  self->numFeatures = 0 ;

  /* prepare buffers ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
  cgeom = geom ;
  if (self->method == VL_COVDET_METHOD_DOG) {
    cgeom.octaveLastSubdivision -= 1 ;
  }
  if (!self->css ||
      !vl_scalespacegeometry_is_equal(cgeom,
                                      vl_scalespace_get_geometry(self->css)))
  {
    if (self->css) vl_scalespace_delete(self->css) ;
    self->css = vl_scalespace_new_with_geometry(cgeom) ;
  }
  if (self->method == VL_COVDET_METHOD_HARRIS_LAPLACE ||
      self->method == VL_COVDET_METHOD_MULTISCALE_HARRIS) {
    VlScaleSpaceOctaveGeometry oct = vl_scalespace_get_octave_geometry(self->gss, geom.firstOctave) ;
    levelxx = vl_malloc(oct.width * oct.height * sizeof(float)) ;
    levelyy = vl_malloc(oct.width * oct.height * sizeof(float)) ;
    levelxy = vl_malloc(oct.width * oct.height * sizeof(float)) ;
  }

  /* compute cornerness ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
  for (o = cgeom.firstOctave ; o <= cgeom.lastOctave ; ++o) {
    VlScaleSpaceOctaveGeometry oct = vl_scalespace_get_octave_geometry(self->css, o) ;

    for (s = cgeom.octaveFirstSubdivision ; s <= cgeom.octaveLastSubdivision ; ++s) {
      float * level = vl_scalespace_get_level(self->gss, o, s) ;
      float * clevel = vl_scalespace_get_level(self->css, o, s) ;
      double sigma = vl_scalespace_get_level_sigma(self->css, o, s) ;
      switch (self->method) {
        case VL_COVDET_METHOD_DOG:
          _vl_dog_response(clevel,
                           vl_scalespace_get_level(self->gss, o, s + 1),
                           level,
                           oct.width, oct.height) ;
          break ;

        case VL_COVDET_METHOD_HARRIS_LAPLACE:
        case VL_COVDET_METHOD_MULTISCALE_HARRIS:
          _vl_harris_response(clevel,
                              level, oct.width, oct.height, oct.step,
                              sigma, 1.4 * sigma, 0.05) ;
          break ;

        case VL_COVDET_METHOD_HESSIAN:
        case VL_COVDET_METHOD_HESSIAN_LAPLACE:
        case VL_COVDET_METHOD_MULTISCALE_HESSIAN:
          _vl_det_hessian_response(clevel, level, oct.width, oct.height, oct.step, sigma) ;
          break ;

        default:
          assert(0) ;
      }
    }
  }

  /* find and refine local maxima ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
  {
    vl_index * extrema = NULL ;
    vl_size extremaBufferSize = 0 ;
    vl_size numExtrema ;
    vl_size index ;
    for (o = cgeom.firstOctave ; o <= cgeom.lastOctave ; ++o) {
      VlScaleSpaceOctaveGeometry octgeom = vl_scalespace_get_octave_geometry(self->css, o) ;
      double step = octgeom.step ;
      vl_size width = octgeom.width ;
      vl_size height = octgeom.height ;
      vl_size depth = cgeom.octaveLastSubdivision - cgeom.octaveFirstSubdivision + 1 ;

      switch (self->method) {
        case VL_COVDET_METHOD_DOG:
        case VL_COVDET_METHOD_HESSIAN:
        {
          /* scale-space extrema */
          float const * octave =
          vl_scalespace_get_level(self->css, o, cgeom.octaveFirstSubdivision) ;
          numExtrema = vl_find_local_extrema_3(&extrema, &extremaBufferSize,
                                               octave, width, height, depth,
                                               0.8 * self->peakThreshold);
          for (index = 0 ; index < numExtrema ; ++index) {
            VlCovDetExtremum3 refined ;
            VlCovDetFeature feature ;
            vl_bool ok ;
            memset(&feature, 0, sizeof(feature)) ;
            ok = vl_refine_local_extreum_3(&refined,
                                           octave, width, height, depth,
                                           extrema[3*index+0],
                                           extrema[3*index+1],
                                           extrema[3*index+2]) ;
            ok &= fabs(refined.peakScore) > self->peakThreshold ;
            ok &= refined.edgeScore < self->edgeThreshold ;
            if (ok) {
              double sigma = cgeom.baseScale *
              pow(2.0, o + (refined.z + cgeom.octaveFirstSubdivision)
                  / cgeom.octaveResolution) ;
              feature.frame.x = refined.x * step ;
              feature.frame.y = refined.y * step ;
              feature.frame.a11 = sigma ;
              feature.frame.a12 = 0.0 ;
              feature.frame.a21 = 0.0 ;
              feature.frame.a22 = sigma ;
              feature.peakScore = refined.peakScore ;
              feature.edgeScore = refined.edgeScore ;
              vl_covdet_append_feature(self, &feature) ;
            }
          }
          break ;
        }

        default:
        {
          for (s = cgeom.octaveFirstSubdivision ; s < cgeom.octaveLastSubdivision ; ++s) {
            /* space extrema */
            float const * level = vl_scalespace_get_level(self->css,o,s) ;
            numExtrema = vl_find_local_extrema_2(&extrema, &extremaBufferSize,
                                                 level,
                                                 width, height,
                                                 0.8 * self->peakThreshold);
            for (index = 0 ; index < numExtrema ; ++index) {
              VlCovDetExtremum2 refined ;
              VlCovDetFeature feature ;
              vl_bool ok ;
              memset(&feature, 0, sizeof(feature)) ;
              ok = vl_refine_local_extreum_2(&refined,
                                             level, width, height,
                                             extrema[2*index+0],
                                             extrema[2*index+1]);
              ok &= fabs(refined.peakScore) > self->peakThreshold ;
              ok &= refined.edgeScore < self->edgeThreshold ;
              if (ok) {
                double sigma = cgeom.baseScale *
                pow(2.0, o + (double)s / cgeom.octaveResolution) ;
                feature.frame.x = refined.x * step ;
                feature.frame.y = refined.y * step ;
                feature.frame.a11 = sigma ;
                feature.frame.a12 = 0.0 ;
                feature.frame.a21 = 0.0 ;
                feature.frame.a22 = sigma ;
                feature.peakScore = refined.peakScore ;
                feature.edgeScore = refined.edgeScore ;
                vl_covdet_append_feature(self, &feature) ;
              }
            }
          }
          break ;
        }
      }
    } /* next octave */

    if (extrema) { vl_free(extrema) ; extrema = 0 ; }
  }

  /* Laplacian scale selection for certain methods */
  switch (self->method) {
    case VL_COVDET_METHOD_HARRIS_LAPLACE :
    case VL_COVDET_METHOD_HESSIAN_LAPLACE :
      vl_covdet_extract_laplacian_scales (self) ;
      break ;
    default:
      break ;
  }

  if (self->nonExtremaSuppression) {
    vl_index i, j ;
    double tol = self->nonExtremaSuppression ;
    self->numNonExtremaSuppressed = 0 ;
    for (i = 0 ; i < (signed)self->numFeatures ; ++i) {
      double x = self->features[i].frame.x ;
      double y = self->features[i].frame.y ;
      double sigma = self->features[i].frame.a11 ;
      double score = self->features[i].peakScore ;

      for (j = 0 ; j < (signed)self->numFeatures ; ++j) {
        double dx_ = self->features[j].frame.x - x ;
        double dy_ = self->features[j].frame.y - y ;
        double sigma_ = self->features[j].frame.a11 ;
        double score_ = self->features[j].peakScore ;
        if (score_ == 0) continue ;
        if (sigma < (1+tol) * sigma_ &&
            sigma_ < (1+tol) * sigma &&
            vl_abs_d(dx_) < tol * sigma &&
            vl_abs_d(dy_) < tol * sigma &&
            vl_abs_d(score) > vl_abs_d(score_)) {
          self->features[j].peakScore = 0 ;
          self->numNonExtremaSuppressed ++ ;
        }
      }
    }
    j = 0 ;
    for (i = 0 ; i < (signed)self->numFeatures ; ++i) {
      VlCovDetFeature feature = self->features[i] ;
      if (self->features[i].peakScore != 0) {
        self->features[j++] = feature ;
      }
    }
    self->numFeatures = j ;
  }

  if (levelxx) vl_free(levelxx) ;
  if (levelyy) vl_free(levelyy) ;
  if (levelxy) vl_free(levelxy) ;
}

/* ---------------------------------------------------------------- */
/*                                                  Extract patches */
/* ---------------------------------------------------------------- */

vl_bool
vl_covdet_extract_patch_helper (VlCovDet * self,
                                double * sigma1,
                                double * sigma2,
                                float * patch,
                                vl_size resolution,
                                double extent,
                                double sigma,
                                double A_ [4],
                                double T_ [2],
                                double d1, double d2)
{
  vl_index o, s ;
  double factor ;
  double sigma_ ;
  float const * level ;
  vl_size width, height ;
  double step ;

  double A [4] = {A_[0], A_[1], A_[2], A_[3]} ;
  double T [2] = {T_[0], T_[1]} ;

  VlScaleSpaceGeometry geom = vl_scalespace_get_geometry(self->gss) ;
  VlScaleSpaceOctaveGeometry oct ;

  /* Starting from a pre-smoothed image at scale sigma_
     because of the mapping A the resulting smoothing in
     the warped patch is S, where

        sigma_^2 I = A S A',

        S = sigma_^2 inv(A) inv(A)' = sigma_^2 V D^-2 V',

        A = U D V'.

     Thus we rotate A by V to obtain an axis-aligned smoothing:

        A = U*D,

        S = sigma_^2 D^-2.

     Then we search the scale-space for the best sigma_ such
     that the target smoothing is approximated from below:

        max sigma_(o,s) :    simga_(o,s) factor <= sigma,
        factor = max{abs(D11), abs(D22)}.
   */


  /*
   Determine the best level (o,s) such that sigma_(o,s) factor <= sigma.
   This can be obtained by scanning octaves from smallest to largest
   and stopping when no level in the octave satisfies the relation.

   Given the range of octave availables, do the best you can.
   */

  factor = 1.0 / VL_MIN(d1, d2) ;

  for (o = geom.firstOctave + 1 ; o <= geom.lastOctave ; ++o) {
    s = vl_floor_d(vl_log2_d(sigma / (factor * geom.baseScale)) - o) ;
    s = VL_MAX(s, geom.octaveFirstSubdivision) ;
    s = VL_MIN(s, geom.octaveLastSubdivision) ;
    sigma_ = geom.baseScale * pow(2.0, o + (double)s / geom.octaveResolution) ;
    /*VL_PRINTF(".. %d D=%g %g; sigma_=%g factor*sigma_=%g\n", o, d1, d2, sigma_, factor* sigma_) ;*/
    if (factor * sigma_ > sigma) {
      o -- ;
      break ;
    }
  }
  o = VL_MIN(o, geom.lastOctave) ;
  s = vl_floor_d(vl_log2_d(sigma / (factor * geom.baseScale)) - o) ;
  s = VL_MAX(s, geom.octaveFirstSubdivision) ;
  s = VL_MIN(s, geom.octaveLastSubdivision) ;
  sigma_ = geom.baseScale * pow(2.0, o + (double)s / geom.octaveResolution) ;
  if (sigma1) *sigma1 = sigma_ / d1 ;
  if (sigma2) *sigma2 = sigma_ / d2 ;

  /*VL_PRINTF("%d %d %g %g %g %g\n", o, s, factor, sigma_, factor * sigma_, sigma) ;*/

  /*
   Now the scale space level to be used for this warping has been
   determined.

   If the patch is partially or completely out of the image boundary,
   create a padded copy of the required region first.
   */

  level = vl_scalespace_get_level(self->gss, o, s) ;
  oct = vl_scalespace_get_octave_geometry(self->gss, o) ;
  width = oct.width ;
  height = oct.height ;
  step = oct.step ;

  A[0] /= step ;
  A[1] /= step ;
  A[2] /= step ;
  A[3] /= step ;
  T[0] /= step ;
  T[1] /= step ;

  {
    /*
     Warp the patch domain [x0hat,y0hat,x1hat,y1hat] to the image domain/
     Obtain a box [x0,y0,x1,y1] enclosing that wrapped box, and then
     an integer vertexes version [x0i, y0i, x1i, y1i], making room
     for one pixel at the boundary to simplify bilinear interpolation
     later on.
     */
    vl_index x0i, y0i, x1i, y1i ;
    double x0 = +VL_INFINITY_D ;
    double x1 = -VL_INFINITY_D ;
    double y0 = +VL_INFINITY_D ;
    double y1 = -VL_INFINITY_D ;
    double boxx [4] = {extent, extent, -extent, -extent} ;
    double boxy [4] = {-extent, extent, extent, -extent} ;
    int i ;
    for (i = 0 ; i < 4 ; ++i) {
      double x = A[0] * boxx[i] + A[2] * boxy[i] + T[0] ;
      double y = A[1] * boxx[i] + A[3] * boxy[i] + T[1] ;
      x0 = VL_MIN(x0, x) ;
      x1 = VL_MAX(x1, x) ;
      y0 = VL_MIN(y0, y) ;
      y1 = VL_MAX(y1, y) ;
    }

    /* Leave one pixel border for bilinear interpolation. */
    x0i = floor(x0) - 1 ;
    y0i = floor(y0) - 1 ;
    x1i = ceil(x1) + 1 ;
    y1i = ceil(y1) + 1 ;

    /*
     If the box [x0i,y0i,x1i,y1i] is not fully contained in the
     image domain, then create a copy of this region by padding
     the image. The image is extended by continuity.
     */

    if (x0i < 0 || x1i > (signed)width-1 ||
        y0i < 0 || y1i > (signed)height-1) {
      vl_index xi, yi ;

      /* compute the amount of l,r,t,b padding needed to complete the patch */
      vl_index padx0 = VL_MAX(0, - x0i) ;
      vl_index pady0 = VL_MAX(0, - y0i) ;
      vl_index padx1 = VL_MAX(0, x1i - ((signed)width - 1)) ;
      vl_index pady1 = VL_MAX(0, y1i - ((signed)height - 1)) ;

      /* make enough room for the patch */
      vl_index patchWidth = x1i - x0i + 1 ;
      vl_index patchHeight = y1i - y0i + 1 ;
      vl_size patchBufferSize = patchWidth * patchHeight * sizeof(float) ;
      if (patchBufferSize > self->patchBufferSize) {
        int err = _vl_resize_buffer((void**)&self->patch, &self->patchBufferSize, patchBufferSize) ;
        if (err) return vl_set_last_error(VL_ERR_ALLOC, NULL) ;
      }

      if (pady0 < patchHeight - pady1) {
        /* start by filling the central horizontal band */
        for (yi = y0i + pady0 ; yi < y0i + patchHeight - pady1 ; ++ yi) {
          float *dst = self->patch + (yi - y0i) * patchWidth ;
          float const *src = level + yi * width + VL_MIN(VL_MAX(0, x0i),(signed)width-1) ;
          for (xi = x0i ; xi < x0i + padx0 ; ++xi) *dst++ = *src ;
          for ( ; xi < x0i + patchWidth - padx1 - 2 ; ++xi) *dst++ = *src++ ;
          for ( ; xi < x0i + patchWidth ; ++xi) *dst++ = *src ;
        }
        /* now extend the central band up and down */
        for (yi = 0 ; yi < pady0 ; ++yi) {
          memcpy(self->patch + yi * patchWidth,
                 self->patch + pady0 * patchWidth,
                 patchWidth * sizeof(float)) ;
        }
        for (yi = patchHeight - pady1 ; yi < patchHeight ; ++yi) {
          memcpy(self->patch + yi * patchWidth,
                 self->patch + (patchHeight - pady1 - 1) * patchWidth,
                 patchWidth * sizeof(float)) ;
        }
      } else {
        /* should be handled better! */
        memset(self->patch, 0, self->patchBufferSize) ;
      }
#if 0
      {
        char name [200] ;
        snprintf(name, 200, "/tmp/%20.0f-ext.pgm", 1e10*vl_get_cpu_time()) ;
        vl_pgm_write_f(name, patch, patchWidth, patchWidth) ;
      }
#endif

      level = self->patch ;
      width = patchWidth ;
      height = patchHeight ;
      T[0] -= x0i ;
      T[1] -= y0i ;
    }
  }

  /*
   Resample by using bilinear interpolation.
   */
  {
    float * pt = patch ;
    double yhat = -extent ;
    vl_index xxi ;
    vl_index yyi ;
    double stephat = extent / resolution ;

    for (yyi = 0 ; yyi < 2 * (signed)resolution + 1 ; ++yyi) {
      double xhat = -extent ;
      double rx = A[2] * yhat + T[0] ;
      double ry = A[3] * yhat + T[1] ;
      for (xxi = 0 ; xxi < 2 * (signed)resolution + 1 ; ++xxi) {
        double x = A[0] * xhat + rx ;
        double y = A[1] * xhat + ry ;
        vl_index xi = vl_floor_d(x) ;
        vl_index yi = vl_floor_d(y) ;
        double i00 = level[yi * width + xi] ;
        double i10 = level[yi * width + xi + 1] ;
        double i01 = level[(yi + 1) * width + xi] ;
        double i11 = level[(yi + 1) * width + xi + 1] ;
        double wx = x - xi ;
        double wy = y - yi ;

        assert(xi >= 0 && xi <= (signed)width - 1) ;
        assert(yi >= 0 && yi <= (signed)height - 1) ;

        *pt++ =
        (1.0 - wy) * ((1.0 - wx) * i00 + wx * i10) +
        wy * ((1.0 - wx) * i01 + wx * i11) ;

        xhat += stephat ;
      }
      yhat += stephat ;
    }
  }
#if 0
    {
      char name [200] ;
      snprintf(name, 200, "/tmp/%20.0f.pgm", 1e10*vl_get_cpu_time()) ;
      vl_pgm_write_f(name, patch, 2*resolution+1, 2*resolution+1) ;
    }
#endif
  return VL_ERR_OK ;
}

vl_bool
vl_covdet_extract_patch_for_frame (VlCovDet * self,
                                   float * patch,
                                   vl_size resolution,
                                   double extent,
                                   double sigma,
                                   VlFrameOrientedEllipse frame)
{
  double A[2*2] = {frame.a11, frame.a21, frame.a12, frame.a22} ;
  double T[2] = {frame.x, frame.y} ;
  double D[4], U[4], V[4] ;

  vl_svd2(D, U, V, A) ;

  return vl_covdet_extract_patch_helper
  (self, NULL, NULL, patch, resolution, extent, sigma, A, T, D[0], D[3]) ;
}

/* ---------------------------------------------------------------- */
/*                                                     Affine shape */
/* ---------------------------------------------------------------- */

int
vl_covdet_extract_affine_shape_for_frame (VlCovDet * self,
                                          VlFrameOrientedEllipse * adapted,
                                          VlFrameOrientedEllipse frame)
{
  vl_index iter = 0 ;

  double A [2*2] = {frame.a11, frame.a21, frame.a12, frame.a22} ;
  double T [2] = {frame.x, frame.y} ;
  double U [2*2] ;
  double V [2*2] ;
  double D [2*2] ;
  double M [2*2] ;
  double P [2*2] ;
  double P_ [2*2] ;
  double Q [2*2] ;
  double sigma1, sigma2 ;
  double sigmaD = VL_COVDET_AA_RELATIVE_DERIVATIVE_SIGMA ;
  double factor ;
  double anisotropy ;
  double referenceScale ;
  vl_size const resolution = VL_COVDET_AA_PATCH_RESOLUTION ;
  vl_size const side = 2*VL_COVDET_AA_PATCH_RESOLUTION + 1 ;
  double const extent = VL_COVDET_AA_PATCH_EXTENT ;


  *adapted = frame ;

  while (1) {
    double lxx = 0, lxy = 0, lyy = 0 ;
    vl_index k ;
    int err ;

    /* A = U D V' */
    vl_svd2(D, U, V, A) ;
    anisotropy = VL_MAX(D[0]/D[3], D[3]/D[0]) ;

    /* VL_PRINTF("anisot: %g\n", anisotropy); */

    if (anisotropy > VL_COVDET_AA_MAX_ANISOTROPY) {
      /* diverged, give up with current solution */
      break ;
    }

    /* make sure that the smallest singluar value stays fixed
       after the first iteration */
    if (iter == 0) {
      referenceScale = VL_MIN(D[0], D[3]) ;
      factor = 1.0 ;
    } else {
      factor = referenceScale / VL_MIN(D[0],D[3]) ;
    }

    D[0] *= factor ;
    D[3] *= factor ;

    A[0] = U[0] * D[0] ;
    A[1] = U[1] * D[0] ;
    A[2] = U[2] * D[3] ;
    A[3] = U[3] * D[3] ;

    adapted->a11 = A[0] ;
    adapted->a21 = A[1] ;
    adapted->a12 = A[2] ;
    adapted->a22 = A[3] ;

    if (++iter >= VL_COVDET_AA_MAX_NUM_ITERATIONS) break ;

    err = vl_covdet_extract_patch_helper(self,
                                         &sigma1, &sigma2,
                                         self->aaPatch,
                                         resolution,
                                         extent,
                                         sigmaD,
                                         A, T, D[0], D[3]) ;
    if (err) return err ;

    if (self->aaAccurateSmoothing ) {
      double deltaSigma1 = sqrt(VL_MAX(sigmaD*sigmaD - sigma1*sigma1,0)) ;
      double deltaSigma2 = sqrt(VL_MAX(sigmaD*sigmaD - sigma2*sigma2,0)) ;
      double stephat = extent / resolution ;
      vl_imsmooth_f(self->aaPatch, side,
                    self->aaPatch, side, side, side,
                    deltaSigma1 / stephat, deltaSigma2 / stephat) ;
    }

    /* compute second moment matrix */
    vl_imgradient_f (self->aaPatchX, self->aaPatchY, 1, side,
                     self->aaPatch, side, side, side) ;

    for (k = 0 ; k < (signed)(side*side) ; ++k) {
      double lx = self->aaPatchX[k] ;
      double ly = self->aaPatchY[k] ;
      lxx += lx * lx * self->aaMask[k] ;
      lyy += ly * ly * self->aaMask[k] ;
      lxy += lx * ly * self->aaMask[k] ;
    }
    M[0] = lxx ;
    M[1] = lxy ;
    M[2] = lxy ;
    M[3] = lyy ;

    if (lxx == 0 || lyy == 0) {
      *adapted = frame ;
      break ;
    }

    /* decompose M = P * Q * P' */
    vl_svd2 (Q, P, P_, M) ;

    /*
     Setting A <- A * dA results in M to change approximatively as

     M --> dA'  M dA = dA' P Q P dA

     To make this proportional to the identity, we set

     dA ~= P Q^1/2

     we also make it so the smallest singular value of A is unchanged.
     */

    if (Q[3]/Q[0] < VL_COVDET_AA_CONVERGENCE_THRESHOLD &&
        Q[0]/Q[3] < VL_COVDET_AA_CONVERGENCE_THRESHOLD) {
      break ;
    }

    {
      double Ap [4] ;
      double q0 = sqrt(Q[0]) ;
      double q1 = sqrt(Q[3]) ;
      Ap[0] = (A[0] * P[0] + A[2] * P[1]) / q0 ;
      Ap[1] = (A[1] * P[0] + A[3] * P[1]) / q0 ;
      Ap[2] = (A[0] * P[2] + A[2] * P[3]) / q1 ;
      Ap[3] = (A[1] * P[2] + A[3] * P[3]) / q1 ;
      memcpy(A,Ap,4*sizeof(double)) ;
    }

  } /* next iteration */

  /*
   Make upright.

   Shape adaptation does not estimate rotation. This is fixed by default
   so that a selected axis is not rotated at all (usually this is the
   vertical axis for upright features). To do so, the frame is rotated
   as follows.
   */
  {
    double A [2*2] = {adapted->a11, adapted->a21, adapted->a12, adapted->a22} ;
    double ref [2] ;
    double ref_ [2] ;
    double angle ;
    double angle_ ;
    double dangle ;
    double r1, r2 ;

    if (self->transposed) {
      /* up is the x axis */
      ref[0] = 1 ;
      ref[1] = 0 ;
    } else {
      /* up is the y axis */
      ref[0] = 0 ;
      ref[1] = 1 ;
    }

    vl_solve_linear_system_2 (ref_, A, ref) ;
    angle = atan2(ref[1], ref[0]) ;
    angle_ = atan2(ref_[1], ref_[0]) ;
    dangle = angle_ - angle ;
    r1 = cos(dangle) ;
    r2 = sin(dangle) ;
    adapted->a11 = + A[0] * r1 + A[2] * r2 ;
    adapted->a21 = + A[1] * r1 + A[3] * r2 ;
    adapted->a12 = - A[0] * r2 + A[2] * r1 ;
    adapted->a22 = - A[1] * r2 + A[3] * r1 ;
  }

  return VL_ERR_OK ;
}

void
vl_covdet_extract_affine_shape (VlCovDet * self)
{
  vl_index i, j = 0 ;
  vl_size numFeatures = vl_covdet_get_num_features(self) ;
  VlCovDetFeature * feature = vl_covdet_get_features(self);
  for (i = 0 ; i < (signed)numFeatures ; ++i) {
    int status ;
    VlFrameOrientedEllipse adapted ;
    status = vl_covdet_extract_affine_shape_for_frame(self, &adapted, feature[i].frame) ;
    if (status == VL_ERR_OK) {
      feature[j] = feature[i] ;
      feature[j].frame = adapted ;
      ++ j ;
    }
  }
  self->numFeatures = j ;
}

/* ---------------------------------------------------------------- */
/*                                                      Orientation */
/* ---------------------------------------------------------------- */

static int
_vl_covdet_compare_orientations_descending (void const * a_,
                                            void const * b_)
{
  VlCovDetFeatureOrientation const * a = a_ ;
  VlCovDetFeatureOrientation const * b = b_ ;
  if (a->score > b->score) return -1 ;
  if (a->score < b->score) return +1 ;
  return 0 ;
}

VlCovDetFeatureOrientation *
vl_covdet_extract_orientations_for_frame (VlCovDet * self,
                                          vl_size * numOrientations,
                                          VlFrameOrientedEllipse frame)
{
  int err ;
  vl_index k, i ;
  vl_index iter ;

  double extent = VL_COVDET_AA_PATCH_EXTENT ;
  vl_size resolution = VL_COVDET_AA_PATCH_RESOLUTION ;
  vl_size side = 2 * resolution + 1  ;

  vl_size const numBins = VL_COVDET_OR_NUM_ORIENTATION_HISTOGAM_BINS ;
  double hist [VL_COVDET_OR_NUM_ORIENTATION_HISTOGAM_BINS] ;
  double const binExtent = 2 * VL_PI / VL_COVDET_OR_NUM_ORIENTATION_HISTOGAM_BINS ;
  double const peakRelativeSize = VL_COVDET_OR_ADDITIONAL_PEAKS_RELATIVE_SIZE ;
  double maxPeakValue ;

  double A [2*2] = {frame.a11, frame.a21, frame.a12, frame.a22} ;
  double T [2] = {frame.x, frame.y} ;
  double U [2*2] ;
  double V [2*2] ;
  double D [2*2] ;
  double sigma1, sigma2 ;
  double sigmaD = 1.0 ;
  double theta0 ;

  assert(self);
  assert(numOrientations) ;

  /*
   The goal is to estimate a rotation R(theta) such that the patch given
   by the transformation A R(theta) has the strongest average
   gradient pointing right (or down for transposed conventions).

   To compensate for tha anisotropic smoothing due to warping,
   A is decomposed as A = U D V' and the patch is warped by
   U D only, meaning that the matrix R_(theta) will be estimated instead,
   where:

      A R(theta) = U D V' R(theta) = U D R_(theta)

   such that R(theta) = V R(theta). That is an extra rotation of
   theta0 = atan2(V(2,1),V(1,1)).
   */

  /* axis aligned anisotropic smoothing for easier compensation */
  vl_svd2(D, U, V, A) ;

  A[0] = U[0] * D[0] ;
  A[1] = U[1] * D[0] ;
  A[2] = U[2] * D[3] ;
  A[3] = U[3] * D[3] ;

  theta0 = atan2(V[1],V[0]) ;

  err = vl_covdet_extract_patch_helper(self,
                                       &sigma1, &sigma2,
                                       self->aaPatch,
                                       resolution,
                                       extent,
                                       sigmaD,
                                       A, T, D[0], D[3]) ;

  if (err) {
    *numOrientations = 0 ;
    return NULL ;
  }

  if (1) {
    double deltaSigma1 = sqrt(VL_MAX(sigmaD*sigmaD - sigma1*sigma1,0)) ;
    double deltaSigma2 = sqrt(VL_MAX(sigmaD*sigmaD - sigma2*sigma2,0)) ;
    double stephat = extent / resolution ;
    vl_imsmooth_f(self->aaPatch, side,
                  self->aaPatch, side, side, side,
                  deltaSigma1 / stephat, deltaSigma2 / stephat) ;
  }

  /* histogram of oriented gradients */
  vl_imgradient_polar_f (self->aaPatchX, self->aaPatchY, 1, side,
                         self->aaPatch, side, side, side) ;

  memset (hist, 0, sizeof(double) * numBins) ;

  for (k = 0 ; k < (signed)(side*side) ; ++k) {
    double modulus = self->aaPatchX[k] ;
    double angle = self->aaPatchY[k] ;
    double weight = self->aaMask[k] ;

    double x = angle / binExtent ;
    vl_index bin = vl_floor_d(x) ;
    double w2 = x - bin ;
    double w1 = 1.0 - w2 ;

    hist[(bin + numBins) % numBins] += w1 * (modulus * weight) ;
    hist[(bin + numBins + 1) % numBins] += w2 * (modulus * weight) ;
  }

  /* smooth histogram */
  for (iter = 0; iter < 6; iter ++) {
    double prev = hist [numBins - 1] ;
    double first = hist [0] ;
    vl_index i ;
    for (i = 0; i < (signed)numBins - 1; ++i) {
      double curr = (prev + hist[i] + hist[(i + 1) % numBins]) / 3.0 ;
      prev = hist[i] ;
      hist[i] = curr ;
    }
    hist[i] = (prev + hist[i] + first) / 3.0 ;
  }

  /* find the histogram maximum */
  maxPeakValue = 0 ;
  for (i = 0 ; i < (signed)numBins ; ++i) {
    maxPeakValue = VL_MAX (maxPeakValue, hist[i]) ;
  }

  /* find peaks within 80% from max */
  *numOrientations = 0 ;
  for(i = 0 ; i < (signed)numBins ; ++i) {
    double h0 = hist [i] ;
    double hm = hist [(i - 1 + numBins) % numBins] ;
    double hp = hist [(i + 1 + numBins) % numBins] ;

    /* is this a peak? */
    if (h0 > peakRelativeSize * maxPeakValue && h0 > hm && h0 > hp) {
      /* quadratic interpolation */
      double di = - 0.5 * (hp - hm) / (hp + hm - 2 * h0) ;
      double th = binExtent * (i + di) + theta0 ;
      if (self->transposed) {
        /* the axis to the right is y, measure orientations from this */
        th = th - VL_PI/2 ;
      }
      self->orientations[*numOrientations].angle = th ;
      self->orientations[*numOrientations].score = h0 ;
      *numOrientations += 1 ;
      //VL_PRINTF("%d %g\n", *numOrientations, th) ;

      if (*numOrientations >= VL_COVDET_MAX_NUM_ORIENTATIONS) break ;
    }
  }

  /* sort the orientations by decreasing scores */
  qsort(self->orientations,
        *numOrientations,
        sizeof(VlCovDetFeatureOrientation),
        _vl_covdet_compare_orientations_descending) ;

  return self->orientations ;
}

void
vl_covdet_extract_orientations (VlCovDet * self)
{
  vl_index i, j  ;
  vl_size numFeatures = vl_covdet_get_num_features(self) ;
  for (i = 0 ; i < (signed)numFeatures ; ++i) {
    vl_size numOrientations ;
    VlCovDetFeature feature = self->features[i] ;
    VlCovDetFeatureOrientation* orientations =
    vl_covdet_extract_orientations_for_frame(self, &numOrientations, feature.frame) ;

    for (j = 0 ; j < (signed)numOrientations ; ++j) {
      double A [2*2] = {
        feature.frame.a11,
        feature.frame.a21,
        feature.frame.a12,
        feature.frame.a22} ;
      double r1 = cos(orientations[j].angle) ;
      double r2 = sin(orientations[j].angle) ;
      VlCovDetFeature * oriented ;

      if (j == 0) {
        oriented = & self->features[i] ;
      } else {
        vl_covdet_append_feature(self, &feature) ;
        oriented = & self->features[self->numFeatures -1] ;
      }

      oriented->orientationScore = orientations[j].score ;
      oriented->frame.a11 = + A[0] * r1 + A[2] * r2 ;
      oriented->frame.a21 = + A[1] * r1 + A[3] * r2 ;
      oriented->frame.a12 = - A[0] * r2 + A[2] * r1 ;
      oriented->frame.a22 = - A[1] * r2 + A[3] * r1 ;
    }
  }
}

/* ---------------------------------------------------------------- */
/*                                                 Laplacian scales */
/* ---------------------------------------------------------------- */

VlCovDetFeatureLaplacianScale *
vl_covdet_extract_laplacian_scales_for_frame (VlCovDet * self,
                                              vl_size * numScales,
                                              VlFrameOrientedEllipse frame)
{
  /*
   We try to explore one octave, with the nominal detection scale 1.0
   (in the patch reference frame) in the middle. Thus the goal is to sample
   the response of the tr-Laplacian operator at logarithmically
   spaced scales in 1/sqrt(2), sqrt(2).

   To this end, the patch is warped with a smoothing of at most
   sigmaImage = 1 / sqrt(2) (beginning of the scale), sampled at
   roughly twice the Nyquist frequency (so step = 1 / (2*sqrt(2))).
   This maes it possible to approximate the Laplacian operator at
   that scale by simple finite differences.

   */
  int err ;
  double const sigmaImage = 1.0 / sqrt(2.0) ;
  double const step = 0.5 * sigmaImage ;
  double actualSigmaImage ;
  vl_size const resolution = VL_COVDET_LAP_PATCH_RESOLUTION ;
  vl_size const num = 2 * resolution + 1 ;
  double extent = step * resolution ;
  double scores [VL_COVDET_LAP_NUM_LEVELS] ;
  double factor = 1.0 ;
  float const * pt ;
  vl_index k ;

  double A[2*2] = {frame.a11, frame.a21, frame.a12, frame.a22} ;
  double T[2] = {frame.x, frame.y} ;
  double D[4], U[4], V[4] ;
  double sigma1, sigma2 ;

  assert(self) ;
  assert(numScales) ;

  *numScales = 0 ;

  vl_svd2(D, U, V, A) ;

  err = vl_covdet_extract_patch_helper
  (self, &sigma1, &sigma2, self->lapPatch, resolution, extent, sigmaImage, A, T, D[0], D[3]) ;
  if (err) return NULL ;

  /* the actual smoothing after warping is never the target one */
  if (sigma1 == sigma2) {
    actualSigmaImage = sigma1 ;
  } else {
    /* here we could compensate */
    actualSigmaImage = sqrt(sigma1*sigma2) ;
  }

  /* now multiply by the bank of Laplacians */
  pt = self->laplacians ;
  for (k = 0 ; k < VL_COVDET_LAP_NUM_LEVELS ; ++k) {
    vl_index q ;
    double score = 0 ;
    double sigmaLap = pow(2.0, -0.5 + (double)k / (VL_COVDET_LAP_NUM_LEVELS - 1)) ;
    /* note that the sqrt argument cannot be negative since by construction
     sigmaLap >= sigmaImage */
    sigmaLap = sqrt(sigmaLap*sigmaLap
                    - sigmaImage*sigmaImage
                    + actualSigmaImage*actualSigmaImage) ;

    for (q = 0 ; q < (signed)(num * num) ; ++q) {
      score += (*pt++) * self->lapPatch[q] ;
    }
    scores[k] = score * sigmaLap * sigmaLap ;
  }

  /* find and interpolate maxima */
  for (k = 1 ; k < VL_COVDET_LAP_NUM_LEVELS - 1 ; ++k) {
    double a = scores[k-1] ;
    double b = scores[k] ;
    double c = scores[k+1] ;
    double t = self->lapPeakThreshold ;

    if ((((b > a) && (b > c)) || ((b < a) && (b < c))) && (vl_abs_d(b) >= t)) {
      double dk = - 0.5 * (c - a) / (c + a - 2 * b) ;
      double s = k + dk ;
      double sigmaLap = pow(2.0, -0.5 + s / (VL_COVDET_LAP_NUM_LEVELS - 1)) ;
      double scale ;
      sigmaLap = sqrt(sigmaLap*sigmaLap
                      - sigmaImage*sigmaImage
                      + actualSigmaImage*actualSigmaImage) ;
      scale = sigmaLap / 1.0 ;
      /*
       VL_PRINTF("** k:%d, s:%f, sigmaLapFilter:%f, sigmaLap%f, scale:%f (%f %f %f)\n",
       k,s,sigmaLapFilter,sigmaLap,scale,a,b,c) ;
       */
      if (*numScales < VL_COVDET_MAX_NUM_LAPLACIAN_SCALES) {
        self->scales[*numScales].scale = scale * factor ;
        self->scales[*numScales].score = b + 0.5 * (c - a) * dk ;
        *numScales += 1 ;
      }
    }
  }
  return self->scales ;
}

void
vl_covdet_extract_laplacian_scales (VlCovDet * self)
{
  vl_index i, j  ;
  vl_bool dropFeaturesWithoutScale = VL_TRUE ;
  vl_size numFeatures = vl_covdet_get_num_features(self) ;
  memset(self->numFeaturesWithNumScales, 0,
         sizeof(self->numFeaturesWithNumScales)) ;

  for (i = 0 ; i < (signed)numFeatures ; ++i) {
    vl_size numScales ;
    VlCovDetFeature feature = self->features[i] ;
    VlCovDetFeatureLaplacianScale const * scales =
    vl_covdet_extract_laplacian_scales_for_frame(self, &numScales, feature.frame) ;

    self->numFeaturesWithNumScales[numScales] ++ ;

    if (numScales == 0 && dropFeaturesWithoutScale) {
      self->features[i].peakScore = 0 ;
    }

    for (j = 0 ; j < (signed)numScales ; ++j) {
      VlCovDetFeature * scaled ;

      if (j == 0) {
        scaled = & self->features[i] ;
      } else {
        vl_covdet_append_feature(self, &feature) ;
        scaled = & self->features[self->numFeatures -1] ;
      }

      scaled->laplacianScaleScore = scales[j].score ;
      scaled->frame.a11 *= scales[j].scale ;
      scaled->frame.a21 *= scales[j].scale ;
      scaled->frame.a12 *= scales[j].scale ;
      scaled->frame.a22 *= scales[j].scale ;
    }
  }
  if (dropFeaturesWithoutScale) {
    j = 0 ;
    for (i = 0 ; i < (signed)self->numFeatures ; ++i) {
      VlCovDetFeature feature = self->features[i] ;
      if (feature.peakScore) {
        self->features[j++] = feature ;
      }
    }
    self->numFeatures = j ;
  }

}

/* ---------------------------------------------------------------- */
/*                       Checking that features are inside an image */
/* ---------------------------------------------------------------- */

vl_bool
_vl_covdet_check_frame_inside (VlCovDet * self, VlFrameOrientedEllipse frame, double margin)
{
  double extent = margin ;
  double A [2*2] = {frame.a11, frame.a21, frame.a12, frame.a22} ;
  double T[2] = {frame.x, frame.y} ;
  double x0 = +VL_INFINITY_D ;
  double x1 = -VL_INFINITY_D ;
  double y0 = +VL_INFINITY_D ;
  double y1 = -VL_INFINITY_D ;
  double boxx [4] = {extent, extent, -extent, -extent} ;
  double boxy [4] = {-extent, extent, extent, -extent} ;
  VlScaleSpaceGeometry geom = vl_scalespace_get_geometry(self->gss) ;
  int i ;
  for (i = 0 ; i < 4 ; ++i) {
    double x = A[0] * boxx[i] + A[2] * boxy[i] + T[0] ;
    double y = A[1] * boxx[i] + A[3] * boxy[i] + T[1] ;
    x0 = VL_MIN(x0, x) ;
    x1 = VL_MAX(x1, x) ;
    y0 = VL_MIN(y0, y) ;
    y1 = VL_MAX(y1, y) ;
  }

  return
  0 <= x0 && x1 <= geom.width-1 &&
  0 <= y0 && y1 <= geom.height-1 ;
}

void
vl_covdet_drop_features_outside (VlCovDet * self, double margin)
{
  vl_index i, j = 0 ;
  vl_size numFeatures = vl_covdet_get_num_features(self) ;
  for (i = 0 ; i < (signed)numFeatures ; ++i) {
    vl_bool inside =
    _vl_covdet_check_frame_inside (self, self->features[i].frame, margin) ;
    if (inside) {
      self->features[j] = self->features[i] ;
      ++j ;
    }
  }
  self->numFeatures = j ;
}

/* ---------------------------------------------------------------- */
/*                                              Setters and getters */
/* ---------------------------------------------------------------- */

/* ---------------------------------------------------------------- */
vl_bool
vl_covdet_get_transposed (VlCovDet const  * self)
{
  return self->transposed ;
}

void
vl_covdet_set_transposed (VlCovDet * self, vl_bool t)
{
  self->transposed = t ;
}

/* ---------------------------------------------------------------- */
double
vl_covdet_get_edge_threshold (VlCovDet const * self)
{
  return self->edgeThreshold ;
}

void
vl_covdet_set_edge_threshold (VlCovDet * self, double edgeThreshold)
{
  assert(edgeThreshold >= 0) ;
  self->edgeThreshold = edgeThreshold ;
}

/* ---------------------------------------------------------------- */
double
vl_covdet_get_peak_threshold (VlCovDet const * self)
{
  return self->peakThreshold ;
}

void
vl_covdet_set_peak_threshold (VlCovDet * self, double peakThreshold)
{
  assert(peakThreshold >= 0) ;
  self->peakThreshold = peakThreshold ;
}

/* ---------------------------------------------------------------- */
double
vl_covdet_get_laplacian_peak_threshold (VlCovDet const * self)
{
  return self->lapPeakThreshold ;
}

void
vl_covdet_set_laplacian_peak_threshold (VlCovDet * self, double peakThreshold)
{
  assert(peakThreshold >= 0) ;
  self->lapPeakThreshold = peakThreshold ;
}

/* ---------------------------------------------------------------- */
vl_index
vl_covdet_get_first_octave (VlCovDet const * self)
{
  return self->firstOctave ;
}

void
vl_covdet_set_first_octave (VlCovDet * self, vl_index o)
{
  self->firstOctave = o ;
  vl_covdet_reset(self) ;
}

/* ---------------------------------------------------------------- */
vl_size
vl_covdet_get_octave_resolution (VlCovDet const * self)
{
  return self->octaveResolution ;
}

void
vl_covdet_set_octave_resolution (VlCovDet * self, vl_size r)
{
  self->octaveResolution = r ;
  vl_covdet_reset(self) ;
}

/* ---------------------------------------------------------------- */
vl_bool
vl_covdet_get_aa_accurate_smoothing (VlCovDet const * self)
{
  return self->aaAccurateSmoothing ;
}

void
vl_covdet_set_aa_accurate_smoothing (VlCovDet * self, vl_bool x)
{
  self->aaAccurateSmoothing = x ;
}

/* ---------------------------------------------------------------- */
double
vl_covdet_get_non_extrema_suppression_threshold (VlCovDet const * self)
{
  return self->nonExtremaSuppression ;
}

void
vl_covdet_set_non_extrema_suppression_threshold (VlCovDet * self, double x)
{
  self->nonExtremaSuppression = x ;
}

vl_size
vl_covdet_get_num_non_extrema_suppressed (VlCovDet const * self)
{
  return self->numNonExtremaSuppressed ;
}


/* ---------------------------------------------------------------- */
vl_size
vl_covdet_get_num_features (VlCovDet const * self)
{
  return self->numFeatures ;
}

void *
vl_covdet_get_features (VlCovDet * self)
{
  return self->features ;
}

VlScaleSpace *
vl_covdet_get_gss (VlCovDet const * self)
{
  return self->gss ;
}

VlScaleSpace *
vl_covdet_get_css (VlCovDet const * self)
{
  return self->css ;
}

vl_size const *
vl_covdet_get_laplacian_scales_statistics (VlCovDet const * self,
                                           vl_size * numScales)
{
  *numScales = VL_COVDET_MAX_NUM_LAPLACIAN_SCALES ;
  return self->numFeaturesWithNumScales ;
}