/*********************************************************** 
*  --- OpenSURF ---                                        *
*  This library is distributed under the GNU GPL. Please   *
*  contact chris.evans@irisys.co.uk for more information.  *
*                                                          *
*  C. Evans, Research Into Robust Visual Features,         *
*  MSc University of Bristol, 2008.                        *
*                                                          *
************************************************************/

#include "surf.h"

// This file has slightly changed for compatibility with the OpenGSURF library
// Among the new capabilities these are:
// 1. New descriptors added (Gauge-SURF, M-SURF and standard based spatial derivatives
// See surf.h or README file for more information about the descriptors
// 2. SIFT clipping normalization style available. Set USE_CLIPPING_NORMALIZATION to 1
// OpenGSURF
// Author: Pablo F. Alcantarilla
// Date: 11/11/2011
// Email: pablofdezalc@gmail.com

//**************************************************************************************
//**************************************************************************************

// SIFT clipping normalization style. Set to 1 in case you want to use 
// this normalization style
const bool USE_CLIPPING_NORMALIZATION = false; 
const float CLIPPING_NORMALIZATION_RATIO = 1.8f;
const int CLIPPING_NORMALIZATION_NITER = 5;

//! SURF priors (these need not be done at runtime)
const float pi = 3.14159f;

const double gauss25 [7][7] = {
  0.02350693969273,0.01849121369071,0.01239503121241,0.00708015417522,0.00344628101733,0.00142945847484,0.00050524879060,
  0.02169964028389,0.01706954162243,0.01144205592615,0.00653580605408,0.00318131834134,0.00131955648461,0.00046640341759,
  0.01706954162243,0.01342737701584,0.00900063997939,0.00514124713667,0.00250251364222,0.00103799989504,0.00036688592278,
  0.01144205592615,0.00900063997939,0.00603330940534,0.00344628101733,0.00167748505986,0.00069579213743,0.00024593098864,
  0.00653580605408,0.00514124713667,0.00344628101733,0.00196854695367,0.00095819467066,0.00039744277546,0.00014047800980,
  0.00318131834134,0.00250251364222,0.00167748505986,0.00095819467066,0.00046640341759,0.00019345616757,0.00006837798818,
  0.00131955648461,0.00103799989504,0.00069579213743,0.00039744277546,0.00019345616757,0.00008024231247,0.00002836202103
};

const double gauss33 [11][11] = {
  0.014614763,0.013958917,0.012162744,0.00966788,0.00701053,0.004637568,0.002798657,0.001540738,0.000773799,0.000354525,0.000148179,
  0.013958917,0.013332502,0.011616933,0.009234028,0.006695928,0.004429455,0.002673066,0.001471597,0.000739074,0.000338616,0.000141529,
  0.012162744,0.011616933,0.010122116,0.008045833,0.005834325,0.003859491,0.002329107,0.001282238,0.000643973,0.000295044,0.000123318,
  0.00966788,0.009234028,0.008045833,0.006395444,0.004637568,0.003067819,0.001851353,0.001019221,0.000511879,0.000234524,9.80224E-05,
  0.00701053,0.006695928,0.005834325,0.004637568,0.003362869,0.002224587,0.001342483,0.000739074,0.000371182,0.000170062,7.10796E-05,
  0.004637568,0.004429455,0.003859491,0.003067819,0.002224587,0.001471597,0.000888072,0.000488908,0.000245542,0.000112498,4.70202E-05,
  0.002798657,0.002673066,0.002329107,0.001851353,0.001342483,0.000888072,0.000535929,0.000295044,0.000148179,6.78899E-05,2.83755E-05,
  0.001540738,0.001471597,0.001282238,0.001019221,0.000739074,0.000488908,0.000295044,0.00016243,8.15765E-05,3.73753E-05,1.56215E-05,
  0.000773799,0.000739074,0.000643973,0.000511879,0.000371182,0.000245542,0.000148179,8.15765E-05,4.09698E-05,1.87708E-05,7.84553E-06,
  0.000354525,0.000338616,0.000295044,0.000234524,0.000170062,0.000112498,6.78899E-05,3.73753E-05,1.87708E-05,8.60008E-06,3.59452E-06,
  0.000148179,0.000141529,0.000123318,9.80224E-05,7.10796E-05,4.70202E-05,2.83755E-05,1.56215E-05,7.84553E-06,3.59452E-06,1.50238E-06
};

//**************************************************************************************
//**************************************************************************************

//! Constructor
Surf::Surf(IplImage *img, IpVec &ipts)
: ipts(ipts)
{
  this->img = img;
  
  // Allocate memory for the SURF descriptor vector
  for( unsigned int i = 0; i < ipts.size(); i++ )
  {
	  ipts[i].Allocate_Descriptor_Memory();
  }
}

//**************************************************************************************
//**************************************************************************************

//! Describe all features in the supplied vector
void Surf::getDescriptors(unsigned char mode, bool upright, bool extended)
{
  struct timeval tod1;
  double t1 = 0.0, t2 = 0.0;
  
  // Check there are Ipoints to be described
  if (!ipts.size()) return;

  // Get the size of the vector for fixed loop bounds
  int ipts_size = (int)ipts.size();

  gettimeofday(&tod1,NULL);
  t1 = tod1.tv_sec + tod1.tv_usec / 1000000.0;
  
  if (upright)
  {
    // U-SURF loop just gets descriptors
    for (int i = 0; i < ipts_size; ++i)
    {
      // Set the Ipoint to be described
      index = i;

      // Extract upright (i.e. not rotation invariant) descriptors
	  // OpenSURF original upright descriptor with Gaussian Weighting (64,128)
      if( mode == 0 )
      {
        if( extended == false )
		{
			getUprightDescriptor_64();	
		}
		else
		{
			getUprightExtendedDescriptor_128();
		}
      }
	  // M-SURF upright descriptor (64,128)
      else if( mode == 1 )
      {
        if( extended == false )
		{
			getMUprightDescriptor_64();
		}
		else
		{
			getMUprightExtendedDescriptor_128();
		}
      }
	  // G-SURF upright descriptor (16,32)
      else if( mode == 2 )
      {
        if( extended == false )
		{
			getGaugeUprightDescriptor_16();
		}
		else
		{
			getGaugeUprightExtendedDescriptor_32();
		}
      }
	  // G-SURF upright descriptor (36,72)
	  else if( mode == 3 )
      {
		if( extended == false )
		{
			getGaugeUprightDescriptor_36();  
		}
		else
		{
			getGaugeUprightExtendedDescriptor_72();  
		}
      }
	  // G-SURF upright descriptor (64,128)
	  else if( mode == 4 )
      {
		if( extended == false )
		{
			getGaugeUprightDescriptor_64();  
		}
		else
		{
			getGaugeUprightExtendedDescriptor_128();  
		}
      }
	  // G-SURF upright descriptor (100,200)
	  else if( mode == 5 )
      {
		if( extended == false )
		{
			getGaugeUprightDescriptor_100();  
		}
		else
		{
			getGaugeUprightExtendedDescriptor_200();  
		}
      }
	  // G-SURF upright descriptor (144,288)
	  else if( mode == 6 )
      {
		if( extended == false )
		{
			getGaugeUprightDescriptor_144();  
		}
		else
		{
			getGaugeUprightExtendedDescriptor_288();  
		}
      }
	  // SURF no Gaussian Weighting upright descriptor (16,32)
	  else if( mode == 7 )
      {
		if( extended == false )
		{
			getUprightDescriptor_NW_16();  
		}
		else
		{
			getUprightExtendedDescriptor_NW_32();   
		}
      }
	  // SURF no Gaussian Weighting upright descriptor (36,72)
	  else if( mode == 8 )
      {
		if( extended == false )
		{
			getUprightDescriptor_NW_36();  
		}
		else
		{
			getUprightExtendedDescriptor_NW_72();   
		}
      }	  
	  // SURF no Gaussian Weighting upright descriptor (64,128)
	  else if( mode == 9 )
      {
		if( extended == false )
		{
			getUprightDescriptor_NW_64();  
		}
		else
		{
			getUprightExtendedDescriptor_NW_128();   
		}
      }	  
	  // SURF no Gaussian Weighting upright descriptor (100,200)
	  else if( mode == 10 )
      {
		if( extended == false )
		{
			getUprightDescriptor_NW_100();  
		}
		else
		{
			getUprightExtendedDescriptor_NW_200();   
		}
      }
	  // SURF no Gaussian Weighting upright descriptor (144,288)
	  else if( mode == 11 )
      {
		if( extended == false )
		{
			getUprightDescriptor_NW_144();  
		}
		else
		{
			getUprightExtendedDescriptor_NW_288();   
		}
      }
	  else if( mode == 12 )
	  {
		  if( extended == false )
		  {
			  getMGaugeUprightDescriptor_64();
		  }
		  else
		  {
			  getMGaugeUprightExtendedDescriptor_128();
		  }
	  }
    }
  }
  else
  {
    // Main SURF-64 loop assigns orientations and gets descriptors
    for (int i = 0; i < ipts_size; ++i)
    {
      // Set the Ipoint to be described
      index = i;

      // Assign Orientations and extract rotation invariant descriptors
      getOrientation();

      // Extract upright (i.e. not rotation invariant) descriptors
	  // OpenSURF original descriptor with Gaussian Weighting (64,128)
      if( mode == 0 )
      {
        if( extended == false )
		{
			getDescriptor_64();	
		}
		else
		{
			getExtendedDescriptor_128();
		}
      }
	  // M-SURF descriptor (64,128)
      else if( mode == 1 )
      {
        if( extended == false )
		{
			getMDescriptor_64();
		}
		else
		{
			getMExtendedDescriptor_128();
		}
      }
	  // G-SURF descriptor (16,32)
      else if( mode == 2 )
      {
        if( extended == false )
		{
			getGaugeDescriptor_16();
		}
		else
		{
			getGaugeExtendedDescriptor_32();
		}
      }
	  // G-SURF descriptor (36,72)
	  else if( mode == 3 )
      {
		if( extended == false )
		{
			getGaugeDescriptor_36();  
		}
		else
		{
			getGaugeExtendedDescriptor_72();  
		}
      }
	  // G-SURF descriptor (64,128)
	  else if( mode == 4 )
      {
		if( extended == false )
		{
			getGaugeDescriptor_64();  
		}
		else
		{
			getGaugeExtendedDescriptor_128();  
		}
      }
	  // G-SURF descriptor (100,200)
	  else if( mode == 5 )
      {
		if( extended == false )
		{
			getGaugeDescriptor_100();  
		}
		else
		{
			getGaugeUprightExtendedDescriptor_200();  
		}
      }
	  // G-SURF descriptor (144,288)
	  else if( mode == 6 )
      {
		if( extended == false )
		{
			getGaugeDescriptor_144();  
		}
		else
		{
			getGaugeExtendedDescriptor_288();  
		}
      }
	  // SURF no Gaussian Weighting descriptor (16,32)
	  else if( mode == 7 )
      {
		if( extended == false )
		{
			getDescriptor_NW_16();  
		}
		else
		{
			getExtendedDescriptor_NW_32();   
		}
      }
	  // SURF no Gaussian Weighting descriptor (36,72)
	  else if( mode == 8 )
      {
		if( extended == false )
		{
			getDescriptor_NW_36();  
		}
		else
		{
			getExtendedDescriptor_NW_72();   
		}
      }	  
	  // SURF no Gaussian Weighting descriptor (64,128)
	  else if( mode == 9 )
      {
		if( extended == false )
		{
			getDescriptor_NW_64();  
		}
		else
		{
			getExtendedDescriptor_NW_128();   
		}
      }	  
	  // SURF no Gaussian Weighting descriptor (100,200)
	  else if( mode == 10 )
      {
		if( extended == false )
		{
			getDescriptor_NW_100();  
		}
		else
		{
			getExtendedDescriptor_NW_200();   
		}
      }
	  // SURF no Gaussian Weighting descriptor (144,288)
	  else if( mode == 11 )
      {
		if( extended == false )
		{
			getDescriptor_NW_144();  
		}
		else
		{
			getExtendedDescriptor_NW_288();   
		}
      }
	  else if( mode == 12 )
	  {
		  if( extended == false )
		  {
			  getMGaugeDescriptor_64();
		  }
		  else
		  {
			  getMGaugeExtendedDescriptor_128();
		  }
	  }
    }
  }
  
  gettimeofday(&tod1,NULL);
  t2 = tod1.tv_sec + tod1.tv_usec / 1000000.0;
  t2 = 1000.0*(t2 - t1);
  
}

//**************************************************************************************
//**************************************************************************************

