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;
static public int OctavesCount = 2;
static public int scalesPerOctave = 4;
static double minDoG = 0.01;
static public double norm;
static double edgeRatio = 10;
public static double[][,] BuildGaussianBlurPyramid(Bitmap image)
{ int ind=1;
int imgCount = scalesPerOctave * OctavesCount;
double deviation = initDev;
double[][,] imgMatrix=new double[imgCount][,];
imgMatrix[0]=new double[image.Width,image.Height];
imgMatrix[0]=ImageProcessing.BitmapToGrayscaleMatrix(image);
for (int i = 0; i < OctavesCount; i++){
for (int j = 1; j < scalesPerOctave; j++){
imgMatrix[ind]=ImageProcessing.GaussianBlurBWMat(imgMatrix[ind-1], deviation);
deviation *= 2;
ind++;
}
if (ind<imgCount-1)
imgMatrix[ind] = ImageProcessing.SubSample(imgMatrix[ind-2]);
deviation /= 3.2;
ind++;
}
return imgMatrix;
}
public static double[][,] BuildDoGPyramid(double[][,] gaussianPyr)
{ int ind = 0;
double[][,] resDoG = new double[(scalesPerOctave - 1) * OctavesCount][,];
for (int i = 0; i < gaussianPyr.GetLength(0) - 1; i += 1){
if ((i + 1) % (scalesPerOctave) == 0){
continue;
}
resDoG[ind] = ImageProcessing.CalculateDoGMatrix(gaussianPyr[i], gaussianPyr[i + 1]);
ind++;
}
return resDoG;
}
public static ArrayList FindKeypoints(double[][,] imgMatrixPyr)
{ ArrayList KpL = new ArrayList();
double[][][,] imgMatrix = new double[2][][,];
norm = 3.0 / (initDev * 2);
int DoGCount = (scalesPerOctave - 1) * OctavesCount;
imgMatrix[0] = new double[3][,];
imgMatrix[1] = new double[3][,];
int pyrL = 0;
int pyrIm = 0;
for (int i = 0; i < DoGCount; i++){
if (i == 3){
pyrL++;
pyrIm = 0;
}
imgMatrix[pyrL][pyrIm] = imgMatrixPyr[i];
pyrIm++;
}
int min = 1, max = 1;
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];
if (Math.Abs(point) <= minDoG)
continue;
for (int i = 0; i < imgMatrix[l].GetLength(0); i++){
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;
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){
Keypoint kp = new Keypoint();
kp.x = x;
kp.y = y;
kp.value = point;
kp.layer = l;
KpL.Add(kp); }
min = 1; max = 1;
}
}
return KpL;
}
public static ArrayList KeypointLocalization(ArrayList initKpList, double[][,] DoGMatrix, int edgeRatio)
{ Keypoint initKp;
ArrayList retKp = new ArrayList();
double Dxx, Dyy, Dxy, Trace, Det;
double r = ((edgeRatio + 1) * (edgeRatio + 1)) / 2;
for (int i = 0; i < initKpList.Count; i++){
initKp = (Keypoint)initKpList[i];
int l = initKp.layer == 0 ? 1 : 4;
if (initKp.x <= 0 || initKp.y <= 0 ||
initKp.x >= DoGMatrix[l].GetLength(0) - 1
|| initKp.y >= DoGMatrix[l].GetLength(0) - 1)
continue;
Dxx = DoGMatrix[l][initKp.x + 1, initKp.y]
+ DoGMatrix[l][initKp.x - 1, initKp.y]
- 2.0 * DoGMatrix[l][initKp.x, initKp.y];
Dyy = DoGMatrix[l][initKp.x, initKp.y + 1]
+ 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]
- 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;
Det = Dxx * Dyy - (Dxy * Dxy);
if ((Trace * Trace) / Det > r)
continue;
else
retKp.Add(initKp);
}
return retKp;
}
public static ArrayList AssignOrientations(ArrayList finKpList, double[][,] blurredMatrix)
{ Keypoint finKp;
ArrayList orKp=new ArrayList();
double binVal = 2 * Math.PI / 36;
for (int i = 0; i < finKpList.Count; i++){
finKp = (Keypoint)finKpList[i];
finKp.scale = initDev*(finKp.layer + 1);
double[] histBin = new double[36];
int kpReg=(int)(1.5*finKp.scale/2+0.5);
double SigmaFact=2*2.25*finKp.scale*finKp.scale;
int layer=finKp.layer==0?1:5;
int regXSt=finKp.x-kpReg<1?1:finKp.x-kpReg;
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;
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++){
double gaussianWight = Math.Exp(-(Math.Pow(x - finKp.x, 2) + Math.Pow(y - finKp.y, 2) / SigmaFact));
double portInHist = (CalMagnitudeAndDirection(x,y, blurredMatrix[layer]).direction + Math.PI) / (2 * Math.PI);
int ind = (int)(portInHist * 36 == 36 ? 0 : portInHist * 36);
histBin[ind]+=CalMagnitudeAndDirection(x,y,blurredMatrix[layer]).magnitude*gaussianWight;
}
double max = 0.0;
int maxBin = 0;
for (int bi = 0; bi < histBin.GetLength(0); ++bi){
if (histBin[bi] > max){
max = histBin[bi];
maxBin = bi; }
}
max = maxBin * binVal - Math.PI;
max=max<-Math.PI?2*Math.PI+max:max;
max = max > Math.PI ? max - 2 * Math.PI : max;
finKp.orientation = max;
orKp.Add(finKp);
}
return orKp;
}
static KpOrientPar CalMagnitudeAndDirection(int x, int y, double[,] blurredMatrix)
{ KpOrientPar kpO=new KpOrientPar();
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));
kpO.direction = Math.Atan2(blurredMatrix[x, y + 1] - blurredMatrix[x, y - 1],
blurredMatrix[x + 1, y] - blurredMatrix[x - 1, y]);
return kpO;
}
public static ArrayList CreateDescriptor(ArrayList finKpList, double[][,] blurredMatrix)
{ int kpReg = 16;
Keypoint finKp;
ArrayList descKp=new ArrayList();
double SigmaFact = 2*(kpReg / 2) * (kpReg / 2);
for (int kpi = 0; kpi < finKpList.Count; kpi++){
finKp = (Keypoint)finKpList[kpi];
int regXSt = finKp.x - kpReg < 1 ? 1 : finKp.x - kpReg;
int regYSt = finKp.y - kpReg < 1 ? 1 : finKp.y - kpReg;
int regXEnd = regXSt + 4;
int regYEnd = regXSt + 4;
double[,][] desc = new double[4, 4][];
for (int jj = 0; jj < 4; jj++)
for (int ii = 0; ii < 4; ii++ )
desc[ii,jj] = new double[8];
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++){
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);
double xRot = Math.Cos(-finKp.orientation) * (x - finKp.x) -
Math.Sin(-finKp.orientation) * (y - finKp.y);
double gaussianWight = Math.Exp(-(Math.Pow(xRot, 2) + Math.Pow(yRot, 2) / SigmaFact));
double direction = CalMagnitudeAndDirection(x, y, blurredMatrix[layer]).direction - finKp.orientation;
if (direction < 0)
direction += 2.0 * Math.PI;
else if (direction > 2.0 * Math.PI)
direction -= 2.0 * Math.PI;
double portInHist = direction / (2 * Math.PI);
int ind = (int)(portInHist * 8 == 8 ? 0 : portInHist * 8);
desc[i, j][ind] += CalMagnitudeAndDirection(x, y, blurredMatrix[layer]).magnitude * gaussianWight;
}
finKp.descriptor=desc;
descKp.Add(finKp);
}
return descKp;
}
}
class KpOrientPar
{
public double magnitude;
public double direction;
}
|