C#

/* SIFT.cs
 * The SIFT algorithm class
 * where all the algorithm specific methods
 * are implemented 
 */
using System;
using System.Collections.Generic;
using System.Text;
using System.Drawing;
using System.Collections;
using System.Windows.Forms;
class SIFT
{   static double initDev = 1.6;                                        // initial standard deviation variable
    static public int OctavesCount = 2;                                 // Octaves count for pyramid
    static public int scalesPerOctave = 4;                              // scales(blurred images) count per octave
    static double minDoG = 0.01;                                        // minimal threshold of keypoint value
    static public double norm;                                          // normalization factor for kp coordinates
    static double edgeRatio = 10;                                       // minimal edge ratio for kp localization
    public static double[][,] BuildGaussianBlurPyramid(Bitmap image)    // method to build the Gaussian pyramid with bitmap input
    {   int ind=1;                                                      // index variable to help for image indexing in pyramid
        int imgCount = scalesPerOctave * OctavesCount;                  // calculate total image count 
        double deviation = initDev;                                     // local initial deviation variable
        double[][,] imgMatrix=new double[imgCount][,];                  // array to store values of each image matrix in the pyramid 
        imgMatrix[0]=new double[image.Width,image.Height];              // image matrix initial dimension  
        imgMatrix[0]=ImageProcessing.BitmapToGrayscaleMatrix(image);    // create grayscale image matrix from bitmap coloured image
        for (int i = 0; i < OctavesCount; i++){             
         for (int j = 1; j < scalesPerOctave; j++){                     // Build Gaussian blur with given image mtrx and deviation 
          imgMatrix[ind]=ImageProcessing.GaussianBlurBWMat(imgMatrix[ind-1], deviation);
          deviation *= 2;                                               // double the previous deviation
          ind++;                                                        // increment current image index
                                                   }
          if (ind<imgCount-1)                                           // check for current image index and build next octave
          imgMatrix[ind] = ImageProcessing.SubSample(imgMatrix[ind-2]); // down-sample the image by factor of 2 for the next octave 
          deviation /= 3.2;                                             // divide the deviation to start from twice the initial
          ind++;                                                        // increment the current image index
              
                                               }
        return imgMatrix;                                               // return the resulting image matrix pyramid array
      }                                                                 // method to build the DoG pyramid
      public static double[][,] BuildDoGPyramid(double[][,] gaussianPyr) 
      {   int ind = 0;                                                  // index value for controlling the image indxes 
          double[][,] resDoG = new double[(scalesPerOctave - 1) * OctavesCount][,];
          for (int i = 0; i < gaussianPyr.GetLength(0) - 1; i += 1){
              if ((i + 1) % (scalesPerOctave) == 0){                    // go to next octave
                 continue;
                                                   }                    // calculate the DoG matrix
              resDoG[ind] = ImageProcessing.CalculateDoGMatrix(gaussianPyr[i], gaussianPyr[i + 1]);
              ind++;                                                    // increment the current image index
          }

          return resDoG;                                                // return resulting DoG matrix pyramid
      }
      public static ArrayList FindKeypoints(double[][,] imgMatrixPyr)   // method to find keypoints
      {   ArrayList KpL = new ArrayList();                              // create arraylist to store keypoints
          double[][][,] imgMatrix = new double[2][][,];                 // create array to store three level pyr img matrix
          norm = 3.0 / (initDev * 2);                                   // calculate the normalization factor
          int DoGCount = (scalesPerOctave - 1) * OctavesCount;          // count the total number of DoG images 
          imgMatrix[0] = new double[3][,];                              // array to store first octave image matrix
          imgMatrix[1] = new double[3][,];                              // array to store second octave image matrix
          int pyrL = 0;                                                 // layer index 
          int pyrIm = 0;                                                // image index in layer
          for (int i = 0; i < DoGCount; i++){                            
              if (i == 3){
                  pyrL++;                                               // increment layer index
                  pyrIm = 0;                                            // assign zero to image index to start next dimension 
                         } 
              imgMatrix[pyrL][pyrIm] = imgMatrixPyr[i];                 // copy DoG pyr. image matrix from 1 dim array to 2 dim.
              pyrIm++;                                                  // increment image index
                                            }
          int min = 1, max = 1;                                         // min and max initialization
          for (int l = 0; l < imgMatrix.GetLength(0); l++){
           for (int y = 1; y < imgMatrix[l][1].GetLength(1) - 1; y++)
            for (int x = 1; x < imgMatrix[l][1].GetLength(0) - 1; x++){
             double point = imgMatrix[l][1][x, y];                      // pick one point of the central DoG image matrix
             if (Math.Abs(point) <= minDoG)                             // check for the keypoint threshold
             continue;
             for (int i = 0; i < imgMatrix[l].GetLength(0); i++){
                                                                        // check the point for maximum 
             if (max == 1 & imgMatrix[l][i][x - 1, y - 1] < point & imgMatrix[l][i][x - 1, y] < point & 
             imgMatrix[l][i][x, y - 1] < point & imgMatrix[l][i][x, y] <= point & imgMatrix[l][i][x + 1, y + 1] < point & 
             imgMatrix[l][i][x + 1, y] < point & imgMatrix[l][i][x, y + 1] < point & imgMatrix[l][i][x - 1, y + 1] < point & 
             imgMatrix[l][i][x + 1, y - 1] < point){
             max = 1;
                                                   }
             else max = 0;     
                                                                        // check the point for minimum
             if (min == 1 & imgMatrix[l][i][x - 1, y - 1] > point & imgMatrix[l][i][x - 1, y] > point & 
             imgMatrix[l][i][x, y - 1] > point & imgMatrix[l][i][x, y] >= point & imgMatrix[l][i][x + 1, y + 1] > point & 
             imgMatrix[l][i][x + 1, y] > point & imgMatrix[l][i][x, y + 1] > point & imgMatrix[l][i][x - 1, y + 1] > point & 
             imgMatrix[l][i][x + 1, y - 1] > point){
             min = 1;                              }
             else min = 0;
                                                   }
             if (max == 1 || min == 1){                                 // check if the point is max or min
                Keypoint kp = new Keypoint();                           // create keypoint object to 
                kp.x = x;                                               // assign x coordinate of extreme point
                kp.y = y;                                               // assign y coordinate of extreme point
                kp.value = point;                                       // assign extreme point value
                kp.layer = l;                                           // assign layer of extreme point
                KpL.Add(kp);          }                                 // add keypoint to arraylist
             min = 1; max = 1;                                          // reinitialize min and max values  
           }
          }
          return KpL;                                                   // return keypoint list
      }
                                                                        // method to localize keypoints for the given edge ratio
    public static ArrayList KeypointLocalization(ArrayList initKpList, double[][,] DoGMatrix, int edgeRatio)
    {   Keypoint initKp;                                                // Keypoint object to store initial keypoints
        ArrayList retKp = new ArrayList();                              // create Arraylist to store refined keypoints 
        double Dxx, Dyy, Dxy, Trace, Det;                               // variables to store elements, trace and Det. of Hessian
        double r = ((edgeRatio + 1) * (edgeRatio + 1)) / 2;             // calculate the threshold from the edgratio
        for (int i = 0; i < initKpList.Count; i++){
          initKp = (Keypoint)initKpList[i];                             // cast the value from Arraylist object to the Keypoint
          int l = initKp.layer == 0 ? 1 : 4;                            // identify the current layer
          if (initKp.x <= 0 || initKp.y <= 0 ||                         // drop keypoints which are on the borders
              initKp.x >= DoGMatrix[l].GetLength(0) - 1 
              || initKp.y >= DoGMatrix[l].GetLength(0) - 1)
          continue;
          Dxx = DoGMatrix[l][initKp.x + 1, initKp.y]                    // calculate the x,x derivative by pixel differences
          + DoGMatrix[l][initKp.x - 1, initKp.y] 
          - 2.0 * DoGMatrix[l][initKp.x, initKp.y];
          Dyy = DoGMatrix[l][initKp.x, initKp.y + 1]                    //calculate the y,y derivative by pixel differences
          + DoGMatrix[l][initKp.x, initKp.y - 1] 
          - 2.0 * DoGMatrix[l][initKp.x, initKp.y];
          Dxy = 0.25 * ((DoGMatrix[l][initKp.x + 1, initKp.y + 1]       // calculate the x,y derivative by pixel differences
          - DoGMatrix[l][initKp.x + 1, initKp.y - 1]) -
          (DoGMatrix[l][initKp.x - 1, initKp.y + 1] 
          - DoGMatrix[l][initKp.x - 1, initKp.y - 1]));
          Trace = Dxx + Dyy;                                            // calculate the trace of Hessian
          Det = Dxx * Dyy - (Dxy * Dxy);                                // calculate the determinant of Hessian
          if ((Trace * Trace) / Det > r)                                // drop the keypoints which contradict the given ratio
          continue;
          else
          retKp.Add(initKp);                                            // add the remaining keypoints
                                                  }
       return retKp;                                                    // return the resulting Arraylist with localized keypoints
      }
                                                                        // method for orientation assignment
      public static ArrayList AssignOrientations(ArrayList finKpList, double[][,] blurredMatrix)
      {   Keypoint finKp;                                               // keypoint variable to store localized keyoints
          ArrayList orKp=new ArrayList();                               // create arraylist to store keypoints with orientations 
          double binVal = 2 * Math.PI / 36;                             // calculate the angle value of one histogram bin
          for (int i = 0; i < finKpList.Count; i++){                  
           finKp = (Keypoint)finKpList[i];                              // cast from array list object to keypoint object
           finKp.scale = initDev*(finKp.layer + 1);                     // caclulate the keypoint scale
           double[] histBin = new double[36];                           // create the double array as a colour histogram
           int kpReg=(int)(1.5*finKp.scale/2+0.5);                      // calculate the radius of the region around the keypoint 
           double SigmaFact=2*2.25*finKp.scale*finKp.scale;             // caclutate the sigma factor for gaussian weight
           int layer=finKp.layer==0?1:5;                                // identify the current keypoint layer
           int regXSt=finKp.x-kpReg<1?1:finKp.x-kpReg;                  // calculate the x coordinate of the region start and end 
           int regXEnd = finKp.x + kpReg > blurredMatrix[layer].GetLength(0) ? blurredMatrix[layer].GetLength(0) : finKp.x + kpReg;
           int regYSt = finKp.y - kpReg < 1 ? 1 : finKp.y - kpReg;      // calculate the x coordinate of the region start and end 
           int regYEnd = finKp.y + kpReg > blurredMatrix[layer].GetLength(1) ? blurredMatrix[layer].GetLength(1) : finKp.y + kpReg;
           for (int y=regYSt; y<regYEnd; y++)
            for (int x=regXSt; x<regXEnd; x++){                         // calculate the gaussian wight 
             double gaussianWight = Math.Exp(-(Math.Pow(x - finKp.x, 2) + Math.Pow(y - finKp.y, 2) / SigmaFact));
                                                                        // calculate the portion of the direction angle in histogram
             double portInHist = (CalMagnitudeAndDirection(x,y, blurredMatrix[layer]).direction + Math.PI) / (2 * Math.PI);
                                                                        // identify the bin index for current keypoint direction 
             int ind = (int)(portInHist * 36 == 36 ? 0 : portInHist * 36); 
                                                                        // add magnitude weightened by gaussian to the respective bin
             histBin[ind]+=CalMagnitudeAndDirection(x,y,blurredMatrix[layer]).magnitude*gaussianWight;
                                              }
            double max = 0.0;                                           // initialize the maximal bin value
            int maxBin = 0;                                             // initialize the maximal bin index
            for (int bi = 0; bi < histBin.GetLength(0); ++bi){
             if (histBin[bi] > max){                                    // check if the next bin is bigger
              max = histBin[bi];                                        // store new bin as the maximal
              maxBin = bi;         }                                    // mark the new bin index as the index of the maximal bin  
                                                              }
            max = maxBin * binVal - Math.PI;                            // change angle range from [0, 2PI] to [-PI,PI]
            max=max<-Math.PI?2*Math.PI+max:max;                         // check if angle is outside [-PI,PI]
            max = max > Math.PI ? max - 2 * Math.PI : max;              // check if angle is outside [-PI,PI]
            finKp.orientation = max;                                    // assign maximal bin direction as keypoint orientation
            orKp.Add(finKp);                                            // add keypoint with assign orientation
                                                    }
          return orKp;                                                  // return keypoint list of keypoints with their orientations
        }
                                                                        // method to calculate keypoint magnitude and direction
        static KpOrientPar CalMagnitudeAndDirection(int x, int y, double[,] blurredMatrix)
        {   KpOrientPar kpO=new KpOrientPar();                          // create orientation object for magnitude and direction pair
                                                                        // calculate the keypoint magnitude by pixel differences
            kpO.magnitude = Math.Sqrt (Math.Pow (blurredMatrix[x + 1, y] - blurredMatrix[x - 1, y], 2.0) +
			Math.Pow (blurredMatrix[x, y + 1] -	blurredMatrix[x, y - 1], 2.0));
                                                                        // calculate the keypoint direction by pixel differences 
            kpO.direction = Math.Atan2(blurredMatrix[x, y + 1] - blurredMatrix[x, y - 1],
			blurredMatrix[x + 1, y] - blurredMatrix[x - 1, y]);
            return kpO;                                                 // return orientation object with magnitude and direction
        }
        public static ArrayList CreateDescriptor(ArrayList finKpList, double[][,] blurredMatrix)
        {   int kpReg = 16;                                             // region size around the keypoint
            Keypoint finKp;                                             // keypoint object to store input keypoint
            ArrayList descKp=new ArrayList();                           // arraylist keypoint to store descriptor assigned keypoints 
            double SigmaFact = 2*(kpReg / 2) * (kpReg / 2);             // calculate the sigma factor for gaussian weight
            for (int kpi = 0; kpi < finKpList.Count; kpi++){    
             finKp = (Keypoint)finKpList[kpi];                          // cast from arraylist object to keypoint object
             int regXSt = finKp.x - kpReg < 1 ? 1 : finKp.x - kpReg;    // calculate x coordinate of keypoint region start
             int regYSt = finKp.y - kpReg < 1 ? 1 : finKp.y - kpReg;    // calculate y coordinate of keypoint region start
             int regXEnd = regXSt + 4;                                  // calculate x coordinate of keypoint region end
             int regYEnd = regXSt + 4;                                  // calculate y coordinate of keypoint region end
             double[,][] desc = new double[4, 4][];                     // create double array to store descriptor
             for (int jj = 0; jj < 4; jj++)
              for (int ii = 0; ii < 4; ii++ )
               desc[ii,jj] = new double[8];                             // initialize double array of descriptor
               int layer = finKp.layer == 0 ? 1 : 5;
               for (int i = 0; i <= 3; i++)
                for (int j =0 ; j <= 3; j++)
                 for (int xi = regXSt; xi < regXEnd; xi++)
                  for (int yi = regYSt; yi < regYEnd; yi++){            // calculate current value of x and y coordinates
                   int x =xi+j*4 >= blurredMatrix[layer].GetLength(0)-1 ? blurredMatrix[layer].GetLength(0)-2 : xi+j*4;
                   int y = yi+i*4 >= blurredMatrix[layer].GetLength(1)-1 ? blurredMatrix[layer].GetLength(1)-2 : yi+i*4;
                   double yRot = Math.Sin(-finKp.orientation) * (x - finKp.x) +
                   Math.Cos(-finKp.orientation) * (y - finKp.y);        // rotate x coordinate by keypoint orientation value
                   double xRot = Math.Cos(-finKp.orientation) * (x - finKp.x) -
                   Math.Sin(-finKp.orientation) * (y - finKp.y);        // rotate y coordinate by keypoint orientation value
                                                                        // calculate gaussian weight 
                   double gaussianWight = Math.Exp(-(Math.Pow(xRot, 2) + Math.Pow(yRot, 2) / SigmaFact));
                                                                        // calculate direction and rotate by keypoint orientation
                   double direction = CalMagnitudeAndDirection(x, y, blurredMatrix[layer]).direction - finKp.orientation;
                   if (direction < 0)                                   // check if direction degree is out of [0,2PI] range
                    direction += 2.0 * Math.PI;
                   else if (direction > 2.0 * Math.PI)
                         direction -= 2.0 * Math.PI;
                   double portInHist = direction / (2 * Math.PI);       // calculate the protion of the direction in the histogram
                                                                        // idetify the index of orientation bin
                   int ind = (int)(portInHist * 8 == 8 ? 0 : portInHist * 8);
                                                                        // add magn. weighted by gaussian to the resp.direction bin
                   desc[i, j][ind] += CalMagnitudeAndDirection(x, y, blurredMatrix[layer]).magnitude * gaussianWight;
                                                          }
                finKp.descriptor=desc;                                  // assign descriptor to keypoint object
                descKp.Add(finKp);                                      // add keypoint to the resulting arraylist object
            }
            return descKp;                                              // return arraylist of keypoints with descriptors
        }
}
class KpOrientPar                                                       // helper class to store magnitude and direction
{
    public double magnitude;                                            // double variable to store keypoint magnitude
    public double direction;                                            // double variable to store keypoint direction
}