//! Assign the supplied Ipoint an orientation
void Surf::getOrientation()
{
  Ipoint *ipt = &ipts[index];
  float gauss = 0.f, scale = ipt->scale;
  const int s = fRound(scale), r = fRound(ipt->y), c = fRound(ipt->x);
  std::vector<float> resX(109), resY(109), Ang(109);
  const int id[] = {6,5,4,3,2,1,0,1,2,3,4,5,6};
  int idx = 0;
  
  // calculate haar responses for points within radius of 6*scale
  for(int i = -6; i <= 6; ++i) 
  {
    for(int j = -6; j <= 6; ++j) 
    {
      if(i*i + j*j < 36) 
      {
        gauss = static_cast<float>(gauss25[id[i+6]][id[j+6]]);
        resX[idx] = gauss * haarX(r+j*s, c+i*s, 4*s);
        resY[idx] = gauss * haarY(r+j*s, c+i*s, 4*s);
        Ang[idx] = getAngle(resX[idx], resY[idx]);
        ++idx;
      }
    }
  }

  // calculate the dominant direction 
  float sumX=0.f, sumY=0.f;
  float max=0.f, orientation = 0.f;
  float ang1=0.f, ang2=0.f;

  // loop slides pi/3 window around feature point
  for(ang1 = 0; ang1 < 2*pi;  ang1+=0.15f)
  {
    ang2 = ( ang1+pi/3.0f > 2*pi ? ang1-5.0f*pi/3.0f : ang1+pi/3.0f);
    sumX = sumY = 0.f; 
	
    for(unsigned int k = 0; k < Ang.size(); ++k) 
    {
      // get angle from the x-axis of the sample point
      const float & ang = Ang[k];

      // determine whether the point is within the window
      if (ang1 < ang2 && ang1 < ang && ang < ang2) 
      {
        sumX+=resX[k];  
        sumY+=resY[k];
      } 
      else if (ang2 < ang1 && 
        ((ang > 0 && ang < ang2) || (ang > ang1 && ang < 2*pi) )) 
      {
        sumX+=resX[k];  
        sumY+=resY[k];
      }
    }

    // if the vector produced from this window is longer than all 
    // previous vectors then this forms the new dominant direction
    if (sumX*sumX + sumY*sumY > max) 
    {
      // store largest orientation
      max = sumX*sumX + sumY*sumY;
      orientation = getAngle(sumX, sumY);
    }
  }

  // assign orientation of the dominant response vector
  ipt->orientation = orientation;
}

//**************************************************************************************
//**************************************************************************************

/** Get the descriptor vector of the provided Ipoint								   */
/** Original OpenSURF descriptor with Gaussian Weighting as proposed in Bay et al. 06  */
/** Descriptor Dimension 64															   */
void Surf::getDescriptor_64()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0, co = 0.0, si = 0.0;
  float gauss = 0.0, rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len=0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  dsize = ipt->size;
  desc = ipt->descriptor;
  sample_step = 5;
  pattern_size = 10;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dx=dy=mdx=mdy=0.0;
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          // Get the gaussian weighted x and y responses
          gauss = static_cast<float>(gauss33[abs(k)][abs(l)]);
		  
          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          // Get the gaussian weighted x and y responses on rotated axis
          rrx = -rx*si + ry*co;
          rry = rx*co + ry*si;

          rrx = gauss*rrx;
          rry = gauss*rry;

          dx += rrx;
          dy += rry;
          mdx += fabs(rrx);
          mdy += fabs(rry);
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dx;
      desc[count++] = dy;
      desc[count++] = mdx;
      desc[count++] = mdy;

      // store the current length^2 of the vector
      len += dx*dx + dy*dy + mdx*mdx + mdy*mdy;
    }
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the upright descriptor vector of the provided Ipoint						   */
/** Original OpenSURF descriptor with Gaussian Weighting as proposed in Bay et al. 06  */
/** Descriptor Dimension 64															   */
void Surf::getUprightDescriptor_64()
{
  int y = 0, x = 0, dsize = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0;
  float gauss = 0.0, rx = 0.0, ry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  y = fRound(ipt->y);
  x = fRound(ipt->x);
  dsize = ipt->size;
  desc = ipt->descriptor;  
  sample_step = 5;
  pattern_size = 10;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dx=dy=mdx=mdy=0;
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // get Gaussian weighted x and y responses
          gauss = static_cast<float>(gauss33[abs(k)][abs(l)]);
		  
          rx = gauss*haarX(fRound(k*scale+y), fRound(l*scale+x), 2*fRound(scale));
          ry = gauss*haarY(fRound(k*scale+y), fRound(l*scale+x), 2*fRound(scale));

          dx += rx;
          dy += ry;
          mdx += fabs(rx);
          mdy += fabs(ry);
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dx;
      desc[count++] = dy;
      desc[count++] = mdx;
      desc[count++] = mdy;

      // store the current length^2 of the vector
      len += dx*dx + dy*dy + mdx*mdx + mdy*mdy;
    }
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the 128 descriptor vector of the provided Ipoint						       */
/** Original OpenSURF descriptor with Gaussian Weighting as proposed in Bay et al. 06  */
/** Descriptor Dimension 128														   */
void Surf::getExtendedDescriptor_128()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float co = 0.0, si = 0.0, dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float gauss = 0.0, rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  sample_step = 5;
  pattern_size = 10;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dxp=dyp=mdxp=mdyp=0.0f;
      dxn=dyn=mdxn=mdyn=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          // Get the gaussian weighted x and y responses
          gauss = static_cast<float>(gauss33[abs(k)][abs(l)]);

          rx = gauss*haarX(sample_y, sample_x, 2*fRound(scale));
          ry = gauss*haarY(sample_y, sample_x, 2*fRound(scale));

          // Get the gaussian weighted x and y responses on rotated axis
          rrx = -rx*si + ry*co;
          rry = rx*co + ry*si;

          if( rry >= 0.0 )
          {
              dxp += rrx;
              mdxp += fabs(rrx);
          }
          else
          {
              dxn += rrx;
              mdxn += fabs(rrx);
          }

          if( rrx >= 0.0 )
          {
              dyp += rry;
              mdyp += fabs(rry);
          }
          else
          {
              dyn += rry;
              mdyn += fabs(rry);
          }
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dxp;
      desc[count++] = mdxp;
      desc[count++] = dxn;
      desc[count++] = mdxn;

      desc[count++] = dyp;
      desc[count++] = mdyp;
      desc[count++] = dyn;
      desc[count++] = mdyn;
      }

      // store the current length^2 of the vector
      len += dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn;

  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the 128 upright descriptor vector of the provided Ipoint					   */
/** Original OpenSURF descriptor with Gaussian Weighting as proposed in Bay et al. 06  */
/** Descriptor Dimension 128														   */
void Surf::getUprightExtendedDescriptor_128()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float gauss = 0.0, rx = 0.0, ry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  desc = ipt->descriptor;
  dsize = ipt->size;
  sample_step = 5;
  pattern_size = 10;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dxp=dyp=mdxp=mdyp=0.0;
      dxn=dyn=mdxn=mdyn=0.0;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + k*scale);
          sample_y = fRound(y + l*scale);

          // Get the gaussian weighted x and y responses
          gauss = static_cast<float>(gauss33[abs(k)][abs(l)]);

          rx = gauss*haarX(sample_y, sample_x, 2*fRound(scale));
          ry = gauss*haarY(sample_y, sample_x, 2*fRound(scale));

          if( ry >= 0.0 )
          {
              dxp += rx;
              mdxp += fabs(rx);
          }
          else
          {
              dxn += rx;
              mdxn += fabs(rx);
          }

          if( rx >= 0.0 )
          {
              dyp += ry;
              mdyp += fabs(ry);
          }
          else
          {
              dyn += ry;
              mdyn += fabs(ry);
          }
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dxp;
      desc[count++] = mdxp;
      desc[count++] = dxn;
      desc[count++] = mdxn;

      desc[count++] = dyp;
      desc[count++] = mdyp;
      desc[count++] = dyn;
      desc[count++] = mdyn;
      }

      // store the current length^2 of the vector
      len += dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn;
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the modified SURF descriptor vector of the provided Ipoint					   */
/** See Agrawal et al. (ECCV08) for more details  									   */
/** Descriptor Dimension 64															   */
void Surf::getMDescriptor_64()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count=0;
  int sample_step = 0, pattern_size = 0, i = 0, ix = 0, j = 0, jx = 0, xs = 0, ys = 0;
  float *desc, scale = 0.0, dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0; 
  float gauss_s1 = 0.0, gauss_s2 = 0.0, co = 0.0, si = 0.0;
  float rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len = 0.0;
  float cx = -0.5, cy = 0.5; //Subregion centers for the 4x4 gaussian weighting

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size ;
  sample_step = 5;
  pattern_size = 12;
  
  i = -8;

  // Calculate descriptor for this interest point
  // Area of size 24 s x 24 s
  while(i < pattern_size)
  {
    j = -8;
    i = i-4;

    cx += 1.f;
    cy = -0.5f;

    while(j < 12)
    {
      dx=dy=mdx=mdy=0.f;
      cy += 1.f;

      j = j - 4;

      ix = i + sample_step;
      jx = j + sample_step;

      xs = fRound(x + ( -jx*scale*si + ix*scale*co));
      ys = fRound(y + ( jx*scale*co + ix*scale*si));

      for (int k = i; k < i + 9; ++k)
      {
        for (int l = j; l < j + 9; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          // Get the gaussian weighted x and y responses
          gauss_s1 = gaussian(xs-sample_x,ys-sample_y,2.5f*scale);

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          // Get the gaussian weighted x and y responses on rotated axis
          rrx = -rx*si + ry*co;
          rry = rx*co + ry*si;

          rrx = gauss_s1*rrx;
          rry = gauss_s1*rry;

          dx += rrx;
          dy += rry;
          mdx += fabs(rrx);
          mdy += fabs(rry);

        }
      }

      // Add the values to the descriptor vector
      gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f);

      desc[count++] = dx*gauss_s2;
      desc[count++] = dy*gauss_s2;
      desc[count++] = mdx*gauss_s2;
      desc[count++] = mdy*gauss_s2;

      len += (dx*dx + dy*dy + mdx*mdx + mdy*mdy)*gauss_s2*gauss_s2;

      j += 9;
    }
    i += 9;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the modified upright SURF descriptor vector of the provided Ipoint			   */
/** See Agrawal et al. (ECCV08) for more details  									   */
/** Descriptor Dimension 64															   */
void Surf::getMUprightDescriptor_64()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count=0;
  int sample_step = 0, pattern_size = 0, i = 0, ix = 0, j = 0, jx = 0, xs = 0, ys = 0;
  float *desc, scale = 0.0, dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0; 
  float gauss_s1 = 0.0, gauss_s2 = 0.0;
  float rx = 0.0, ry = 0.0, len = 0.0;
  float cx = -0.5, cy = 0.5; // Subregion centers for the 4x4 gaussian weighting

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  desc = ipt->descriptor;
  dsize = ipt->size;
  sample_step = 5;
  pattern_size = 12;
  
  i = -8;

  // Calculate descriptor for this interest point
  // Area of size 24 s x 24 s
  while(i < pattern_size)
  {
    j = -8;
    i = i-4;

    cx += 1.0f;
    cy = -0.5f;

    while(j < pattern_size)
    {
      dx=dy=mdx=mdy=0.0f;
      cy += 1.0f;

      j = j - 4;

      ix = i + sample_step;
      jx = j + sample_step;

      xs = fRound(x + (ix*scale));
      ys = fRound(y + (jx*scale));

      for (int k = i; k < i + 9; ++k)
      {
        for (int l = j; l < j + 9; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + k*scale);
          sample_y = fRound(y + l*scale);

          // Get the gaussian weighted x and y responses
          gauss_s1 = gaussian(xs-sample_x,ys-sample_y,2.5f*scale);

          rx = gauss_s1*haarX(sample_y, sample_x, 2*fRound(scale));
          ry = gauss_s1*haarY(sample_y, sample_x, 2*fRound(scale));

          dx += rx;
          dy += ry;
          mdx += fabs(rx);
          mdy += fabs(ry);
        }
      }
	  
      // Add the values to the descriptor vector
      gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f);

      desc[count++] = dx*gauss_s2;
      desc[count++] = dy*gauss_s2;
      desc[count++] = mdx*gauss_s2;
      desc[count++] = mdy*gauss_s2;

      len += (dx*dx + dy*dy + mdx*mdx + mdy*mdy)*gauss_s2*gauss_s2;

      j += 9;
    }

      i += 9;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the modified 128 SURF descriptor vector of the provided Ipoint			   	   */
/** See Agrawal et al. (ECCV08) for more details  									   */
/** Descriptor Dimension 128														   */
void Surf::getMExtendedDescriptor_128()
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count=0, dsize = 0;
  int i = 0, ix = 0, j = 0, jx = 0, xs = 0, ys = 0;
  float sample_step = 0, pattern_size = 0, scale = 0.0, *desc, co = 0.0, si = 0.0;
  float dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float gauss_s1 = 0.f, gauss_s2 = 0.f;
  float rx = 0.f, ry = 0.f, rrx = 0.f, rry = 0.f, len = 0.f;
  float cx = -0.5f, cy = 0.f; //Subregion centers for the 4x4 gaussian weighting

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  sample_step = 5;
  pattern_size = 12;
  
  i = -8;

  // Calculate descriptor for this interest point
  // Area of size 24 s x 24 s
  while(i < pattern_size)
  {
    j = -8;
    i = i-4;

    cx += 1.f;
    cy = -0.5f;

    while(j < pattern_size)
    {
      dxp=dyp=mdxp=mdyp=0.0f;
      dxn=dyn=mdxn=mdyn=0.0f;

      cy += 1.f;

      j = j - 4;

      ix = i + sample_step;
      jx = j + sample_step;

      xs = fRound(x + ( -jx*scale*si + ix*scale*co));
      ys = fRound(y + ( jx*scale*co + ix*scale*si));

      for (int k = i; k < i + 9; ++k)
      {
        for (int l = j; l < j + 9; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          // Get the gaussian weighted x and y responses
          gauss_s1 = gaussian(xs-sample_x,ys-sample_y,2.5f*scale);

          rx = gauss_s1*haarX(sample_y, sample_x, 2*fRound(scale));
          ry = gauss_s1*haarY(sample_y, sample_x, 2*fRound(scale));

          // Get the gaussian weighted x and y responses on rotated axis
          rrx = -rx*si + ry*co;
          rry = rx*co + ry*si;

          if( rry >= 0.0 )
          {
              dxp += rrx;
              mdxp += fabs(rrx);
          }
          else
          {
              dxn += rrx;
              mdxn += fabs(rrx);
          }

          if( rrx >= 0.0 )
          {
              dyp += rry;
              mdyp += fabs(rry);
          }
          else
          {
              dyn += rry;
              mdyn += fabs(rry);
          }
        }
      }

      // Add the values to the descriptor vector
      gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f);

      desc[count++] = dxp*gauss_s2;
      desc[count++] = mdxp*gauss_s2;
      desc[count++] = dxn*gauss_s2;
      desc[count++] = mdxn*gauss_s2;

      desc[count++] = dyp*gauss_s2;
      desc[count++] = mdyp*gauss_s2;
      desc[count++] = dyn*gauss_s2;
      desc[count++] = mdyn*gauss_s2;

      // store the current length^2 of the vector
      len += (dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn)*gauss_s2*gauss_s2;

      j += 9;
    }
    i += 9;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the modified 128 upright SURF descriptor vector of the provided Ipoint		   */
/** See Agrawal et al. (ECCV08) for more details  									   */
/** Descriptor Dimension 128														   */
void Surf::getMUprightExtendedDescriptor_128()
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count=0, dsize = 0;
  int i = 0, ix = 0, j = 0, jx = 0, xs = 0, ys = 0, sample_step = 0, pattern_size = 0;
  float scale = 0.0, *desc;
  float dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float gauss_s1 = 0.f, gauss_s2 = 0.f;
  float rx = 0.f, ry = 0.f, len = 0.f;
  float cx = -0.5f, cy = 0.f; // Subregion centers for the 4x4 gaussian weighting

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  desc = ipt->descriptor;
  dsize = ipt->size;
  sample_step = 5;
  pattern_size = 12;
  
  i = -8;

  // Calculate descriptor for this interest point
  // Area of size 24 s x 24 s
  while(i < pattern_size)
  {
    j = -8;
    i = i-4;

    cx += 1.f;
    cy = -0.5f;

    while(j < pattern_size)
    {
      dxp=dyp=mdxp=mdyp=0.0f;
      dxn=dyn=mdxn=mdyn=0.0f;

      cy += 1.f;

      j = j - 4;

      ix = i + sample_step;
      jx = j + sample_step;

      xs = fRound(x + ix*scale);
      ys = fRound(y + jx*scale);

      for (int k = i; k < i + 9; ++k)
      {
        for (int l = j; l < j + 9; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + k*scale);
          sample_y = fRound(y + l*scale);

          // Get the gaussian weighted x and y responses
          gauss_s1 = gaussian(xs-sample_x,ys-sample_y,2.5f*scale);

          rx = gauss_s1*haarX(sample_y, sample_x, 2*fRound(scale));
          ry = gauss_s1*haarY(sample_y, sample_x, 2*fRound(scale));

          if( ry >= 0.0 )
          {
              dxp += rx;
              mdxp += fabs(rx);
          }
          else
          {
              dxn += rx;
              mdxn += fabs(rx);
          }

          if( rx >= 0.0 )
          {
              dyp += ry;
              mdyp += fabs(ry);
          }
          else
          {
              dyn += ry;
              mdyn += fabs(ry);
          }
        }
      }

      // Add the values to the descriptor vector
      gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f);

      desc[count++] = dxp*gauss_s2;
      desc[count++] = mdxp*gauss_s2;
      desc[count++] = dxn*gauss_s2;
      desc[count++] = mdxn*gauss_s2;

      desc[count++] = dyp*gauss_s2;
      desc[count++] = mdyp*gauss_s2;
      desc[count++] = dyn*gauss_s2;
      desc[count++] = mdyn*gauss_s2;

      // store the current length^2 of the vector
      len += (dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn)*gauss_s2*gauss_s2;

      j += 9;
    }
    i += 9;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;
	
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the Gauge SURF descriptor vector of the provided Ipoint		   				  */
/** Descriptor Dimension 64															  */
void Surf::getGaugeDescriptor_64()
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, sample_step = 0, pattern_size = 0;
  float scale = 0.0, *desc, co = 0.0, si = 0.0;
  float Lvv = 0.0, Lww = 0.0, mLww = 0.0, mLvv = 0.0, lvv = 0.0, lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0, modg = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 10;
  sample_step = 5;

  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 20 s x 20 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvv=Lww=mLvv=mLww=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

		  rx = haarX(sample_y,sample_x,2*fRound(scale));
          ry = haarY(sample_y,sample_x,2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);

          if( modg != 0.0 )
          {
               rxx = haarXX(sample_y,sample_x,2*fRound(scale));
               rxy = haarXY(sample_y,sample_x,2*fRound(scale));
               ryy = haarYY(sample_y,sample_x,2*fRound(scale));

			   // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
			   lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

			   // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
			   lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

                Lww += lww;
                Lvv += lvv;
                mLww += fabs(lww);
                mLvv += fabs(lvv);
        }

      }

      // Add the values to the descriptor vector
      desc[count++] = Lww;
      desc[count++] = Lvv;
      desc[count++] = mLww;
      desc[count++] = mLvv;

      len += Lww*Lww + Lvv*Lvv + mLww*mLww + mLvv*mLvv;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; ++i)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the Gauge SURF upright descriptor vector of the provided Ipoint		 		  */
/** Descriptor Dimension 64															  */
void Surf::getGaugeUprightDescriptor_64()
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, sample_step = 0, pattern_size = 0;
  float scale = 0.0, *desc;
  float Lvv = 0.0, Lww = 0.0, mLww = 0.0, mLvv = 0.0, lvv = 0.0, lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0, modg = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 10;
  sample_step = 5;
  
  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 20 s x 20 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvv=Lww=mLvv=mLww=0.0f;
	  
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + k*scale);
          sample_y = fRound(y + l*scale);
		  
          rx = haarX(sample_y,sample_x,2.0*fRound(scale));
          ry = haarY(sample_y,sample_x,2.0*fRound(scale));
		  
          modg = pow(rx,2) + pow(ry,2);
		  
          if( modg != 0.0 )
          {    
			   rxx =  haarXX(sample_y,sample_x,2*fRound(scale));
			   rxy =  haarXY(sample_y,sample_x,2*fRound(scale));
			   ryy =  haarYY(sample_y,sample_x,2*fRound(scale));

			   // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
			   lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

			   // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
			   lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }
	
           Lww += lww;
           Lvv += lvv;
		   mLww += 1.0*fabs(lww);
           mLvv += 1.0*fabs(lvv);
        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lww;
      desc[count++] = Lvv;
      desc[count++] = mLww;
      desc[count++] = mLvv;
	  
      len += Lww*Lww + Lvv*Lvv + mLww*mLww + mLvv*mLvv;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; ++i)
    desc[i] /= len;
	
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the Gauge SURF extended descriptor vector of the provided Ipoint		  	   */
/** Descriptor Dimension 128														   */
void Surf::getGaugeExtendedDescriptor_128(void)
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, pattern_size = 0, sample_step = 0;
  float scale = 0.0, *desc, co = 0.0, si = 0.0;
  float Lvvp = 0.0, Lwwp = 0.0, mLwwp = 0.0, mLvvp = 0.0;
  float Lvvn = 0.0, Lwwn = 0.0, mLwwn = 0.0, mLvvn = 0.0;
  float lvv = 0.0, lww = 0.0, modg = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  sample_step = 5;
  pattern_size = 10;
  
  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 20 s x 20 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvvp=Lwwp=mLvvp=mLwwp=0.0f;
      Lvvn=Lwwn=mLvvn=mLwwn=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y,sample_x,2*fRound(scale));
          ry = haarY(sample_y,sample_x,2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);

          if( modg != 0.0 )
          {    
                rxx = haarXX(sample_y,sample_x,2*fRound(scale));
                rxy = haarXY(sample_y,sample_x,2*fRound(scale));
                ryy = haarYY(sample_y,sample_x,2*fRound(scale));

                // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
                lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

                // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
                lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

          if( lww >= 0.0 )
          {
              Lvvp += lvv;
              mLvvp += fabs(lvv);
          }
          else
          {
              Lvvn += lvv;
              mLvvn += fabs(lvv);
          }

          if( lvv >= 0.0 )
          {
              Lwwp += lww;
              mLwwp += fabs(lww);
          }
          else
          {
              Lwwn += lww;
              mLwwn += fabs(lww);
          }

        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lwwn;
      desc[count++] = Lvvn;
      desc[count++] = mLwwn;
      desc[count++] = mLvvn;

      desc[count++] = Lwwp;
      desc[count++] = Lvvp;
      desc[count++] = mLwwp;
      desc[count++] = mLvvp;

      len += Lwwp*Lwwp + Lvvp*Lvvp + mLwwp*mLwwp + mLvvp*mLvvp + Lwwn*Lwwn + 
			 Lvvn*Lvvn + mLwwn*mLwwn + mLvvn*mLvvn;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;
		
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//************************************************************************
//************************************************************************

/** Get the Gauge SURF upright extended descriptor vector of the provided Ipoint	   */
/** Descriptor Dimension 128														   */
void Surf::getGaugeUprightExtendedDescriptor_128(void)
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, sample_step = 0, pattern_size = 0;
  float scale = 0.0, *desc, co = 0.0, si = 0.0;
  float Lvvp = 0.0, Lwwp = 0.0, mLwwp = 0.0, mLvvp = 0.0;
  float Lvvn = 0.0, Lwwn = 0.0, mLwwn = 0.0, mLvvn = 0.0;
  float lvv = 0.0, lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0, modg = 0.0;
  
  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  sample_step = 5;
  pattern_size = 10;
  
  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 20 s x 20 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvvp=Lwwp=mLvvp=mLwwp=0.0f;
      Lvvn=Lwwn=mLvvn=mLwwn=0.0f;
	  
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y,sample_x,2*fRound(scale));
          ry = haarY(sample_y,sample_x,2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);

          if( modg != 0.0 )
          {
                rxx = haarXX(sample_y,sample_x,2*fRound(scale));
                rxy = haarXY(sample_y,sample_x,2*fRound(scale));
                ryy = haarYY(sample_y,sample_x,2*fRound(scale));

                // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
                lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

                // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
                lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

          if( lww >= 0.0 )
          {
              Lvvp += lvv;
              mLvvp += fabs(lvv);
          }
          else
          {
              Lvvn += lvv;
              mLvvn += fabs(lvv);
          }

          if( lvv >= 0.0 )
          {
              Lwwp += lww;
              mLwwp += fabs(lww);
          }
          else
          {
              Lwwn += lww;
              mLwwn += fabs(lww);
          }

        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lwwn;
      desc[count++] = Lvvn;
      desc[count++] = mLwwn;
      desc[count++] = mLvvn;

      desc[count++] = Lwwp;
      desc[count++] = Lvvp;
      desc[count++] = mLwwp;
      desc[count++] = mLvvp;

      len += Lwwp*Lwwp + Lvvp*Lvvp + mLwwp*mLwwp + mLvvp*mLvvp + Lwwn*Lwwn + 
			 Lvvn*Lvvn + mLwwn*mLwwn + mLvvn*mLvvn;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;
  
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the Gauge SURF descriptor vector of the provided Ipoint		   				  */
/** Descriptor Dimension 16															  */
void Surf::getGaugeDescriptor_16()
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, sample_step = 0, pattern_size = 0;
  float scale = 0.0, *desc, co = 0.0, si = 0.0;
  float Lvv = 0.0, Lww = 0.0, mLww = 0.0, mLvv = 0.0, lvv = 0.0, lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0, modg = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 5;
  sample_step = 5;

  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 10 s x 10 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvv=Lww=mLvv=mLww=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

		  rx = haarX(sample_y,sample_x,2*fRound(scale));
          ry = haarY(sample_y,sample_x,2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);

          if( modg != 0.0 )
          {
               rxx = haarXX(sample_y,sample_x,2*fRound(scale));
               rxy = haarXY(sample_y,sample_x,2*fRound(scale));
               ryy = haarYY(sample_y,sample_x,2*fRound(scale));

			   // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
			   lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

			   // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
			   lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

                Lww += lww;
                Lvv += lvv;
                mLww += fabs(lww);
                mLvv += fabs(lvv);
        }

      }

      // Add the values to the descriptor vector
      desc[count++] = Lww;
      desc[count++] = Lvv;
      desc[count++] = mLww;
      desc[count++] = mLvv;

      len += Lww*Lww + Lvv*Lvv + mLww*mLww + mLvv*mLvv;

      j += sample_step;
    }

      i += sample_step;
  }

  //Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; ++i)
    desc[i] /= len;
	
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the Gauge SURF upright descriptor vector of the provided Ipoint		 		  */
/** Descriptor Dimension 16															  */
void Surf::getGaugeUprightDescriptor_16()
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, sample_step = 0, pattern_size = 0;
  float scale = 0.0, *desc;
  float Lvv = 0.0, Lww = 0.0, mLww = 0.0, mLvv = 0.0, lvv = 0.0, lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0, modg = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 5;
  sample_step = 5;
  
  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 20 s x 20 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvv=Lww=mLvv=mLww=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + k*scale);
          sample_y = fRound(y + l*scale);

          rx = haarX(sample_y,sample_x,2.0*fRound(scale));
          ry = haarY(sample_y,sample_x,2.0*fRound(scale));
		  
          modg = pow(rx,2) + pow(ry,2);
		  
          if( modg != 0.0 )
          {    
			   rxx =  haarXX(sample_y,sample_x,2.0*fRound(scale));
			   rxy =  haarXY(sample_y,sample_x,2.0*fRound(scale));
			   ryy =  haarYY(sample_y,sample_x,2.0*fRound(scale));

			   // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
			   lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

			   // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
			   lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }
	
           Lww += lww;
           Lvv += lvv;
		   mLww += 1.0*fabs(lww);
           mLvv += 1.0*fabs(lvv);
        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lww;
      desc[count++] = Lvv;
      desc[count++] = mLww;
      desc[count++] = mLvv;
	  
      len += Lww*Lww + Lvv*Lvv + mLww*mLww + mLvv*mLvv;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; ++i)
    desc[i] /= len;
	
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the Gauge SURF extended descriptor vector of the provided Ipoint		  	   */
/** Descriptor Dimension 32														   	   */
void Surf::getGaugeExtendedDescriptor_32(void)
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, pattern_size = 0, sample_step = 0;
  float scale = 0.0, *desc, co = 0.0, si = 0.0;
  float Lvvp = 0.0, Lwwp = 0.0, mLwwp = 0.0, mLvvp = 0.0;
  float Lvvn = 0.0, Lwwn = 0.0, mLwwn = 0.0, mLvvn = 0.0;
  float lvv = 0.0, lww = 0.0, modg = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  sample_step = 10;
  pattern_size = 10;
  
  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 20 s x 20 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvvp=Lwwp=mLvvp=mLwwp=0.0f;
      Lvvn=Lwwn=mLvvn=mLwwn=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y,sample_x,2*fRound(scale));
          ry = haarY(sample_y,sample_x,2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);

          if( modg != 0.0 )
          {    
                rxx = haarXX(sample_y,sample_x,2*fRound(scale));
                rxy = haarXY(sample_y,sample_x,2*fRound(scale));
                ryy = haarYY(sample_y,sample_x,2*fRound(scale));

                // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
                lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

                // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
                lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

          if( lww >= 0.0 )
          {
              Lvvp += lvv;
              mLvvp += fabs(lvv);
          }
          else
          {
              Lvvn += lvv;
              mLvvn += fabs(lvv);
          }

          if( lvv >= 0.0 )
          {
              Lwwp += lww;
              mLwwp += fabs(lww);
          }
          else
          {
              Lwwn += lww;
              mLwwn += fabs(lww);
          }

        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lwwn;
      desc[count++] = Lvvn;
      desc[count++] = mLwwn;
      desc[count++] = mLvvn;

      desc[count++] = Lwwp;
      desc[count++] = Lvvp;
      desc[count++] = mLwwp;
      desc[count++] = mLvvp;

      len += Lwwp*Lwwp + Lvvp*Lvvp + mLwwp*mLwwp + mLvvp*mLvvp + Lwwn*Lwwn + 
			 Lvvn*Lvvn + mLwwn*mLwwn + mLvvn*mLvvn;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;
 
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//************************************************************************
//************************************************************************

/** Get the Gauge SURF upright extended descriptor vector of the provided Ipoint	   */
/** Descriptor Dimension 32														   	   */
void Surf::getGaugeUprightExtendedDescriptor_32(void)
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, sample_step = 0, pattern_size = 0;
  float scale = 0.0, *desc, co = 0.0, si = 0.0;
  float Lvvp = 0.0, Lwwp = 0.0, mLwwp = 0.0, mLvvp = 0.0;
  float Lvvn = 0.0, Lwwn = 0.0, mLwwn = 0.0, mLvvn = 0.0;
  float lvv = 0.0, lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0, modg = 0.0;
  
  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  sample_step = 10;
  pattern_size = 10;
  
  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 20 s x 20 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvvp=Lwwp=mLvvp=mLwwp=0.0f;
      Lvvn=Lwwn=mLvvn=mLwwn=0.0f;
	  
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y,sample_x,2*fRound(scale));
          ry = haarY(sample_y,sample_x,2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);

          if( modg != 0.0 )
          {
                rxx = haarXX(sample_y,sample_x,2*fRound(scale));
                rxy = haarXY(sample_y,sample_x,2*fRound(scale));
                ryy = haarYY(sample_y,sample_x,2*fRound(scale));

                // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
                lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

                // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
                lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

          if( lww >= 0.0 )
          {
              Lvvp += lvv;
              mLvvp += fabs(lvv);
          }
          else
          {
              Lvvn += lvv;
              mLvvn += fabs(lvv);
          }

          if( lvv >= 0.0 )
          {
              Lwwp += lww;
              mLwwp += fabs(lww);
          }
          else
          {
              Lwwn += lww;
              mLwwn += fabs(lww);
          }

        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lwwn;
      desc[count++] = Lvvn;
      desc[count++] = mLwwn;
      desc[count++] = mLvvn;

      desc[count++] = Lwwp;
      desc[count++] = Lvvp;
      desc[count++] = mLwwp;
      desc[count++] = mLvvp;

      len += Lwwp*Lwwp + Lvvp*Lvvp + mLwwp*mLwwp + mLvvp*mLvvp + Lwwn*Lwwn + 
			 Lvvn*Lvvn + mLwwn*mLwwn + mLvvn*mLvvn;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}
//**************************************************************************************
//**************************************************************************************

/** Get the Gauge SURF descriptor vector of the provided Ipoint		   				  */
/** Descriptor Dimension 144														  */
void Surf::getGaugeDescriptor_144()
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, sample_step = 0, pattern_size = 0;
  float scale = 0.0, *desc, co = 0.0, si = 0.0;
  float Lvv = 0.0, Lww = 0.0, mLww = 0.0, mLvv = 0.0, lvv = 0.0, lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0, modg = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 12;
  sample_step = 4;

  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 24 s x 24 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvv=Lww=mLvv=mLww=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

		  rx = haarX(sample_y,sample_x,2*fRound(scale));
          ry = haarY(sample_y,sample_x,2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);

          if( modg != 0.0 )
          {
               rxx = haarXX(sample_y,sample_x,2*fRound(scale));
               rxy = haarXY(sample_y,sample_x,2*fRound(scale));
               ryy = haarYY(sample_y,sample_x,2*fRound(scale));

			   // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
			   lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

			   // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
			   lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

                Lww += lww;
                Lvv += lvv;
                mLww += fabs(lww);
                mLvv += fabs(lvv);
        }

      }

      // Add the values to the descriptor vector
      desc[count++] = Lww;
      desc[count++] = Lvv;
      desc[count++] = mLww;
      desc[count++] = mLvv;

      len += Lww*Lww + Lvv*Lvv + mLww*mLww + mLvv*mLvv;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; ++i)
    desc[i] /= len;
	
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the Gauge SURF upright descriptor vector of the provided Ipoint		 		  */
/** Descriptor Dimension 144														  */
void Surf::getGaugeUprightDescriptor_144()
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, sample_step = 0, pattern_size = 0;
  float scale = 0.0, *desc;
  float Lvv = 0.0, Lww = 0.0, mLww = 0.0, mLvv = 0.0, lvv = 0.0, lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0, modg = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 12;
  sample_step = 4;
  
  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 24 s x 24 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvv=Lww=mLvv=mLww=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + k*scale);
          sample_y = fRound(y + l*scale);

          rx = haarX(sample_y,sample_x,2.0*fRound(scale));
          ry = haarY(sample_y,sample_x,2.0*fRound(scale));
		  
          modg = pow(rx,2) + pow(ry,2);
		  
          if( modg != 0.0 )
          {    
			   rxx =  haarXX(sample_y,sample_x,2*fRound(scale));
			   rxy =  haarXY(sample_y,sample_x,2*fRound(scale));
			   ryy =  haarYY(sample_y,sample_x,2*fRound(scale));

			   // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
			   lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

			   // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
			   lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }
	
           Lww += lww;
           Lvv += lvv;
		   mLww += 1.0*fabs(lww);
           mLvv += 1.0*fabs(lvv);
        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lww;
      desc[count++] = Lvv;
      desc[count++] = mLww;
      desc[count++] = mLvv;
	  
      len += Lww*Lww + Lvv*Lvv + mLww*mLww + mLvv*mLvv;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; ++i)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the Gauge SURF extended descriptor vector of the provided Ipoint		  	   */
/** Descriptor Dimension 288														   */
void Surf::getGaugeExtendedDescriptor_288(void)
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, pattern_size = 0, sample_step = 0;
  float scale = 0.0, *desc, co = 0.0, si = 0.0;
  float Lvvp = 0.0, Lwwp = 0.0, mLwwp = 0.0, mLvvp = 0.0;
  float Lvvn = 0.0, Lwwn = 0.0, mLwwn = 0.0, mLvvn = 0.0;
  float lvv = 0.0, lww = 0.0, modg = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 12;
  sample_step = 4;
  
  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 24 s x 24 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvvp=Lwwp=mLvvp=mLwwp=0.0f;
      Lvvn=Lwwn=mLvvn=mLwwn=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y,sample_x,2*fRound(scale));
          ry = haarY(sample_y,sample_x,2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);

          if( modg != 0.0 )
          {    
                rxx = haarXX(sample_y,sample_x,2*fRound(scale));
                rxy = haarXY(sample_y,sample_x,2*fRound(scale));
                ryy = haarYY(sample_y,sample_x,2*fRound(scale));

                // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
                lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

                // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
                lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

          if( lww >= 0.0 )
          {
              Lvvp += lvv;
              mLvvp += fabs(lvv);
          }
          else
          {
              Lvvn += lvv;
              mLvvn += fabs(lvv);
          }

          if( lvv >= 0.0 )
          {
              Lwwp += lww;
              mLwwp += fabs(lww);
          }
          else
          {
              Lwwn += lww;
              mLwwn += fabs(lww);
          }

        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lwwn;
      desc[count++] = Lvvn;
      desc[count++] = mLwwn;
      desc[count++] = mLvvn;

      desc[count++] = Lwwp;
      desc[count++] = Lvvp;
      desc[count++] = mLwwp;
      desc[count++] = mLvvp;

      len += Lwwp*Lwwp + Lvvp*Lvvp + mLwwp*mLwwp + mLvvp*mLvvp + Lwwn*Lwwn + 
			 Lvvn*Lvvn + mLwwn*mLwwn + mLvvn*mLvvn;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//************************************************************************
//************************************************************************

/** Get the Gauge SURF upright extended descriptor vector of the provided Ipoint	   */
/** Descriptor Dimension 288														   */
void Surf::getGaugeUprightExtendedDescriptor_288(void)
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0;
  int i = 0, j = 0, sample_step = 0, pattern_size = 0, dsize = 0;
  float scale = 0.0, *desc;
  float Lvvp = 0.0, Lwwp = 0.0, mLwwp = 0.0, mLvvp = 0.0;
  float Lvvn = 0.0, Lwwn = 0.0, mLwwn = 0.0, mLvvn = 0.0;
  float lvv = 0.0, lww = 0.0, co = 0.0, si = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0, modg = 0.0;
  
  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 12;
  sample_step = 4;
  
  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 24 s x 24 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvvp=Lwwp=mLvvp=mLwwp=0.0f;
      Lvvn=Lwwn=mLvvn=mLwwn=0.0f;
	  
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y,sample_x,2*fRound(scale));
          ry = haarY(sample_y,sample_x,2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);

          if( modg != 0.0 )
          {
                rxx = haarXX(sample_y,sample_x,2*fRound(scale));
                rxy = haarXY(sample_y,sample_x,2*fRound(scale));
                ryy = haarYY(sample_y,sample_x,2*fRound(scale));

                // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
                lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

                // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
                lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

          if( lww >= 0.0 )
          {
              Lvvp += lvv;
              mLvvp += fabs(lvv);
          }
          else
          {
              Lvvn += lvv;
              mLvvn += fabs(lvv);
          }

          if( lvv >= 0.0 )
          {
              Lwwp += lww;
              mLwwp += fabs(lww);
          }
          else
          {
              Lwwn += lww;
              mLwwn += fabs(lww);
          }

        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lwwn;
      desc[count++] = Lvvn;
      desc[count++] = mLwwn;
      desc[count++] = mLvvn;

      desc[count++] = Lwwp;
      desc[count++] = Lvvp;
      desc[count++] = mLwwp;
      desc[count++] = mLvvp;

      len += Lwwp*Lwwp + Lvvp*Lvvp + mLwwp*mLwwp + mLvvp*mLvvp + Lwwn*Lwwn + 
			 Lvvn*Lvvn + mLwwn*mLwwn + mLvvn*mLvvn;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the Gauge SURF descriptor vector of the provided Ipoint		   				  */
/** Descriptor Dimension 36														      */
void Surf::getGaugeDescriptor_36()
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0;
  int i = 0, j = 0, sample_step = 0, pattern_size = 0, dsize = 0;
  float scale = 0.0, *desc, co = 0.0, si = 0.0;
  float Lvv = 0.0, Lww = 0.0, mLww = 0.0, mLvv = 0.0, lvv = 0.0, lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0, modg = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 9;
  sample_step = 6;

  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 18 s x 18 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvv=Lww=mLvv=mLww=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

		  rx = haarX(sample_y,sample_x,2*fRound(scale));
          ry = haarY(sample_y,sample_x,2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);

          if( modg != 0.0 )
          {
               rxx = haarXX(sample_y,sample_x,2*fRound(scale));
               rxy = haarXY(sample_y,sample_x,2*fRound(scale));
               ryy = haarYY(sample_y,sample_x,2*fRound(scale));

			   // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
			   lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

			   // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
			   lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

                Lww += lww;
                Lvv += lvv;
                mLww += fabs(lww);
                mLvv += fabs(lvv);
        }

      }

      // Add the values to the descriptor vector
      desc[count++] = Lww;
      desc[count++] = Lvv;
      desc[count++] = mLww;
      desc[count++] = mLvv;

      len += Lww*Lww + Lvv*Lvv + mLww*mLww + mLvv*mLvv;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; ++i)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the Gauge SURF upright descriptor vector of the provided Ipoint		 		  */
/** Descriptor Dimension 36														      */
void Surf::getGaugeUprightDescriptor_36()
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, sample_step = 0, pattern_size = 0;
  float scale = 0.0, *desc;
  float Lvv = 0.0, Lww = 0.0, mLww = 0.0, mLvv = 0.0, lvv = 0.0, lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0, modg = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 9;
  sample_step = 6;
  
  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 18 s x 18 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvv=Lww=mLvv=mLww=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + k*scale);
          sample_y = fRound(y + l*scale);

          rx = haarX(sample_y,sample_x,2.0*fRound(scale));
          ry = haarY(sample_y,sample_x,2.0*fRound(scale));
		  
          modg = pow(rx,2) + pow(ry,2);
		  
          if( modg != 0.0 )
          {    
			   rxx =  haarXX(sample_y,sample_x,2*fRound(scale));
			   rxy =  haarXY(sample_y,sample_x,2*fRound(scale));
			   ryy =  haarYY(sample_y,sample_x,2*fRound(scale));

			   // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
			   lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

			   // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
			   lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }
	
           Lww += lww;
           Lvv += lvv;
		   mLww += 1.0*fabs(lww);
           mLvv += 1.0*fabs(lvv);
        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lww;
      desc[count++] = Lvv;
      desc[count++] = mLww;
      desc[count++] = mLvv;
	  
      len += Lww*Lww + Lvv*Lvv + mLww*mLww + mLvv*mLvv;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; ++i)
    desc[i] /= len;
	
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the Gauge SURF extended descriptor vector of the provided Ipoint		  	   */
/** Descriptor Dimension 72														       */
void Surf::getGaugeExtendedDescriptor_72(void)
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, pattern_size = 0, sample_step = 0;
  float scale = 0.0, *desc, co = 0.0, si = 0.0;
  float Lvvp = 0.0, Lwwp = 0.0, mLwwp = 0.0, mLvvp = 0.0;
  float Lvvn = 0.0, Lwwn = 0.0, mLwwn = 0.0, mLvvn = 0.0;
  float lvv = 0.0, lww = 0.0, modg = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 9;
  sample_step = 6;
  
  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 18 s x 18 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvvp=Lwwp=mLvvp=mLwwp=0.0f;
      Lvvn=Lwwn=mLvvn=mLwwn=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y,sample_x,2*fRound(scale));
          ry = haarY(sample_y,sample_x,2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);

          if( modg != 0.0 )
          {    
                rxx = haarXX(sample_y,sample_x,2*fRound(scale));
                rxy = haarXY(sample_y,sample_x,2*fRound(scale));
                ryy = haarYY(sample_y,sample_x,2*fRound(scale));

                // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
                lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

                // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
                lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

          if( lww >= 0.0 )
          {
              Lvvp += lvv;
              mLvvp += fabs(lvv);
          }
          else
          {
              Lvvn += lvv;
              mLvvn += fabs(lvv);
          }

          if( lvv >= 0.0 )
          {
              Lwwp += lww;
              mLwwp += fabs(lww);
          }
          else
          {
              Lwwn += lww;
              mLwwn += fabs(lww);
          }

        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lwwn;
      desc[count++] = Lvvn;
      desc[count++] = mLwwn;
      desc[count++] = mLvvn;

      desc[count++] = Lwwp;
      desc[count++] = Lvvp;
      desc[count++] = mLwwp;
      desc[count++] = mLvvp;

      len += Lwwp*Lwwp + Lvvp*Lvvp + mLwwp*mLwwp + mLvvp*mLvvp + Lwwn*Lwwn + 
			 Lvvn*Lvvn + mLwwn*mLwwn + mLvvn*mLvvn;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//************************************************************************
//************************************************************************

/** Get the Gauge SURF upright extended descriptor vector of the provided Ipoint	   */
/** Descriptor Dimension 72														       */
void Surf::getGaugeUprightExtendedDescriptor_72(void)
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, sample_step = 0, pattern_size = 0;
  float scale = 0.0, *desc, co = 0.0, si = 0.0;
  float Lvvp = 0.0, Lwwp = 0.0, mLwwp = 0.0, mLvvp = 0.0;
  float Lvvn = 0.0, Lwwn = 0.0, mLwwn = 0.0, mLvvn = 0.0;
  float lvv = 0.0, lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0, modg = 0.0;
  
  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 9;
  sample_step = 6;
  
  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 18 s x 18 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvvp=Lwwp=mLvvp=mLwwp=0.0f;
      Lvvn=Lwwn=mLvvn=mLwwn=0.0f;
	  
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y,sample_x,2*fRound(scale));
          ry = haarY(sample_y,sample_x,2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);

          if( modg != 0.0 )
          {
                rxx = haarXX(sample_y,sample_x,2*fRound(scale));
                rxy = haarXY(sample_y,sample_x,2*fRound(scale));
                ryy = haarYY(sample_y,sample_x,2*fRound(scale));

                // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
                lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

                // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
                lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

          if( lww >= 0.0 )
          {
              Lvvp += lvv;
              mLvvp += fabs(lvv);
          }
          else
          {
              Lvvn += lvv;
              mLvvn += fabs(lvv);
          }

          if( lvv >= 0.0 )
          {
              Lwwp += lww;
              mLwwp += fabs(lww);
          }
          else
          {
              Lwwn += lww;
              mLwwn += fabs(lww);
          }

        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lwwn;
      desc[count++] = Lvvn;
      desc[count++] = mLwwn;
      desc[count++] = mLvvn;

      desc[count++] = Lwwp;
      desc[count++] = Lvvp;
      desc[count++] = mLwwp;
      desc[count++] = mLvvp;

      len += Lwwp*Lwwp + Lvvp*Lvvp + mLwwp*mLwwp + mLvvp*mLvvp + Lwwn*Lwwn + 
			 Lvvn*Lvvn + mLwwn*mLwwn + mLvvn*mLvvn;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;
		
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the Gauge SURF descriptor vector of the provided Ipoint		   				  */
/** Descriptor Dimension 100														  */
void Surf::getGaugeDescriptor_100()
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0;
  int i = 0, j = 0, sample_step = 0, pattern_size = 0, dsize = 0;
  float scale = 0.0, *desc, co = 0.0, si = 0.0;
  float Lvv = 0.0, Lww = 0.0, mLww = 0.0, mLvv = 0.0, lvv = 0.0, lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0, modg = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 10;
  sample_step = 4;

  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 20 s x 20 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvv=Lww=mLvv=mLww=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

		  rx = haarX(sample_y,sample_x,2*fRound(scale));
          ry = haarY(sample_y,sample_x,2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);

          if( modg != 0.0 )
          {
               rxx = haarXX(sample_y,sample_x,2*fRound(scale));
               rxy = haarXY(sample_y,sample_x,2*fRound(scale));
               ryy = haarYY(sample_y,sample_x,2*fRound(scale));

			   // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
			   lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

			   // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
			   lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

                Lww += lww;
                Lvv += lvv;
                mLww += fabs(lww);
                mLvv += fabs(lvv);
        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lww;
      desc[count++] = Lvv;
      desc[count++] = mLww;
      desc[count++] = mLvv;

      len += Lww*Lww + Lvv*Lvv + mLww*mLww + mLvv*mLvv;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; ++i)
    desc[i] /= len;
	
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the Gauge SURF upright descriptor vector of the provided Ipoint		 		  */
/** Descriptor Dimension 100														  */
void Surf::getGaugeUprightDescriptor_100()
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0;
  int i = 0, j = 0, sample_step = 0, pattern_size = 0, dsize = 0;
  float *desc, scale = 0.0;
  float Lvv = 0.0, Lww = 0.0, mLww = 0.0, mLvv = 0.0, lvv = 0.0, lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0, modg = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 10;
  sample_step = 4;
  
  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 20 s x 20 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvv=Lww=mLvv=mLww=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + k*scale);
          sample_y = fRound(y + l*scale);

          rx = haarX(sample_y,sample_x,2.0*fRound(scale));
          ry = haarY(sample_y,sample_x,2.0*fRound(scale));
		  
          modg = pow(rx,2) + pow(ry,2);
		  
          if( modg != 0.0 )
          {    
			   rxx =  haarXX(sample_y,sample_x,2*fRound(scale));
			   rxy =  haarXY(sample_y,sample_x,2*fRound(scale));
			   ryy =  haarYY(sample_y,sample_x,2*fRound(scale));

			   // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
			   lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

			   // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
			   lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }
	
           Lww += lww;
           Lvv += lvv;
		   mLww += 1.0*fabs(lww);
           mLvv += 1.0*fabs(lvv);
        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lww;
      desc[count++] = Lvv;
      desc[count++] = mLww;
      desc[count++] = mLvv;
	  
      len += Lww*Lww + Lvv*Lvv + mLww*mLww + mLvv*mLvv;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; ++i)
    desc[i] /= len;
	
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the Gauge SURF extended descriptor vector of the provided Ipoint		  	   */
/** Descriptor Dimension 200														   */
void Surf::getGaugeExtendedDescriptor_200(void)
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, pattern_size = 0, sample_step = 0;
  float scale = 0.0, *desc, co = 0.0, si = 0.0;
  float Lvvp = 0.0, Lwwp = 0.0, mLwwp = 0.0, mLvvp = 0.0;
  float Lvvn = 0.0, Lwwn = 0.0, mLwwn = 0.0, mLvvn = 0.0;
  float lvv = 0.0, lww = 0.0, modg = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 10;
  sample_step = 4;
  
  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 20 s x 20 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvvp=Lwwp=mLvvp=mLwwp=0.0f;
      Lvvn=Lwwn=mLvvn=mLwwn=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y,sample_x,2*fRound(scale));
          ry = haarY(sample_y,sample_x,2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);

          if( modg != 0.0 )
          {    
                rxx = haarXX(sample_y,sample_x,2*fRound(scale));
                rxy = haarXY(sample_y,sample_x,2*fRound(scale));
                ryy = haarYY(sample_y,sample_x,2*fRound(scale));

                // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
                lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

                // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
                lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

          if( lww >= 0.0 )
          {
              Lvvp += lvv;
              mLvvp += fabs(lvv);
          }
          else
          {
              Lvvn += lvv;
              mLvvn += fabs(lvv);
          }

          if( lvv >= 0.0 )
          {
              Lwwp += lww;
              mLwwp += fabs(lww);
          }
          else
          {
              Lwwn += lww;
              mLwwn += fabs(lww);
          }

        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lwwn;
      desc[count++] = Lvvn;
      desc[count++] = mLwwn;
      desc[count++] = mLvvn;

      desc[count++] = Lwwp;
      desc[count++] = Lvvp;
      desc[count++] = mLwwp;
      desc[count++] = mLvvp;

      len += Lwwp*Lwwp + Lvvp*Lvvp + mLwwp*mLwwp + mLvvp*mLvvp + Lwwn*Lwwn + 
			 Lvvn*Lvvn + mLwwn*mLwwn + mLvvn*mLvvn;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;
		
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//************************************************************************
//************************************************************************

/** Get the Gauge SURF upright extended descriptor vector of the provided Ipoint	   */
/** Descriptor Dimension 200														   */
void Surf::getGaugeUprightExtendedDescriptor_200(void)
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count = 0, dsize = 0;
  int i = 0, j = 0, sample_step = 0, pattern_size = 0;
  float scale = 0.0, *desc, co = 0.0, si = 0.0;
  float Lvvp = 0.0, Lwwp = 0.0, mLwwp = 0.0, mLvvp = 0.0;
  float Lvvn = 0.0, Lwwn = 0.0, mLwwn = 0.0, mLvvn = 0.0;
  float lvv = 0.0, lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, ryy = 0.0, rxy = 0.0, len = 0.0, modg = 0.0;
  
  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 10;
  sample_step = 4;
  
  i = j = -pattern_size;

  // Calculate descriptor for this interest point
  // Area of size 18 s x 18 s
  while(i < pattern_size)
  {
    j = -pattern_size;

    while(j < pattern_size)
    {
      Lvvp=Lwwp=mLvvp=mLwwp=0.0f;
      Lvvn=Lwwn=mLvvn=mLwwn=0.0f;
	  
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y,sample_x,2*fRound(scale));
          ry = haarY(sample_y,sample_x,2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);

          if( modg != 0.0 )
          {
                rxx = haarXX(sample_y,sample_x,2*fRound(scale));
                rxy = haarXY(sample_y,sample_x,2*fRound(scale));
                ryy = haarYY(sample_y,sample_x,2*fRound(scale));

                // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
                lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

                // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
                lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

          if( lww >= 0.0 )
          {
              Lvvp += lvv;
              mLvvp += fabs(lvv);
          }
          else
          {
              Lvvn += lvv;
              mLvvn += fabs(lvv);
          }

          if( lvv >= 0.0 )
          {
              Lwwp += lww;
              mLwwp += fabs(lww);
          }
          else
          {
              Lwwn += lww;
              mLwwn += fabs(lww);
          }

        }
      }

      // Add the values to the descriptor vector
      desc[count++] = Lwwn;
      desc[count++] = Lvvn;
      desc[count++] = mLwwn;
      desc[count++] = mLvvn;

      desc[count++] = Lwwp;
      desc[count++] = Lvvp;
      desc[count++] = mLwwp;
      desc[count++] = mLvvp;

      len += Lwwp*Lwwp + Lvvp*Lvvp + mLwwp*mLwwp + mLvvp*mLvvp + Lwwn*Lwwn + 
			 Lvvn*Lvvn + mLwwn*mLwwn + mLvvn*mLvvn;

      j += sample_step;
    }

      i += sample_step;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;
		
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the descriptor vector of the provided Ipoint								   */
/** SURF descriptor without Gaussian Weighting										   */
/** Descriptor Dimension 64															   */
void Surf::getDescriptor_NW_64()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0, co = 0.0, si = 0.0;
  float rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len=0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  dsize = ipt->size;
  desc = ipt->descriptor;
  sample_step = 5;
  pattern_size = 10;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dx=dy=mdx=mdy=0.0;
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          // Get the gaussian weighted x and y responses on rotated axis
          rrx = -rx*si + ry*co;
          rry = rx*co + ry*si;

          dx += rrx;
          dy += rry;
          mdx += fabs(rrx);
          mdy += fabs(rry);
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dx;
      desc[count++] = dy;
      desc[count++] = mdx;
      desc[count++] = mdy;

      // store the current length^2 of the vector
      len += dx*dx + dy*dy + mdx*mdx + mdy*mdy;
    }
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the upright descriptor vector of the provided Ipoint						   */
/** SURF descriptor without Gaussian Weighting  									   */
/** Descriptor Dimension 64															   */
void Surf::getUprightDescriptor_NW_64()
{
  int y = 0, x = 0, dsize = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0;
  float rx = 0.0, ry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  y = fRound(ipt->y);
  x = fRound(ipt->x);
  dsize = ipt->size;
  desc = ipt->descriptor;  
  sample_step = 5;
  pattern_size = 10;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dx=dy=mdx=mdy=0.0;
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          rx = haarX(fRound(k*scale+y), fRound(l*scale+x), 2*fRound(scale));
          ry = haarY(fRound(k*scale+y), fRound(l*scale+x), 2*fRound(scale));

          dx += rx;
          dy += ry;
          mdx += fabs(rx);
          mdy += fabs(ry);
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dx;
      desc[count++] = dy;
      desc[count++] = mdx;
      desc[count++] = mdy;

      // store the current length^2 of the vector
      len += dx*dx + dy*dy + mdx*mdx + mdy*mdy;
    }
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the 128 descriptor vector of the provided Ipoint						   */
/** SURF descriptor without Gaussian Weighting  								   */
/** Descriptor Dimension 128													   */
void Surf::getExtendedDescriptor_NW_128()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float co = 0.0, si = 0.0, dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 10;
  sample_step = 5;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dxp=dyp=mdxp=mdyp=0.0f;
      dxn=dyn=mdxn=mdyn=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          // Get the gaussian weighted x and y responses on rotated axis
          rrx = -rx*si + ry*co;
          rry = rx*co + ry*si;

          if( rry >= 0.0 )
          {
              dxp += rrx;
              mdxp += fabs(rrx);
          }
          else
          {
              dxn += rrx;
              mdxn += fabs(rrx);
          }

          if( rrx >= 0.0 )
          {
              dyp += rry;
              mdyp += fabs(rry);
          }
          else
          {
              dyn += rry;
              mdyn += fabs(rry);
          }
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dxp;
      desc[count++] = mdxp;
      desc[count++] = dxn;
      desc[count++] = mdxn;

      desc[count++] = dyp;
      desc[count++] = mdyp;
      desc[count++] = dyn;
      desc[count++] = mdyn;
      }

      // store the current length^2 of the vector
      len += dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn;

  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the 128 upright descriptor vector of the provided Ipoint					   */
/** SURF descriptor without Gaussian Weighting  									   */
/** Descriptor Dimension 128														   */
void Surf::getUprightExtendedDescriptor_NW_128()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float rx = 0.0, ry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 10;
  sample_step = 5;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dxp=dyp=mdxp=mdyp=0.0;
      dxn=dyn=mdxn=mdyn=0.0;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + k*scale);
          sample_y = fRound(y + l*scale);

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          if( ry >= 0.0 )
          {
              dxp += rx;
              mdxp += fabs(rx);
          }
          else
          {
              dxn += rx;
              mdxn += fabs(rx);
          }

          if( rx >= 0.0 )
          {
              dyp += ry;
              mdyp += fabs(ry);
          }
          else
          {
              dyn += ry;
              mdyn += fabs(ry);
          }
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dxp;
      desc[count++] = mdxp;
      desc[count++] = dxn;
      desc[count++] = mdxn;

      desc[count++] = dyp;
      desc[count++] = mdyp;
      desc[count++] = dyn;
      desc[count++] = mdyn;
      }

      // store the current length^2 of the vector
      len += dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn;
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the descriptor vector of the provided Ipoint								   */
/** SURF descriptor without Gaussian Weighting										   */
/** Descriptor Dimension 16															   */
void Surf::getDescriptor_NW_16()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0, co = 0.0, si = 0.0;
  float rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len=0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  dsize = ipt->size;
  desc = ipt->descriptor;
  pattern_size = 5;
  sample_step = 5;
	
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dx=dy=mdx=mdy=0.0;
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          // Get the gaussian weighted x and y responses on rotated axis
          rrx = -rx*si + ry*co;
          rry = rx*co + ry*si;

          dx += rrx;
          dy += rry;
          mdx += fabs(rrx);
          mdy += fabs(rry);
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dx;
      desc[count++] = dy;
      desc[count++] = mdx;
      desc[count++] = mdy;

      // store the current length^2 of the vector
      len += dx*dx + dy*dy + mdx*mdx + mdy*mdy;
    }
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the upright descriptor vector of the provided Ipoint						   */
/** SURF descriptor without Gaussian Weighting  									   */
/** Descriptor Dimension 16															   */
void Surf::getUprightDescriptor_NW_16()
{
  int y = 0, x = 0, dsize = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0;
  float rx = 0.0, ry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  y = fRound(ipt->y);
  x = fRound(ipt->x);
  dsize = ipt->size;
  desc = ipt->descriptor;  
  pattern_size = 5;
  sample_step = 5;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dx=dy=mdx=mdy=0.0;
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          rx = haarX(fRound(k*scale+y), fRound(l*scale+x), 2*fRound(scale));
          ry = haarY(fRound(k*scale+y), fRound(l*scale+x), 2*fRound(scale));

          dx += rx;
          dy += ry;
          mdx += fabs(rx);
          mdy += fabs(ry);
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dx;
      desc[count++] = dy;
      desc[count++] = mdx;
      desc[count++] = mdy;

      // store the current length^2 of the vector
      len += dx*dx + dy*dy + mdx*mdx + mdy*mdy;
    }
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the 32 descriptor vector of the provided Ipoint						       */
/** SURF descriptor without Gaussian Weighting  								   */
/** Descriptor Dimension 32													       */
void Surf::getExtendedDescriptor_NW_32()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float co = 0.0, si = 0.0, dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 10;
  sample_step = 10;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dxp=dyp=mdxp=mdyp=0.0f;
      dxn=dyn=mdxn=mdyn=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          // Get the gaussian weighted x and y responses on rotated axis
          rrx = -rx*si + ry*co;
          rry = rx*co + ry*si;

          if( rry >= 0.0 )
          {
              dxp += rrx;
              mdxp += fabs(rrx);
          }
          else
          {
              dxn += rrx;
              mdxn += fabs(rrx);
          }

          if( rrx >= 0.0 )
          {
              dyp += rry;
              mdyp += fabs(rry);
          }
          else
          {
              dyn += rry;
              mdyn += fabs(rry);
          }
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dxp;
      desc[count++] = mdxp;
      desc[count++] = dxn;
      desc[count++] = mdxn;

      desc[count++] = dyp;
      desc[count++] = mdyp;
      desc[count++] = dyn;
      desc[count++] = mdyn;
      }

      // store the current length^2 of the vector
      len += dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn;

  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

}

//**************************************************************************************
//**************************************************************************************

/** Get the 32 upright descriptor vector of the provided Ipoint					       */
/** SURF descriptor without Gaussian Weighting  									   */
/** Descriptor Dimension 32														       */
void Surf::getUprightExtendedDescriptor_NW_32()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float rx = 0.0, ry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 10;
  sample_step = 10;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dxp=dyp=mdxp=mdyp=0.0;
      dxn=dyn=mdxn=mdyn=0.0;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + k*scale);
          sample_y = fRound(y + l*scale);

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          if( ry >= 0.0 )
          {
              dxp += rx;
              mdxp += fabs(rx);
          }
          else
          {
              dxn += rx;
              mdxn += fabs(rx);
          }

          if( rx >= 0.0 )
          {
              dyp += ry;
              mdyp += fabs(ry);
          }
          else
          {
              dyn += ry;
              mdyn += fabs(ry);
          }
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dxp;
      desc[count++] = mdxp;
      desc[count++] = dxn;
      desc[count++] = mdxn;

      desc[count++] = dyp;
      desc[count++] = mdyp;
      desc[count++] = dyn;
      desc[count++] = mdyn;
      }

      // store the current length^2 of the vector
      len += dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn;
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the descriptor vector of the provided Ipoint								   */
/** SURF descriptor without Gaussian Weighting										   */
/** Descriptor Dimension 36															   */
void Surf::getDescriptor_NW_36()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0, co = 0.0, si = 0.0;
  float rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len=0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  dsize = ipt->size;
  desc = ipt->descriptor;
  pattern_size = 9;
  sample_step = 6;
	
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dx=dy=mdx=mdy=0.0;
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          // Get the gaussian weighted x and y responses on rotated axis
          rrx = -rx*si + ry*co;
          rry = rx*co + ry*si;

          dx += rrx;
          dy += rry;
          mdx += fabs(rrx);
          mdy += fabs(rry);
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dx;
      desc[count++] = dy;
      desc[count++] = mdx;
      desc[count++] = mdy;

      // store the current length^2 of the vector
      len += dx*dx + dy*dy + mdx*mdx + mdy*mdy;
    }
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the upright descriptor vector of the provided Ipoint						   */
/** SURF descriptor without Gaussian Weighting  									   */
/** Descriptor Dimension 36															   */
void Surf::getUprightDescriptor_NW_36()
{
  int y = 0, x = 0, dsize = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0;
  float rx = 0.0, ry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  y = fRound(ipt->y);
  x = fRound(ipt->x);
  dsize = ipt->size;
  desc = ipt->descriptor;  
  pattern_size = 9;
  sample_step = 6;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dx=dy=mdx=mdy=0.0;
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          rx = haarX(fRound(k*scale+y), fRound(l*scale+x), 2*fRound(scale));
          ry = haarY(fRound(k*scale+y), fRound(l*scale+x), 2*fRound(scale));

          dx += rx;
          dy += ry;
          mdx += fabs(rx);
          mdy += fabs(ry);
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dx;
      desc[count++] = dy;
      desc[count++] = mdx;
      desc[count++] = mdy;

      // store the current length^2 of the vector
      len += dx*dx + dy*dy + mdx*mdx + mdy*mdy;
    }
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the 32 descriptor vector of the provided Ipoint						       */
/** SURF descriptor without Gaussian Weighting  								   */
/** Descriptor Dimension 72													       */
void Surf::getExtendedDescriptor_NW_72()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float co = 0.0, si = 0.0, dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 9;
  sample_step = 6;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dxp=dyp=mdxp=mdyp=0.0f;
      dxn=dyn=mdxn=mdyn=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          // Get the gaussian weighted x and y responses on rotated axis
          rrx = -rx*si + ry*co;
          rry = rx*co + ry*si;

          if( rry >= 0.0 )
          {
              dxp += rrx;
              mdxp += fabs(rrx);
          }
          else
          {
              dxn += rrx;
              mdxn += fabs(rrx);
          }

          if( rrx >= 0.0 )
          {
              dyp += rry;
              mdyp += fabs(rry);
          }
          else
          {
              dyn += rry;
              mdyn += fabs(rry);
          }
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dxp;
      desc[count++] = mdxp;
      desc[count++] = dxn;
      desc[count++] = mdxn;

      desc[count++] = dyp;
      desc[count++] = mdyp;
      desc[count++] = dyn;
      desc[count++] = mdyn;
      }

      // store the current length^2 of the vector
      len += dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn;

  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the 72 upright descriptor vector of the provided Ipoint					       */
/** SURF descriptor without Gaussian Weighting  									   */
/** Descriptor Dimension 72														       */
void Surf::getUprightExtendedDescriptor_NW_72()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float rx = 0.0, ry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 9;
  sample_step = 6;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dxp=dyp=mdxp=mdyp=0.0;
      dxn=dyn=mdxn=mdyn=0.0;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + k*scale);
          sample_y = fRound(y + l*scale);

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          if( ry >= 0.0 )
          {
              dxp += rx;
              mdxp += fabs(rx);
          }
          else
          {
              dxn += rx;
              mdxn += fabs(rx);
          }

          if( rx >= 0.0 )
          {
              dyp += ry;
              mdyp += fabs(ry);
          }
          else
          {
              dyn += ry;
              mdyn += fabs(ry);
          }
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dxp;
      desc[count++] = mdxp;
      desc[count++] = dxn;
      desc[count++] = mdxn;

      desc[count++] = dyp;
      desc[count++] = mdyp;
      desc[count++] = dyn;
      desc[count++] = mdyn;
      }

      // store the current length^2 of the vector
      len += dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn;
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the descriptor vector of the provided Ipoint								   */
/** SURF descriptor without Gaussian Weighting										   */
/** Descriptor Dimension 100														   */
void Surf::getDescriptor_NW_100()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0, co = 0.0, si = 0.0;
  float rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len=0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  dsize = ipt->size;
  desc = ipt->descriptor;
  pattern_size = 10;
  sample_step = 4;
	
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dx=dy=mdx=mdy=0.0;
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          // Get the gaussian weighted x and y responses on rotated axis
          rrx = -rx*si + ry*co;
          rry = rx*co + ry*si;

          dx += rrx;
          dy += rry;
          mdx += fabs(rrx);
          mdy += fabs(rry);
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dx;
      desc[count++] = dy;
      desc[count++] = mdx;
      desc[count++] = mdy;

      // store the current length^2 of the vector
      len += dx*dx + dy*dy + mdx*mdx + mdy*mdy;
    }
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the upright descriptor vector of the provided Ipoint						   */
/** SURF descriptor without Gaussian Weighting  									   */
/** Descriptor Dimension 100														   */
void Surf::getUprightDescriptor_NW_100()
{
  int y = 0, x = 0, dsize = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0;
  float rx = 0.0, ry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  y = fRound(ipt->y);
  x = fRound(ipt->x);
  dsize = ipt->size;
  desc = ipt->descriptor;  
  pattern_size = 10;
  sample_step = 4;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dx=dy=mdx=mdy=0.0;
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          rx = haarX(fRound(k*scale+y), fRound(l*scale+x), 2*fRound(scale));
          ry = haarY(fRound(k*scale+y), fRound(l*scale+x), 2*fRound(scale));

          dx += rx;
          dy += ry;
          mdx += fabs(rx);
          mdy += fabs(ry);
        }
      }
	  
      // add the values to the descriptor vector
      desc[count++] = dx;
      desc[count++] = dy;
      desc[count++] = mdx;
      desc[count++] = mdy;

      // store the current length^2 of the vector
      len += dx*dx + dy*dy + mdx*mdx + mdy*mdy;
    }
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;
  
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the 200 descriptor vector of the provided Ipoint						       */
/** SURF descriptor without Gaussian Weighting  								   	   */
/** Descriptor Dimension 200													       */
void Surf::getExtendedDescriptor_NW_200()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float co = 0.0, si = 0.0, dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 10;
  sample_step = 4;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dxp=dyp=mdxp=mdyp=0.0f;
      dxn=dyn=mdxn=mdyn=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          // Get the gaussian weighted x and y responses on rotated axis
          rrx = -rx*si + ry*co;
          rry = rx*co + ry*si;

          if( rry >= 0.0 )
          {
              dxp += rrx;
              mdxp += fabs(rrx);
          }
          else
          {
              dxn += rrx;
              mdxn += fabs(rrx);
          }

          if( rrx >= 0.0 )
          {
              dyp += rry;
              mdyp += fabs(rry);
          }
          else
          {
              dyn += rry;
              mdyn += fabs(rry);
          }
        }
      }

		  // add the values to the descriptor vector
		  desc[count++] = dxp;
		  desc[count++] = mdxp;
		  desc[count++] = dxn;
		  desc[count++] = mdxn;

		  desc[count++] = dyp;
		  desc[count++] = mdyp;
		  desc[count++] = dyn;
		  desc[count++] = mdyn;
      }

      // store the current length^2 of the vector
      len += dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn;

  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the 200 upright descriptor vector of the provided Ipoint					   */
/** SURF descriptor without Gaussian Weighting  									   */
/** Descriptor Dimension 200														   */
void Surf::getUprightExtendedDescriptor_NW_200()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float rx = 0.0, ry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 10;
  sample_step = 4;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dxp=dyp=mdxp=mdyp=0.0;
      dxn=dyn=mdxn=mdyn=0.0;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + k*scale);
          sample_y = fRound(y + l*scale);

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          if( ry >= 0.0 )
          {
              dxp += rx;
              mdxp += fabs(rx);
          }
          else
          {
              dxn += rx;
              mdxn += fabs(rx);
          }

          if( rx >= 0.0 )
          {
              dyp += ry;
              mdyp += fabs(ry);
          }
          else
          {
              dyn += ry;
              mdyn += fabs(ry);
          }
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dxp;
      desc[count++] = mdxp;
      desc[count++] = dxn;
      desc[count++] = mdxn;

      desc[count++] = dyp;
      desc[count++] = mdyp;
      desc[count++] = dyn;
      desc[count++] = mdyn;
      }

      // store the current length^2 of the vector
      len += dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn;
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the descriptor vector of the provided Ipoint								   */
/** SURF descriptor without Gaussian Weighting										   */
/** Descriptor Dimension 144														   */
void Surf::getDescriptor_NW_144()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0, co = 0.0, si = 0.0;
  float rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len=0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  dsize = ipt->size;
  desc = ipt->descriptor;
  pattern_size = 12;
  sample_step = 4;
	
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dx=dy=mdx=mdy=0.0;
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          // Get the gaussian weighted x and y responses on rotated axis
          rrx = -rx*si + ry*co;
          rry = rx*co + ry*si;

          dx += rrx;
          dy += rry;
          mdx += fabs(rrx);
          mdy += fabs(rry);
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dx;
      desc[count++] = dy;
      desc[count++] = mdx;
      desc[count++] = mdy;

      // store the current length^2 of the vector
      len += dx*dx + dy*dy + mdx*mdx + mdy*mdy;
    }
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the upright descriptor vector of the provided Ipoint						   */
/** SURF descriptor without Gaussian Weighting  									   */
/** Descriptor Dimension 144														   */
void Surf::getUprightDescriptor_NW_144()
{
  int y = 0, x = 0, dsize = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0;
  float rx = 0.0, ry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  y = fRound(ipt->y);
  x = fRound(ipt->x);
  dsize = ipt->size;
  desc = ipt->descriptor;  
  pattern_size = 12;
  sample_step = 4;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dx=dy=mdx=mdy=0.0;
      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          rx = haarX(fRound(k*scale+y), fRound(l*scale+x), 2*fRound(scale));
          ry = haarY(fRound(k*scale+y), fRound(l*scale+x), 2*fRound(scale));

          dx += rx;
          dy += ry;
          mdx += fabs(rx);
          mdy += fabs(ry);
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dx;
      desc[count++] = dy;
      desc[count++] = mdx;
      desc[count++] = mdy;

      // store the current length^2 of the vector
      len += dx*dx + dy*dy + mdx*mdx + mdy*mdy;
    }
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the extended descriptor vector of the provided Ipoint						   */
/** SURF descriptor without Gaussian Weighting  								   	   */
/** Descriptor Dimension 288													       */
void Surf::getExtendedDescriptor_NW_288()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float co = 0.0, si = 0.0, dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 12;
  sample_step = 4;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dxp=dyp=mdxp=mdyp=0.0f;
      dxn=dyn=mdxn=mdyn=0.0f;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          // Get the gaussian weighted x and y responses on rotated axis
          rrx = -rx*si + ry*co;
          rry = rx*co + ry*si;

          if( rry >= 0.0 )
          {
              dxp += rrx;
              mdxp += fabs(rrx);
          }
          else
          {
              dxn += rrx;
              mdxn += fabs(rrx);
          }

          if( rrx >= 0.0 )
          {
              dyp += rry;
              mdyp += fabs(rry);
          }
          else
          {
              dyn += rry;
              mdyn += fabs(rry);
          }
        }
      }
		  // add the values to the descriptor vector
		  desc[count++] = dxp;
		  desc[count++] = mdxp;
		  desc[count++] = dxn;
		  desc[count++] = mdxn;

		  desc[count++] = dyp;
		  desc[count++] = mdyp;
		  desc[count++] = dyn;
		  desc[count++] = mdyn;
      }

      // store the current length^2 of the vector
      len += dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn;

  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the extended upright descriptor vector of the provided Ipoint				   */
/** SURF descriptor without Gaussian Weighting  									   */
/** Descriptor Dimension 288														   */
void Surf::getUprightExtendedDescriptor_NW_288()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count = 0;
  int sample_step = 0, pattern_size = 0;
  float *desc, scale = 0.0, dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float rx = 0.0, ry = 0.0, len = 0.0;

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  desc = ipt->descriptor;
  dsize = ipt->size;
  pattern_size = 12;
  sample_step = 4;
  
  // Calculate descriptor for this interest point
  for (int i = -pattern_size; i < pattern_size; i+=sample_step)
  {
    for (int j = -pattern_size; j < pattern_size; j+=sample_step)
    {
      dxp=dyp=mdxp=mdyp=0.0;
      dxn=dyn=mdxn=mdyn=0.0;

      for (int k = i; k < i + sample_step; ++k)
      {
        for (int l = j; l < j + sample_step; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + k*scale);
          sample_y = fRound(y + l*scale);

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          if( ry >= 0.0 )
          {
              dxp += rx;
              mdxp += fabs(rx);
          }
          else
          {
              dxn += rx;
              mdxn += fabs(rx);
          }

          if( rx >= 0.0 )
          {
              dyp += ry;
              mdyp += fabs(ry);
          }
          else
          {
              dyn += ry;
              mdyn += fabs(ry);
          }
        }
      }

      // add the values to the descriptor vector
      desc[count++] = dxp;
      desc[count++] = mdxp;
      desc[count++] = dxn;
      desc[count++] = mdxn;

      desc[count++] = dyp;
      desc[count++] = mdyp;
      desc[count++] = dyn;
      desc[count++] = mdyn;
      }

      // store the current length^2 of the vector
      len += dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn;
  }

  // convert to unit vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
        desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}


//**************************************************************************************
//**************************************************************************************

/** Get the modified gauge SURF descriptor vector of the provided Ipoint			   */
/** See Agrawal et al. (ECCV08) for more details about the M-SURF square grid		   */
/** Descriptor Dimension 64															   */
void Surf::getMGaugeDescriptor_64()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count=0;
  int sample_step = 0, pattern_size = 0, i = 0, ix = 0, j = 0, jx = 0, xs = 0, ys = 0;
  float *desc, scale = 0.0, gauss_s1 = 0.0, gauss_s2 = 0.0, co = 0.0, si = 0.0;
  float lvv = 0.0, lww = 0.0, mLvv = 0.0, mLww = 0.0, Lvv = 0.0, Lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, rxy = 0.0, ryy = 0.0, modg = 0.0, len = 0.0;
  float cx = -0.5, cy = 0.5; //Subregion centers for the 4x4 gaussian weighting

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size ;
  sample_step = 5;
  pattern_size = 12;
  
  i = -8;

  // Calculate descriptor for this interest point
  // Area of size 24 s x 24 s
  while(i < pattern_size)
  {
    j = -8;
    i = i-4;

    cx += 1.f;
    cy = -0.5f;

    while(j < 12)
    {
      Lvv=Lww=mLvv=mLww=0.0f;
      cy += 1.f;

      j = j - 4;

      ix = i + sample_step;
      jx = j + sample_step;

      xs = fRound(x + ( -jx*scale*si + ix*scale*co));
      ys = fRound(y + ( jx*scale*co + ix*scale*si));

      for (int k = i; k < i + 9; ++k)
      {
        for (int l = j; l < j + 9; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          // Get the gaussian weighted x and y responses
          gauss_s1 = gaussian(xs-sample_x,ys-sample_y,2.5f*scale);

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

		  modg = pow(rx,2) + pow(ry,2);
		  
          if( modg != 0.0 )
          {    
			   rxx =  haarXX(sample_y,sample_x,2*fRound(scale));
			   rxy =  haarXY(sample_y,sample_x,2*fRound(scale));
			   ryy =  haarYY(sample_y,sample_x,2*fRound(scale));

			   // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
			   lww = gauss_s1*(pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

			   // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
			   lvv = gauss_s1*(-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }
		  
		  Lww += lww;
          Lvv += lvv;
		  mLww += 1.0*fabs(lww);
          mLvv += 1.0*fabs(lvv);
        }
      }

      // Add the values to the descriptor vector
      gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f);

      desc[count++] = Lww*gauss_s2;
      desc[count++] = Lvv*gauss_s2;
      desc[count++] = mLww*gauss_s2;
      desc[count++] = mLvv*gauss_s2;

      len += (Lww*Lww + Lvv*Lvv + mLww*mLww + mLvv*mLvv)*gauss_s2*gauss_s2;

      j += 9;
    }
    i += 9;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the modified upright G-SURF descriptor vector of the provided Ipoint		   */
/** See Agrawal et al. (ECCV08) for more details about the M-SURF square grid		   */
/** Descriptor Dimension 64															   */
void Surf::getMGaugeUprightDescriptor_64()
{
  int y = 0, x = 0, dsize = 0, sample_x = 0, sample_y = 0, count=0;
  int sample_step = 0, pattern_size = 0, i = 0, ix = 0, j = 0, jx = 0, xs = 0, ys = 0;
  float *desc, scale = 0.0; 
  float gauss_s1 = 0.0, gauss_s2 = 0.0;
  float lvv = 0.0, lww = 0.0, mLvv = 0.0, mLww = 0.0, Lvv = 0.0, Lww = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, rxy = 0.0, ryy = 0.0, modg = 0.0, len = 0.0;
  float cx = -0.5, cy = 0.5; // Subregion centers for the 4x4 gaussian weighting
  
  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  desc = ipt->descriptor;
  dsize = ipt->size;
  sample_step = 5;
  pattern_size = 12;
  
  i = -8;

  // Calculate descriptor for this interest point
  // Area of size 24 s x 24 s
  while(i < pattern_size)
  {
    j = -8;
    i = i-4;

    cx += 1.0f;
    cy = -0.5f;

    while(j < pattern_size)
    {
      Lvv=Lww=mLvv=mLww=0.0f;
      cy += 1.0f;

      j = j - 4;

      ix = i + sample_step;
      jx = j + sample_step;

      xs = fRound(x + (ix*scale));
      ys = fRound(y + (jx*scale));

      for (int k = i; k < i + 9; ++k)
      {
        for (int l = j; l < j + 9; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + k*scale);
          sample_y = fRound(y + l*scale);

          // Get the gaussian weighted x and y responses
          gauss_s1 = gaussian(xs-sample_x,ys-sample_y,2.5f*scale);

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

          modg = pow(rx,2) + pow(ry,2);
		  
          if( modg != 0.0 )
          {    
			   rxx =  haarXX(sample_y,sample_x,2*fRound(scale));
			   rxy =  haarXY(sample_y,sample_x,2*fRound(scale));
			   ryy =  haarYY(sample_y,sample_x,2*fRound(scale));

			   // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
			   lww = gauss_s1*(pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

			   // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
			   lvv = gauss_s1*(-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }
	
		  Lww += lww;
          Lvv += lvv;
		  mLww += 1.0*fabs(lww);
          mLvv += 1.0*fabs(lvv);

        }
      }
      // Add the values to the descriptor vector
      gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f);

      desc[count++] = Lww*gauss_s2;
      desc[count++] = Lvv*gauss_s2;
      desc[count++] = mLww*gauss_s2;
      desc[count++] = mLvv*gauss_s2;

      len += (Lww*Lww + Lvv*Lvv + mLww*mLww + mLvv*mLvv)*gauss_s2*gauss_s2;

      j += 9;
    }

      i += 9;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}


//**************************************************************************************
//**************************************************************************************

/** Get the modified gauge 128 SURF descriptor vector of the provided Ipoint		   */
/** See Agrawal et al. (ECCV08) for more details  									   */
/** Descriptor Dimension 128														   */
void Surf::getMGaugeExtendedDescriptor_128()
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count=0, dsize = 0;
  int i = 0, ix = 0, j = 0, jx = 0, xs = 0, ys = 0;
  float sample_step = 0, pattern_size = 0, scale = 0.0, *desc, co = 0.0, si = 0.0;
  float dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float gauss_s1 = 0.f, gauss_s2 = 0.f, lvv = 0.0, lww = 0.0, modg = 0.0;
  float rx = 0.f, ry = 0.f, rxx = 0.f, ryy = 0.f, rxy = 0.f, len = 0.f;
  float cx = -0.5f, cy = 0.f; // Subregion centers for the 4x4 gaussian weighting

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  co = cos(ipt->orientation);
  si = sin(ipt->orientation);
  desc = ipt->descriptor;
  dsize = ipt->size;
  sample_step = 5;
  pattern_size = 12;
  
  i = -8;

  // Calculate descriptor for this interest point
  // Area of size 24 s x 24 s
  while(i < pattern_size)
  {
    j = -8;
    i = i-4;

    cx += 1.f;
    cy = -0.5f;

    while(j < pattern_size)
    {
      dxp=dyp=mdxp=mdyp=0.0f;
      dxn=dyn=mdxn=mdyn=0.0f;

      cy += 1.f;

      j = j - 4;

      ix = i + sample_step;
      jx = j + sample_step;

      xs = fRound(x + ( -jx*scale*si + ix*scale*co));
      ys = fRound(y + ( jx*scale*co + ix*scale*si));

      for (int k = i; k < i + 9; ++k)
      {
        for (int l = j; l < j + 9; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + (-l*scale*si + k*scale*co));
          sample_y = fRound(y + ( l*scale*co + k*scale*si));

          // Get the gaussian weighted x and y responses
          gauss_s1 = gaussian(xs-sample_x,ys-sample_y,2.5f*scale);

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

		  modg = pow(rx,2) + pow(ry,2);
		  
          if( modg != 0.0 )
          {    
			   rxx =  haarXX(sample_y,sample_x,2*fRound(scale));
			   rxy =  haarXY(sample_y,sample_x,2*fRound(scale));
			   ryy =  haarYY(sample_y,sample_x,2*fRound(scale));

			   // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
			   lww = gauss_s1*(pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

			   // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
			   lvv = gauss_s1*(-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

          if( lvv >= 0.0 )
          {
              dxp += lww;
              mdxp += fabs(lww);
          }
          else
          {
              dxn += lww;
              mdxn += fabs(lww);
          }

          if( lww >= 0.0 )
          {
              dyp += lvv;
              mdyp += fabs(lvv);
          }
          else
          {
              dyn += lvv;
              mdyn += fabs(lvv);
          }
        }
      }

      // Add the values to the descriptor vector
      gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f);

      desc[count++] = dxp*gauss_s2;
      desc[count++] = mdxp*gauss_s2;
      desc[count++] = dxn*gauss_s2;
      desc[count++] = mdxn*gauss_s2;

      desc[count++] = dyp*gauss_s2;
      desc[count++] = mdyp*gauss_s2;
      desc[count++] = dyn*gauss_s2;
      desc[count++] = mdyn*gauss_s2;

      // store the current length^2 of the vector
      len += (dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn)*gauss_s2*gauss_s2;

      j += 9;
    }
    i += 9;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;

  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

/** Get the modified gauge 128 upright SURF descriptor vector of the provided Ipoint   */
/** See Agrawal et al. (ECCV08) for more details  									   */
/** Descriptor Dimension 128														   */
void Surf::getMGaugeUprightExtendedDescriptor_128()
{
  int y = 0, x = 0, sample_x = 0, sample_y = 0, count=0, dsize = 0;
  int i = 0, ix = 0, j = 0, jx = 0, xs = 0, ys = 0, sample_step = 0, pattern_size = 0;
  float scale = 0.0, *desc;
  float dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0;
  float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0;
  float gauss_s1 = 0.0, gauss_s2 = 0.0, lww = 0.0, lvv = 0.0, modg = 0.0;
  float rx = 0.0, ry = 0.0, rxx = 0.0, rxy = 0.0, ryy = 0.0, len = 0.0;
  float cx = -0.5f, cy = 0.f; //Subregion centers for the 4x4 gaussian weighting

  Ipoint *ipt = &ipts[index];
  scale = ipt->scale;
  x = fRound(ipt->x);
  y = fRound(ipt->y);
  desc = ipt->descriptor;
  dsize = ipt->size;
  sample_step = 5;
  pattern_size = 12;
  
  i = -8;

  // Calculate descriptor for this interest point
  // Area of size 24 s x 24 s
  while(i < pattern_size)
  {
    j = -8;
    i = i-4;

    cx += 1.f;
    cy = -0.5f;

    while(j < pattern_size)
    {
      dxp=dyp=mdxp=mdyp=0.0f;
      dxn=dyn=mdxn=mdyn=0.0f;

      cy += 1.f;

      j = j - 4;

      ix = i + sample_step;
      jx = j + sample_step;

      xs = fRound(x + ix*scale);
      ys = fRound(y + jx*scale);

      for (int k = i; k < i + 9; ++k)
      {
        for (int l = j; l < j + 9; ++l)
        {
          // Get coords of sample point on the rotated axis
          sample_x = fRound(x + k*scale);
          sample_y = fRound(y + l*scale);

          // Get the gaussian weighted x and y responses
          gauss_s1 = gaussian(xs-sample_x,ys-sample_y,2.5f*scale);

          rx = haarX(sample_y, sample_x, 2*fRound(scale));
          ry = haarY(sample_y, sample_x, 2*fRound(scale));

		  modg = pow(rx,2) + pow(ry,2);
		  
          if( modg != 0.0 )
          {    
			   rxx =  haarXX(sample_y,sample_x,2*fRound(scale));
			   rxy =  haarXY(sample_y,sample_x,2*fRound(scale));
			   ryy =  haarYY(sample_y,sample_x,2*fRound(scale));

			   // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2)
			   lww = gauss_s1*(pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg);

			   // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2)
			   lvv = gauss_s1*(-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg);
          }
          else
          {
                lww = 0.0;
                lvv = 0.0;
          }

          if( lvv >= 0.0 )
          {
              dxp += lww;
              mdxp += fabs(lww);
          }
          else
          {
              dxn += lww;
              mdxn += fabs(lww);
          }

          if( lww >= 0.0 )
          {
              dyp += lvv;
              mdyp += fabs(lvv);
          }
          else
          {
              dyn += lvv;
              mdyn += fabs(lvv);
          }
        }
      }

      // Add the values to the descriptor vector
      gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f);

      desc[count++] = dxp*gauss_s2;
      desc[count++] = mdxp*gauss_s2;
      desc[count++] = dxn*gauss_s2;
      desc[count++] = mdxn*gauss_s2;

      desc[count++] = dyp*gauss_s2;
      desc[count++] = mdyp*gauss_s2;
      desc[count++] = dyn*gauss_s2;
      desc[count++] = mdyn*gauss_s2;

      // store the current length^2 of the vector
      len += (dxp*dxp + mdxp*mdxp + dxn*dxn + mdxn*mdxn + dyp*dyp + mdyp*mdyp + dyn*dyn + mdyn*mdyn)*gauss_s2*gauss_s2;

      j += 9;
    }
    i += 9;
  }

  // Convert to Unit Vector
  len = sqrt(len);

  for(int i = 0; i < dsize; i++)
    desc[i] /= len;
	
  if( USE_CLIPPING_NORMALIZATION == true )
  {
	  clippingDescriptor(CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO);
  }
}

//**************************************************************************************
//**************************************************************************************

//! Calculate the value of the 2d gaussian at x,y
inline float Surf::gaussian(int x, int y, float sig)
{
  return exp( -(x*x+y*y)/(2.0f*sig*sig));
}

//**************************************************************************************
//**************************************************************************************

//! Calculate the value of the 2d gaussian at x,y
inline float Surf::gaussian(float x, float y, float sig)
{
  return 1.0f/(2.0f*pi*sig*sig) * exp( -(x*x+y*y)/(2.0f*sig*sig));
}

//**************************************************************************************
//**************************************************************************************

//! Calculate Haar wavelet responses in x direction
inline float Surf::haarX(int row, int column, int s)
{
  return BoxIntegral(img, row-s/2, column, s, s/2) 
    -1 * BoxIntegral(img, row-s/2, column-s/2, s, s/2);
}

//**************************************************************************************
//**************************************************************************************

//! Calculate Haar wavelet responses in y direction
inline float Surf::haarY(int row, int column, int s)
{
  return BoxIntegral(img, row, column-s/2, s/2, s) 
    -1 * BoxIntegral(img, row-s/2, column-s/2, s/2, s);
}

//**************************************************************************************
//**************************************************************************************

//! Calculate Haar wavelet responses in xy direction
inline float Surf::haarXY(int row, int column, int s)
{
  s = s+1;

  return BoxIntegral(img, row-s,column-s,s,s) + BoxIntegral(img,row+1,column+1,s,s)
        -BoxIntegral(img, row+1,column-s,s,s) - BoxIntegral(img,row-s,column+1,s,s);
}

//**************************************************************************************
//**************************************************************************************

//! Calculate Haar wavelet responses in xx direction
inline float Surf::haarXX(int row, int column, int s)
{
  s = s+1;
  int r = row;
  int c = column;
  int l = s;
  int w = 3*s;
  int b = w/2;

  return BoxIntegral(img, r - l + 1, c - b, 2*l - 1, w)
              - BoxIntegral(img, r - l + 1, c - l / 2, 2*l - 1, l)*3;	
}

//**************************************************************************************
//**************************************************************************************

//! Calculate Haar wavelet responses in yy direction
inline float Surf::haarYY(int row, int column, int s)
{
  s = s+1;
  int r = row;
  int c = column;
  int l = s;
  int w = 3*s;
  int b = w/2;

  return BoxIntegral(img, r - b, c - l + 1, w, 2*l - 1)
              - BoxIntegral(img, r - l / 2, c - l + 1, l, 2*l - 1)*3;
}

//**************************************************************************************
//**************************************************************************************

//! Get the angle from the +ve x-axis of the vector given by (X Y)
float Surf::getAngle(float X, float Y)
{
  if(X >= 0 && Y >= 0)
    return atan(Y/X);

  if(X < 0 && Y >= 0)
    return pi - atan(-Y/X);

  if(X < 0 && Y < 0)
    return pi + atan(Y/X);

  if(X >= 0 && Y < 0)
    return 2*pi - atan(-Y/X);

  return 0;

}

//**************************************************************************************
//**************************************************************************************

/** This function performs descriptor clipping											*/
void Surf::clippingDescriptor(int niter, float ratio)
{
	Ipoint *ipt = &ipts[index];
	float *desc =  ipt->descriptor;
	int dsize = ipt->size;
	float cratio = ratio / sqrt(dsize);
	float len = 0.0;

	for( int i = 0; i < niter; i++ )
	{
		len = 0.0;
		for( int j = 0; j < dsize; j++ )
		{
			if( desc[j] > cratio )
			{
				desc[j] = cratio;
			}
			else if( desc[j] < -cratio )
			{
				desc[j] = -cratio;
			}
			len += desc[j]*desc[j];
		}
		
		// Normalize again
		len = sqrt(len);

		for( int j = 0; j < dsize; j++ )
		desc[j] = desc[j] / len;
	}
}