From f618466c25d43f3bae9e40920273bf77de1e1149 Mon Sep 17 00:00:00 2001 From: leochanj105 Date: Mon, 19 Oct 2020 23:09:30 -0400 Subject: initial sd-vbs initial sd-vbs add sd-vbs sd-vbs --- SD-VBS/benchmarks/sift/src/c/diffss.c | 53 +++ SD-VBS/benchmarks/sift/src/c/doubleSize.c | 55 +++ .../benchmarks/sift/src/c/filterBoundaryPoints.c | 61 +++ SD-VBS/benchmarks/sift/src/c/gaussianss.c | 185 ++++++++ SD-VBS/benchmarks/sift/src/c/halveSize.c | 35 ++ SD-VBS/benchmarks/sift/src/c/imsmooth.c | 104 ++++ SD-VBS/benchmarks/sift/src/c/script_sift.c | 83 ++++ SD-VBS/benchmarks/sift/src/c/sift.c | 314 ++++++++++++ SD-VBS/benchmarks/sift/src/c/sift.h | 27 ++ SD-VBS/benchmarks/sift/src/c/siftlocalmax.c | 247 ++++++++++ SD-VBS/benchmarks/sift/src/c/siftrefinemx.c | 196 ++++++++ SD-VBS/benchmarks/sift/src/matlab/diffss.m | 72 +++ SD-VBS/benchmarks/sift/src/matlab/gaussianss.m | 216 +++++++++ SD-VBS/benchmarks/sift/src/matlab/imreadbw.m | 52 ++ SD-VBS/benchmarks/sift/src/matlab/imsmooth.c | 115 +++++ SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexa64 | Bin 0 -> 9570 bytes SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexglx | Bin 0 -> 7450 bytes SD-VBS/benchmarks/sift/src/matlab/imsmooth_.c | 103 ++++ SD-VBS/benchmarks/sift/src/matlab/mexutils.c | 98 ++++ SD-VBS/benchmarks/sift/src/matlab/plotmatches.m | 124 +++++ .../sift/src/matlab/plotsiftdescriptor.m | 169 +++++++ SD-VBS/benchmarks/sift/src/matlab/plotsiftframe.m | 119 +++++ SD-VBS/benchmarks/sift/src/matlab/plotss.m | 68 +++ .../sift/src/matlab/script_run_profile.m | 23 + SD-VBS/benchmarks/sift/src/matlab/sift.m | 110 +++++ SD-VBS/benchmarks/sift/src/matlab/sift_compile.m | 60 +++ SD-VBS/benchmarks/sift/src/matlab/sift_demo2.m | 110 +++++ SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.css | 159 +++++++ SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.m | 50 ++ SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.pl | 296 ++++++++++++ SD-VBS/benchmarks/sift/src/matlab/sift_overview.m | 33 ++ SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.c | 524 +++++++++++++++++++++ SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.m | 79 ++++ .../sift/src/matlab/siftdescriptor.mexa64 | Bin 0 -> 15700 bytes .../sift/src/matlab/siftdescriptor.mexglx | Bin 0 -> 13332 bytes SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.c | 291 ++++++++++++ SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.m | 33 ++ .../benchmarks/sift/src/matlab/siftlocalmax.mexa64 | Bin 0 -> 10836 bytes .../benchmarks/sift/src/matlab/siftlocalmax.mexglx | Bin 0 -> 8518 bytes SD-VBS/benchmarks/sift/src/matlab/siftmatch.c | 250 ++++++++++ SD-VBS/benchmarks/sift/src/matlab/siftmatch.m | 60 +++ SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexa64 | Bin 0 -> 11877 bytes SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexglx | Bin 0 -> 9607 bytes SD-VBS/benchmarks/sift/src/matlab/siftormx.c | 251 ++++++++++ SD-VBS/benchmarks/sift/src/matlab/siftormx.mexa64 | Bin 0 -> 12860 bytes SD-VBS/benchmarks/sift/src/matlab/siftormx.mexglx | Bin 0 -> 10524 bytes SD-VBS/benchmarks/sift/src/matlab/siftread.m | 101 ++++ SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.c | 342 ++++++++++++++ SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.m | 61 +++ .../benchmarks/sift/src/matlab/siftrefinemx.mexa64 | Bin 0 -> 12392 bytes .../benchmarks/sift/src/matlab/siftrefinemx.mexglx | Bin 0 -> 9007 bytes SD-VBS/benchmarks/sift/src/matlab/tightsubplot.m | 151 ++++++ 52 files changed, 5480 insertions(+) create mode 100644 SD-VBS/benchmarks/sift/src/c/diffss.c create mode 100644 SD-VBS/benchmarks/sift/src/c/doubleSize.c create mode 100644 SD-VBS/benchmarks/sift/src/c/filterBoundaryPoints.c create mode 100644 SD-VBS/benchmarks/sift/src/c/gaussianss.c create mode 100644 SD-VBS/benchmarks/sift/src/c/halveSize.c create mode 100644 SD-VBS/benchmarks/sift/src/c/imsmooth.c create mode 100644 SD-VBS/benchmarks/sift/src/c/script_sift.c create mode 100644 SD-VBS/benchmarks/sift/src/c/sift.c create mode 100644 SD-VBS/benchmarks/sift/src/c/sift.h create mode 100644 SD-VBS/benchmarks/sift/src/c/siftlocalmax.c create mode 100644 SD-VBS/benchmarks/sift/src/c/siftrefinemx.c create mode 100644 SD-VBS/benchmarks/sift/src/matlab/diffss.m create mode 100644 SD-VBS/benchmarks/sift/src/matlab/gaussianss.m create mode 100644 SD-VBS/benchmarks/sift/src/matlab/imreadbw.m create mode 100644 SD-VBS/benchmarks/sift/src/matlab/imsmooth.c create mode 100755 SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexa64 create mode 100755 SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexglx create mode 100644 SD-VBS/benchmarks/sift/src/matlab/imsmooth_.c create mode 100644 SD-VBS/benchmarks/sift/src/matlab/mexutils.c create mode 100644 SD-VBS/benchmarks/sift/src/matlab/plotmatches.m create mode 100644 SD-VBS/benchmarks/sift/src/matlab/plotsiftdescriptor.m create mode 100644 SD-VBS/benchmarks/sift/src/matlab/plotsiftframe.m create mode 100644 SD-VBS/benchmarks/sift/src/matlab/plotss.m create mode 100644 SD-VBS/benchmarks/sift/src/matlab/script_run_profile.m create mode 100644 SD-VBS/benchmarks/sift/src/matlab/sift.m create mode 100644 SD-VBS/benchmarks/sift/src/matlab/sift_compile.m create mode 100644 SD-VBS/benchmarks/sift/src/matlab/sift_demo2.m create mode 100644 SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.css create mode 100644 SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.m create mode 100755 SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.pl create mode 100644 SD-VBS/benchmarks/sift/src/matlab/sift_overview.m create mode 100644 SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.c create mode 100644 SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.m create mode 100755 SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.mexa64 create mode 100755 SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.mexglx create mode 100644 SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.c create mode 100644 SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.m create mode 100755 SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.mexa64 create mode 100755 SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.mexglx create mode 100644 SD-VBS/benchmarks/sift/src/matlab/siftmatch.c create mode 100644 SD-VBS/benchmarks/sift/src/matlab/siftmatch.m create mode 100755 SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexa64 create mode 100755 SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexglx create mode 100644 SD-VBS/benchmarks/sift/src/matlab/siftormx.c create mode 100755 SD-VBS/benchmarks/sift/src/matlab/siftormx.mexa64 create mode 100755 SD-VBS/benchmarks/sift/src/matlab/siftormx.mexglx create mode 100644 SD-VBS/benchmarks/sift/src/matlab/siftread.m create mode 100644 SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.c create mode 100644 SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.m create mode 100755 SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.mexa64 create mode 100755 SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.mexglx create mode 100644 SD-VBS/benchmarks/sift/src/matlab/tightsubplot.m (limited to 'SD-VBS/benchmarks/sift/src') diff --git a/SD-VBS/benchmarks/sift/src/c/diffss.c b/SD-VBS/benchmarks/sift/src/c/diffss.c new file mode 100644 index 0000000..683cbe7 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/c/diffss.c @@ -0,0 +1,53 @@ +/******************************** +Author: Sravanthi Kota Venkata +********************************/ + +#include "sift.h" + +/** + DIFFSS Difference of scale space + Returns a scale space DSS obtained by subtracting + consecutive levels of the scale space SS. + + In SIFT, this function is used to compute the difference of + Gaussian scale space from the Gaussian scale space of an image. +**/ + +F2D** diffss(F2D** ss, int num, int intervals) +{ + F2D** dss; + int o, sizeM, sizeN, s, i, j; + F2D *current, *in1, *in2; + + dss = malloc(num*intervals*sizeof(F2D*)); + + for(o=0; oheight; + sizeN = ss[o*intervals+s]->width; + + dss[o*intervals+s] = fMallocHandle(sizeM, sizeN); + + current = dss[o*intervals+s]; + in1 = ss[o*intervals+s+1]; + in2 = ss[o*intervals+s]; + + for(i=0; iheight; + N = I->width; + J = fSetArray(2*M, 2*N, 0); + + for(i=0; iwidth; i++) + { + if(asubsref(oframes,i)>3 && asubsref(oframes,i)<(N-3) && + subsref(oframes,1,i)>3 && subsref(oframes,1,i)<(M-3)) + { + cnt++; + } + } + + sel = iSetArray(cnt, 1, 0); + for(i=0; iwidth; i++) + { + if(asubsref(oframes,i)>3 && asubsref(oframes,i)<(N-3) && + subsref(oframes,1,i)>3 && subsref(oframes,1,i)<(M-3)) + { + asubsref(sel,k) = i; + k++; + } + m++; + } + + if( sel->height > 0) + { + ret = fSetArray(oframes->height, sel->height, 0); + { + for(i=0; iheight; i++) + { + subsref(ret,0,i) = subsref(oframes,0,asubsref(sel,i)); + subsref(ret,1,i) = subsref(oframes,1,asubsref(sel,i)); + subsref(ret,2,i) = subsref(oframes,2,asubsref(sel,i)); + } + } + } + else + ret = fSetArray(1,1,0); + + iFreeHandle(sel); + return ret; +} + + + + + + diff --git a/SD-VBS/benchmarks/sift/src/c/gaussianss.c b/SD-VBS/benchmarks/sift/src/c/gaussianss.c new file mode 100644 index 0000000..52dd95b --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/c/gaussianss.c @@ -0,0 +1,185 @@ +/******************************** +Author: Sravanthi Kota Venkata +********************************/ + +#include "sift.h" + +F2D* resizeArray(F2D* array, int omin) +{ + F2D* prev = NULL; + F2D* current = array; + int o; + if(omin<0) + { + for(o=1; o>=-omin; o--) + { + prev = current; + current = doubleSize(current); + fFreeHandle(prev); + } + } + if(omin>0) + { + for(o=1; o<= omin; o++) + { + prev = current; + current = halveSize(current); + fFreeHandle(prev); + } + } + return current; +} + +/** + Returns the Gaussian scale space of image I. Image I is assumed to be + pre-smoothed at level SIGMAN. O,S,OMIN,SMIN,SMAX,SIGMA0 are the + parameters of the scale space. +**/ + +F2D** gaussianss(F2D* array, float sigman, int O, int S, int omin, int smin, int smax, float sigma0) +{ + /* We compute the following items in the function + 1. Smooth input image per octave + 2. Smooth each octave for different intervals + 3. Subtract each "interval-1" smooth image from "interval" image per octave. So, per octave, we have "interval" * DOG images. + 4. So, octave * intervals * image + 5. Note: At each octave, the image is scaled down in both x and y directions + */ + + float k, dsigma0, dsigma; + int s, i, j, o, so, M, N, sbest; + int intervals = smax-smin+1; + float temp, target_sigma, prev_sigma; + F2D *TMP, **gss; + F2D* I = array; + + // Scale multiplicative step + k = pow(2, (1.0/S)); + dsigma0 = sigma0 * sqrt(1-(1.0/pow(k,2))); // Scale step factor + + // If omin < 0, multiply the size of the image. + I = resizeArray(I, omin); + M = I->height; + N = I->width; + so = -smin+1; // Index offset + + gss = malloc(O*intervals*sizeof(F2D*)); + if(gss == NULL) + { + printf("Could not allocate memory\n"); + return NULL; + } + +/** + First octave +-------------------------------------------------------------------- + +The first level of the first octave has scale index (o,s) = +(omin,smin) and scale coordinate + + sigma(omin,smin) = sigma0 2^omin k^smin + +The input image I is at nominal scale sigman. Thus in order to get +the first level of the pyramid we need to apply a smoothing of + + sqrt( (sigma0 2^omin k^smin)^2 - sigman^2 ). + +As we have pre-scaled the image omin octaves (up or down, +depending on the sign of omin), we need to correct this value +by dividing by 2^omin, getting + + sqrt( (sigma0 k^smin)^2 - (sigman/2^omin)^2 ) + +**/ + + temp = sqrt(pow((sigma0*pow(k,smin)),2) - pow((sigman/pow(2,omin)),2)); + + { + gss[0] = fSetArray(I->height, I->width, 0); + imsmooth(I, temp, gss[0] ); + } + + for(s=smin; sheight, gss[s+so-1]->width, 0); + imsmooth( gss[(s+so-1)] , dsigma, gss[(s+so)] ); + } + + /** Other octaves **/ + + for(o=1; o smax + + The amount of extra smoothing we need to apply is then given by + + ( sigma0 2^o 2^(smin/S) )^2 - ( sigma0 2^o 2^(sbest/S - 1) )^2 + + As usual, we divide by 2^o to cancel out the effect of the + downsampling and we get + + ( sigma 0 k^smin )^2 - ( sigma0 2^o k^(sbest - S) )^2 + **/ + + sbest = MIN(smin+S-1, smax-1); + TMP = halveSize( gss[(o-1)*intervals+sbest+so]); + target_sigma = sigma0 * pow(k,smin); + prev_sigma = sigma0 * pow(k, (sbest+1)-S); + + temp = sqrt(pow(target_sigma,2) - pow(prev_sigma, 2)); + if(target_sigma > prev_sigma) + { + gss[o*intervals] = fSetArray(TMP->height, TMP->width, 0); + imsmooth(TMP, temp, gss[o*intervals] ); + } + else + { + int i; + gss[o*intervals] = fSetArray(TMP->height, TMP->width, 0); + for(i=0; i<(TMP->height*TMP->width); i++) + asubsref(gss[o*intervals],i) = asubsref(TMP,i); + } + + M = TMP->height; + N = TMP->width; + + fFreeHandle(TMP); + + for(s=smin; sheight, gss[o*intervals+s-1+so]->width, 0); + imsmooth( gss[o*intervals+s-1+so] , dsigma, gss[o*intervals+s+so]); + } + } + + fFreeHandle(I); + return gss; +} diff --git a/SD-VBS/benchmarks/sift/src/c/halveSize.c b/SD-VBS/benchmarks/sift/src/c/halveSize.c new file mode 100644 index 0000000..fe1e536 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/c/halveSize.c @@ -0,0 +1,35 @@ +/******************************** +Author: Sravanthi Kota Venkata +********************************/ + +#include +#include +#include "sift.h" + +F2D* halveSize(F2D* I) +{ + F2D *J; + int i, j, k; + int hM, hN; + int M, N; + + M = I->height; + N = I->width; + + hM = (M+1)/2; + hN = (N+1)/2; + + J = fSetArray(hM, hN, 0.0); + + k = 0; + for(i=0; i +#include +const double win_factor = 1.5 ; +const int nbins = 36 ; +const float threshold = 0.01; + +/** + This function is similar to imageBlur in common/c folder. + Here, we can specify the sigma value for the gaussian filter + function. +**/ + +void imsmooth(F2D* array, float dsigma, F2D* out) +{ + int M,N ; + int i,j,k; + float s ; + + /* ------------------------------------------------------------------ + ** Check the arguments + ** --------------------------------------------------------------- */ + + M = array->height; + N = array->width; + s = dsigma; + + /* ------------------------------------------------------------------ + ** Do the job + ** --------------------------------------------------------------- */ + if(s > threshold) + { + int W = (int) ceil(4*s) ; + float temp[2*W+1]; + F2D* buffer; + float acc = 0.0; + + buffer = fSetArray(M,N,0); + + for(j = 0 ; j < (2*W+1) ; ++j) + { + temp[j] = (float)(expf(-0.5 * (j - W)*(j - W)/(s*s))) ; + acc += temp[j]; + } + + for(j = 0 ; j < (2*W+1) ; ++j) + { + temp[j] /= acc ; + } + + /* + ** Convolve along the columns + **/ + + for(j = 0 ; j < M ; ++j) + { + for(i = 0 ; i < N ; ++i) + { + int startCol = MAX(i-W,0); + int endCol = MIN(i+W, N-1); + int filterStart = MAX(0, W-i); + + assert(j < array->height); + assert(j < buffer->height); + assert(i < buffer->width); + for(k=startCol; k<=endCol; k++) { + assert(k < array->width); + assert(filterStart < 2*W+1); + subsref(buffer,j,i) += subsref(array, j, k) * temp[filterStart++]; + } + } + } + + /* + ** Convolve along the rows + **/ + for(j = 0 ; j < M ; ++j) + { + for(i = 0 ; i < N ; ++i) + { + int startRow = MAX(j-W,0); + int endRow = MIN(j+W, M-1); + int filterStart = MAX(0, W-j); + for(k=startRow; k<=endRow; k++) + subsref(out,j,i) += subsref(buffer,k,i) * temp[filterStart++]; + } + } + + fFreeHandle(buffer); + + } + else + { + for(i=0;i +#include +#include "sift.h" +#include +#include "extra.h" +#define SIFT_MEM 1<<29 +void normalizeImage(F2D* image) +{ + int i; + int rows; + int cols; + + int tempMin = 10000, tempMax = -1; + rows = image->height; + cols = image->width; + + for(i=0; i<(rows*cols); i++) + if(tempMin > asubsref(image,i)) + tempMin = asubsref(image,i); + + for(i=0; i<(rows*cols); i++) + asubsref(image,i) = asubsref(image,i) - tempMin; + + for(i=0; i<(rows*cols); i++) + if(tempMax < asubsref(image,i)) + tempMax = asubsref(image,i); + + for(i=0; i<(rows*cols); i++) + asubsref(image,i) = ( asubsref(image,i) / (tempMax+0.0) ); +} + +int main(int argc, char* argv[]) +{ + SET_UP + mallopt(M_TOP_PAD, SIFT_MEM); + mallopt(M_MMAP_MAX, 0); + I2D* im; + F2D *image; + int rows, cols, i, j; + F2D* frames; + unsigned int* startTime, *endTime, *elapsed; + + char imSrc[100]; + printf("Input image: "); + scanf("%s", imSrc); + im = readImage(imSrc); + image = fiDeepCopy(im); + rows = image->height; + cols = image->width; + + + printf("start\n"); + for_each_job{ + image = fiDeepCopy(im); + normalizeImage(image); + /** Extract sift features for the normalized image **/ + frames = sift(image); + } + printf("end..\n"); + +#ifdef CHECK + { + int ret=0; + float tol = 0.2; +#ifdef GENERATE_OUTPUT + fWriteMatrix(frames, argv[1]); +#endif + ret = fSelfCheck(frames, "expected_C.txt", tol); + if (ret == -1) + printf("Error in SIFT\n"); + } +#endif + + iFreeHandle(im); + fFreeHandle(frames); + WRITE_TO_FILE + return 0; +} + diff --git a/SD-VBS/benchmarks/sift/src/c/sift.c b/SD-VBS/benchmarks/sift/src/c/sift.c new file mode 100644 index 0000000..ea0d308 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/c/sift.c @@ -0,0 +1,314 @@ +/******************************** +Author: Sravanthi Kota Venkata +********************************/ + +#include +#include "sift.h" + +/** SIFT- Scale Invariant Feature Transform. This algorithm is based on + David Lowe's implementation of sift. So, we will use the parameter + values from Lowe's implementation. + See: http://www.cs.ubc.ca/~lowe/keypoints/ + + SIFT extracts from an image a collection of frames or keypoints. These + are oriented disks attacked to blob-like structures of the image. As + the image translates, rotates and scales, the frames track these blobs + and the deformation. + + 'BoundaryPoint' [bool] + If set to 1, frames too close to the image boundaries are discarded. + + 'Sigma0' [pixels] + Smoothing of the level 0 of octave 0 of the scale space. + (Note that Lowe's 1.6 value refers to the level -1 of octave 0.) + + 'SigmaN' [pixels] + Nominal smoothing of the input image. Typically set to 0.5. + + 'Threshold' + Threshold used to eliminate weak keypoints. Typical values for + intensity images in the range [0,1] are around 0.01. Smaller + values mean more keypoints. + +**/ + +F2D* sift(F2D* I) +{ + int rows, cols, K; + int subLevels, omin, Octaves, r, NBP, NBO, smin, smax, intervals, o; + float sigma, sigman, sigma0, thresh, magnif; + int discardBoundaryPoints; + F2D *descriptors; + F2D **gss, *tfr; + F2D **dogss; + I2D* i_, *j_, *s_; + F2D *oframes, *frames; + int i, j, k, m, startRow, startCol, endRow, endCol, firstIn=1; + float minVal; + I2D* tx1, *ty1, *ts1; + I2D* x1, *y1, *s1, *txys; + + rows = I->height; + cols = I->width; + + /** + Lowe's choices + Octaves - octave + subLevels - sub-level for image + sigma - sigma value for gaussian kernel, for smoothing the image + At each successive octave, the data is spatially downsampled by half + **/ + + subLevels = 3; + omin = -1; + minVal = log2f(MIN(rows,cols)); + Octaves = (int)(floor(minVal))-omin-4; /* Upto 16x16 images */ + sigma0 = 1.6 * pow(2, (1.0/subLevels)); + sigman = 0.5; + thresh = (0.04/subLevels)/2; + r = 10; + +#ifdef test + subLevels = 1; + Octaves = 1; + sigma0 = pow(0.6 * 2, 1); + sigman = 1.0; + thresh = (1/subLevels)/2; + r = 1; +#endif + + NBP = 4; + NBO = 8 ; + magnif = 3.0 ; + discardBoundaryPoints = 1 ; + + smin = -1; + smax = subLevels+1; + intervals = smax - smin + 1; + + /** + We build gaussian pyramid for the input image. Given image I, + we sub-sample the image into octave 'Octaves' number of levels. At + each level, we smooth the image with varying sigman values. + + Gaussiam pyramid can be assumed as a 2-D matrix, where each + element is an image. Number of rows corresponds to the number + of scales of the pyramid (octaves, "Octaves"). Row 0 (scale 0) is the + size of the actual image, Row 1 (scale 1) is half the actual + size and so on. + + At each scale, the image is smoothened with different sigma values. + So, each row has "intervals" number of smoothened images, starting + with least blurred. + + gss holds the entire gaussian pyramid. + **/ + + gss = gaussianss(I, sigman, Octaves, subLevels, omin, -1, subLevels+1, sigma0); + + /** + Once we build the gaussian pyramid, we compute DOG, the + Difference of Gaussians. At every scale, we do: + + dogss[fixedScale][0] = gss[fixedScale][1] - gss[fixedScale][0] + + Difference of gaussian gives us edge information, at each scale. + In order to detect keypoints on the intervals per octave, we + inspect DOG images at highest and lowest scales of the octave, for + extrema detection. + **/ + + dogss = diffss(gss, Octaves, intervals); + + /** The extraction of keypoints is carried one octave per time **/ + for(o=0; oheight; + int sizeCols = dogss[o*intervals]->width; + + { + int i,j,k=0; + temp = fMallocHandle(intervals-1, sizeRows*sizeCols); + negate = fMallocHandle(intervals-1, sizeRows*sizeCols); + + /** + Keypoints are detected as points of local extrema of the + DOG pyramid, at a given octave. In softlocalmax function, the + keypoints are extracted by looking at 9x9x9 neighborhood samples. + + We populate temp and negate arrays with the values of the DOG + pyramid for a given octave, o. Since we are interested in both + local maxima and minima, we compute negate matrix, which is the + negated values of the DOG pyramid. + + **/ + + for(i1=0; i1<(intervals-1); i1++) + { + for(j=0; jheight,idx->width,0); + y = iSetArray(idx->height,idx->width,0); + s_ = iSetArray(idx->height,idx->width,0); + + { + int i, j; + for(i=0; iheight; i++) + { + for(j=0; jwidth; j++) + { + int v, u, w, z; + w = subsref(idx,i,j); + v = ceil((w/(sizeRows*sizeCols)) + 0.5); + u = floor(w/(sizeRows*sizeCols)); + z = w - (sizeRows*sizeCols*u); + + /** v is the interval number, s **/ + subsref(s_,i,j) = v; + /** row number of the index **/ + subsref(y,i,j) = ceil((z / sizeRows)+0.5); + /** col number of the index **/ + subsref(x,i,j) = z - (sizeCols * floor(z / sizeRows)); + } + } + } + + { + + tx1 = isMinus(x, 1); + x1 = iReshape(tx1, 1, (tx1->height*tx1->width)); + + ty1 = isMinus(y, 1); + y1 = iReshape(ty1, 1, (ty1->height*ty1->width)); + + ts1 = isPlus(s_, (smin-1)); + s1 = iReshape(ts1, 1, (ts1->height*ts1->width)); + + txys = iVertcat(y1, x1); + i_ = iVertcat(txys, s1); + } + + /** + Stack all x,y,s into oframes. + Row 0 of oframes = x + Row 1 of oframes = y + Row 2 of oframes = s + **/ + oframes = fiDeepCopy(i_); + + { + int i,j; + F2D* temp; + temp = oframes; + + /** + Remove points too close to the boundary + **/ + + if(discardBoundaryPoints) + oframes = filterBoundaryPoints(sizeRows, sizeCols, temp); + fFreeHandle(temp); + } + + /** + Refine the location, threshold strength and remove points on edges + **/ + if( asubsref(oframes,0) != 0) + { + F2D* temp_; + temp_ = fTranspose(oframes); + fFreeHandle(oframes); + oframes = siftrefinemx(temp_, temp, smin, thresh, r, sizeRows, sizeCols, intervals-1); + fFreeHandle(temp_); + + if( firstIn == 0) + { + tfr = fDeepCopy(frames); + fFreeHandle(frames); + frames = fHorzcat(tfr, oframes); + fFreeHandle(tfr); + } + else + frames = fDeepCopy(oframes); + firstIn = 0; + } + else if(Octaves == 1) + frames = fDeepCopy(oframes); + + fFreeHandle(oframes); + iFreeHandle(y); + iFreeHandle(x); + iFreeHandle(s_); + iFreeHandle(y1); + iFreeHandle(x1); + iFreeHandle(s1); + iFreeHandle(ty1); + iFreeHandle(tx1); + iFreeHandle(ts1); + iFreeHandle(txys); + iFreeHandle(i_); + iFreeHandle(idx); + fFreeHandle(idxf); + fFreeHandle(idxft); + fFreeHandle(temp); + fFreeHandle(t); + fFreeHandle(negate); + } + + { int s; + for(o=0; o + +#define GREATER(a,b) ((a) > (b)) +#define MAX(a,b) (((a)>(b))?(a):(b)) +#define MIN(a,b) (((a)<(b))?(a):(b)) + +F2D* sift(F2D* I); +F2D* halveSize(F2D* I); +F2D** gaussianss(F2D* I, float sigman, int O, int S, int omin, int smin, int smax, float sigma0); +F2D** diffss(F2D** ss, int O, int intervals); +F2D* doubleSize(F2D* I); +void imsmooth(F2D* I_pt, float dsigma, F2D* out); +F2D* siftlocalmax(F2D* in, float thresh, int intervals, int M, int N); +F2D* filterBoundaryPoints(int M, int N, F2D* oframes); +F2D* siftrefinemx(F2D* oframes, F2D* dogss, int smin, float thresh, int rin, int M, int N, int intervals); + +#endif + + diff --git a/SD-VBS/benchmarks/sift/src/c/siftlocalmax.c b/SD-VBS/benchmarks/sift/src/c/siftlocalmax.c new file mode 100644 index 0000000..d75bcbd --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/c/siftlocalmax.c @@ -0,0 +1,247 @@ +#include "sift.h" + +/** + SIFTLOCALMAX Find local maximizers + Returns the indexes of the local maximizers of + the 3-dimensional array F. + + Say, we have a Q-dimensional array F. + A local maximizer is an element whose value is greater than the + value of all its neighbors. The neighbors of an element i1...iQ + are the subscripts j1...jQ such that iq-1 <= jq <= iq (excluding + i1...iQ itself). For example, if Q=1 the neighbors of an element + are its predecessor and successor in the linear order; if Q=2, its + neighbors are the elements immediately to its north, south, west, + est, north-west, north-est, south-west and south-est + (8-neighborhood). + + Points on the boundary of F are ignored (and never selected as + local maximizers). + + SEL=SIFTLOCALMAX(F,THRESH) accepts an element as a mazimizer only + if it is at least THRES greater than all its neighbors. + + SEL=SIFTLOCALMAX(F,THRESH,P) look for neighbors only in the first + P dimensions of the Q-dimensional array F. This is useful to + process F in ``slices''. + + REMARK. Matrices (2-array) with a singleton dimension are + interpreted as vectors (1-array). So for example SIFTLOCALMAX([0 1 + 0]) and SIFTLOCALMAX([0 1 0]') both return 2 as an aswer. However, + if [0 1 0] is to be interpreted as a 1x2 matrix, then the correct + answer is the empty set, as all elements are on the boundary. + Unfortunately MATLAB does not distinguish between vectors and + 2-matrices with a singleton dimension. To forece the + interpretation of all matrices as 2-arrays, use + SIFTLOCALMAX(F,TRESH,2) (but note that in this case the result is + always empty!). +**/ + +#define NDIMS 3 + +F2D* siftlocalmax(F2D* in, float thresh, int intervals, int M, int N) +{ + int ndims, ptoffset=0, maxIter = 0; + int offsets[NDIMS]; + int midx[NDIMS]; + int dims[NDIMS]; + I2D* neighbors ; + int nneighbors ; + F2D* out; + + /** + We pass dogss[o], which has >1 intervals + If >1 intervals, then the dimensions of in[F] will be 1 x intervals + If =1 interval, then the dimensions of in[F] will be the size of the dogss[o] image + We will check for this condition. + **/ + + ndims = NDIMS; /* Since we have intervals number of images of size M x N */ + dims[0] = M; + dims[1] = N; + dims[2] = intervals-1; + + /* + If there are only two dimensions and if one is singleton, then + assume that a vector has been provided as input (and treat this + as a COLUMN matrix with p=1). We do this because Matlab does not + distinguish between vectors and 1xN or Mx1 matrices and because + the cases 1xN and Mx1 are trivial (the result is alway empty). + */ + + + + /* ------------------------------------------------------------------ + * Do the job + * --------------------------------------------------------------- */ + int maxima_size = M*N ; + + I2D* maxima_start = iMallocHandle(1, maxima_size); + int* maxima_iterator = maxima_start->data ; + int* maxima_end = maxima_start->data + maxima_size ; + + int i,h,o,j; + F2D* pt; + + pt = in; + + /* Compute the offsets between dimensions. */ + offsets[0] = 1 ; + for(i = 1 ; i < ndims ; ++i) + { + offsets[i] = offsets[i-1]*dims[i-1] ; + } + + /* Multi-index. */ + for(i = 0 ; i < ndims ; ++i) + midx[i] = 1 ; + + /* Neighbors. */ + nneighbors = 1 ; + o=0 ; + for(i = 0 ; i < ndims ; ++i) + { + nneighbors *= 3 ; + midx[i] = -1 ; + o -= offsets[i] ; + } + nneighbors -= 1 ; + neighbors = iMallocHandle(1,nneighbors); + + /* Precompute offsets from offset(-1,...,-1,0,...0) to + * offset(+1,...,+1,0,...,0). */ + i = 0 ; + + while(1) + { + if(o != 0) + { + asubsref(neighbors, i) = o ; + i++; + } + h = 0 ; + while( o += offsets[h], (++midx[h]) > 1 ) + { + o -= 3*offsets[h] ; + midx[h] = -1 ; + if(++h >= ndims) + goto stop ; + } + } + + stop: ; + + /* Starts at the corner (1,1,...,1,0,0,...0) */ + for(h = 0 ; h < ndims ; ++h) + { + midx[h] = 1 ; + ptoffset += offsets[h] ; + } + + + /* --------------------------------------------------------------- + * Loop + * ------------------------------------------------------------ */ + + /* + If any dimension in the first P is less than 3 elements wide + then just return the empty matrix (if we proceed without doing + anything we break the carry reporting algorithm below). + */ + + for(h=0 ; h < ndims ; ++h) + if(dims[h] < 3) + goto end ; + + while(1) + { + /* Propagate carry along multi index midx */ + + h = 0 ; + while(midx[h] >= dims[h] - 1) + { + midx[h] = 1 ; + if(++h >= ndims) + goto next_layer ; + ++midx[h] ; + } + + /* Scan neighbors */ + { + float v = asubsref(pt, ptoffset); //*pt ; + int is_greater = (v>=thresh) ? 1 : 0; + + assert(ptoffset < pt->width*pt->height); + + i = 0 ; + while(is_greater && i < nneighbors) + { + if(v > asubsref(pt, ptoffset + asubsref(neighbors,i))) + is_greater *= 1; + else + is_greater *= 0; + i = i+1; + } + + /* Add the local maximum */ + if(is_greater) + { + /* Need more space? */ + if( &(maxima_iterator[maxIter]) == maxima_end) + { + int *temp, i, j; + maxima_size += M*N ; + + free(maxima_start); + maxima_start = iMallocHandle(1, maxima_size); + + maxima_end = maxima_start->data + maxima_size ; + maxima_iterator = maxima_end - M*N ; + maxIter = 0; + } + + maxima_iterator[maxIter] = ptoffset + 1; + maxIter++; + } + + /* Go to next element */ + ptoffset += 1 ; + ++midx[0] ; + continue ; + + next_layer: ; + if( h >= ndims ) + goto end ; + + while((++midx[h]) >= dims[h]) + { + midx[h] = 0 ; + if(++h >= ndims) + goto end ; + } + } + } + + end:; + /* Return. */ + { + int i=0; + out = fMallocHandle(1, maxIter); + + for(i=0; iheight; + + /* If the input array is empty, then output an empty array as well. */ + if( K == 0) { + out = fDeepCopy(oframes); + return out; + } + + /* ----------------------------------------------------------------- + * Do the job + * -------------------------------------------------------------- */ + { + F2D* buffer = fMallocHandle(K,3); + int p ; + const int yo = 1 ; + const int xo = M ; + const int so = M*N; + int buffCounter = 0; + int pptcount = 0; + + for(p = 0 ; p < K ; ++p) + { + float tx, ty, ts; + int x,y,s; + int iter ; + float b[3] ; + + tx = asubsref(oframes, pptcount++); + ty = asubsref(oframes, pptcount++); + ts = asubsref(oframes, pptcount++); + + x = ceil(tx); + y = ceil(ty); + s = ceil(ts) - smin; + + + /* Local maxima extracted from the DOG + * have coorrinates 1<=x<=N-2, 1<=y<=M-2 + * and 1<=s-mins<=S-2. This is also the range of the points + * that we can refine. + */ + + if(x < 1 || x > N-2 || + y < 1 || y > M-2 || + s < 1 || s > S-2) { + continue ; + } + +#define at(dx,dy,ds) asubsref(dogss, (dx)*xo + (dy)*yo + (ds)*so) + + { + float Dx=0,Dy=0,Ds=0,Dxx=0,Dyy=0,Dss=0,Dxy=0,Dxs=0,Dys=0 ; + int dx = 0 ; + int dy = 0 ; + + for(iter = 0 ; iter < max_iter ; ++iter) + { + + float A[3*3] ; + +#define Aat(i,j) (A[(i)+(j)*3]) + + x += dx ; + y += dy ; + + /* Compute the gradient. */ + Dx = 0.5 * (at(x+1,y+0,s+0) - at(x-1,y+0,s+0)) ; + Dy = 0.5 * (at(x+0,y+1,s+0) - at(x+0,y-1,s+0)); + Ds = 0.5 * (at(x+0,y+0,s+1) - at(x+0,y+0,s-1)) ; + + /* Compute the Hessian. */ + Dxx = (at(x+1,y,s) + at(x-1,y,s) - 2.0 * at(x,y,s)) ; + Dyy = (at(x,y+1,s) + at(x,y-1,s) - 2.0 * at(x,y,s)) ; + Dss = (at(x,y,s+1) + at(x,y,s-1) - 2.0 * at(x,y,s)) ; + + Dxy = 0.25 * ( at(x+1,y+1,s) + at(x-1,y-1,s) - at(x-1,y+1,s) - at(x+1,y-1,s) ) ; + Dxs = 0.25 * ( at(x+1,y,s+1) + at(x-1,y,s-1) - at(x-1,y,s+1) - at(x+1,y,s-1) ) ; + Dys = 0.25 * ( at(x,y+1,s+1) + at(x,y-1,s-1) - at(x,y-1,s+1) - at(x,y+1,s-1) ) ; + + /* Solve linear system. */ + Aat(0,0) = Dxx ; + Aat(1,1) = Dyy ; + Aat(2,2) = Dss ; + Aat(0,1) = Aat(1,0) = Dxy ; + Aat(0,2) = Aat(2,0) = Dxs ; + Aat(1,2) = Aat(2,1) = Dys ; + + b[0] = - Dx ; + b[1] = - Dy ; + b[2] = - Ds ; + + /* If the translation of the keypoint is big, move the keypoint + * and re-iterate the computation. Otherwise we are all set. + */ + dx= ((b[0] > 0.6 && x < N-2) ? 1 : 0 ) + + ((b[0] < -0.6 && x > 1 ) ? -1 : 0 ) ; + + dy= ((b[1] > 0.6 && y < M-2) ? 1 : 0 ) + + ((b[1] < -0.6 && y > 1 ) ? -1 : 0 ) ; + + if( dx == 0 && dy == 0 ) + break ; + + } + + { + float val = at(x,y,s) + 0.5 * (Dx * b[0] + Dy * b[1] + Ds * b[2]) ; + float score = (Dxx+Dyy)*(Dxx+Dyy) / (Dxx*Dyy - Dxy*Dxy) ; + float xn = x + b[0] ; + float yn = y + b[1] ; + float sn = s + b[2] ; + + if(fabs(val) > thresh && + score < (r+1)*(r+1)/r && + score >= 0 && + fabs(b[0]) < 1.5 && + fabs(b[1]) < 1.5 && + fabs(b[2]) < 1.5 && + xn >= 0 && + xn <= N-1 && + yn >= 0 && + yn <= M-1 && + sn >= 0 && + sn <= S-1) + { + asubsref(buffer,buffCounter++) = xn ; + asubsref(buffer,buffCounter++) = yn ; + asubsref(buffer,buffCounter++) = sn+smin ; + } + } + } + } + + /* Copy the result into an array. */ + { + int i, j, k=0; + int NL = buffCounter/3; + out = fMallocHandle(3, NL); + + for(i=0; i 0 + for o=1:omin + I = halveSize(I) ; + end +end + +[M,N] = size(I) ; + +% Index offset +so = -smin+1 ; + +% -------------------------------------------------------------------- +% First octave +% -------------------------------------------------------------------- +% +% The first level of the first octave has scale index (o,s) = +% (omin,smin) and scale coordinate +% +% sigma(omin,smin) = sigma0 2^omin k^smin +% +% The input image I is at nominal scale sigman. Thus in order to get +% the first level of the pyramid we need to apply a smoothing of +% +% sqrt( (sigma0 2^omin k^smin)^2 - sigman^2 ). +% +% As we have pre-scaled the image omin octaves (up or down, +% depending on the sign of omin), we need to correct this value +% by dividing by 2^omin, getting +%e +% sqrt( (sigma0 k^smin)^2 - (sigman/2^omin)^2 ) +% + +SS.octave{1} = zeros(M,N,smax-smin+1) ; +SS.octave{1}(:,:,1) = imsmooth(I, ... + sqrt((sigma0*k^smin)^2 - (sigman/2^omin)^2)) ; + +temp = sqrt((sigma0*k^smin)^2 - (sigman/2^omin)^2) ; + +for s=smin+1:smax + % Here we go from (omin,s-1) to (omin,s). The extra smoothing + % standard deviation is + % + % (sigma0 2^omin 2^(s/S) )^2 - (simga0 2^omin 2^(s/S-1/S) )^2 + % + % Aftred dividing by 2^omin (to take into account the fact + % that the image has been pre-scaled omin octaves), the + % standard deviation of the smoothing kernel is + % + % dsigma = sigma0 k^s sqrt(1-1/k^2) + % + dsigma = k^s * dsigma0 ; + SS.octave{1}(:,:,s +so) = ... + imsmooth(squeeze(... + SS.octave{1}(:,:,s-1 +so)... + ), dsigma ) ; +end + +% -------------------------------------------------------------------- +% Other octaves +% -------------------------------------------------------------------- + +for o=2:O + % We need to initialize the first level of octave (o,smin) from + % the closest possible level of the previous octave. A level (o,s) + % in this octave corrsponds to the level (o-1,s+S) in the previous + % octave. In particular, the level (o,smin) correspnds to + % (o-1,smin+S). However (o-1,smin+S) might not be among the levels + % (o-1,smin), ..., (o-1,smax) that we have previously computed. + % The closest pick is + % + % / smin+S if smin+S <= smax + % (o-1,sbest) , sbest = | + % \ smax if smin+S > smax + % + % The amount of extra smoothing we need to apply is then given by + % + % ( sigma0 2^o 2^(smin/S) )^2 - ( sigma0 2^o 2^(sbest/S - 1) )^2 + % + % As usual, we divide by 2^o to cancel out the effect of the + % downsampling and we get + % + % ( sigma 0 k^smin )^2 - ( sigma0 2^o k^(sbest - S) )^2 + % + sbest = min(smin + S, smax) ; + TMP = halveSize(squeeze(SS.octave{o-1}(:,:,sbest+so))) ; + + target_sigma = sigma0 * k^smin ; + prev_sigma = sigma0 * k^(sbest - S) ; + + if(target_sigma > prev_sigma) + temp = sqrt(target_sigma^2 - prev_sigma^2); + TMP = imsmooth(TMP, temp ) ; + end + [M,N] = size(TMP) ; + + SS.octave{o} = zeros(M,N,smax-smin+1) ; + SS.octave{o}(:,:,1) = TMP ; + + for s=smin+1:smax + % The other levels are determined as above for the first octave. + dsigma = k^s * dsigma0 ; + SS.octave{o}(:,:,s +so) = ... + imsmooth(squeeze(... + SS.octave{o}(:,:,s-1 +so)... + ), dsigma) ; + end + + +end + +% ------------------------------------------------------------------------- +% Auxiliary functions +% ------------------------------------------------------------------------- +function J = doubleSize(I) +[M,N]=size(I) ; +J = zeros(2*M,2*N) ; +J(1:2:end,1:2:end) = I ; +J(2:2:end-1,2:2:end-1) = ... + 0.25*I(1:end-1,1:end-1) + ... + 0.25*I(2:end,1:end-1) + ... + 0.25*I(1:end-1,2:end) + ... + 0.25*I(2:end,2:end) ; +J(2:2:end-1,1:2:end) = ... + 0.5*I(1:end-1,:) + ... + 0.5*I(2:end,:) ; +J(1:2:end,2:2:end-1) = ... + 0.5*I(:,1:end-1) + ... + 0.5*I(:,2:end) ; + +function J = halveSize(I) +J=I(1:2:end,1:2:end) ; +%[M,N] = size(I) ; +%m=floor((M+1)/2) ; +%n=floor((N+1)/2) ; +%J = I(:,1:2:2*n) + I(:,2:2:2*n+1) ; +%J = 0.25*(J(1:2:2*m,:)+J(2:2:2*m+1,:)) ; diff --git a/SD-VBS/benchmarks/sift/src/matlab/imreadbw.m b/SD-VBS/benchmarks/sift/src/matlab/imreadbw.m new file mode 100644 index 0000000..55cb708 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/imreadbw.m @@ -0,0 +1,52 @@ +function I = imreadbw(file) +% IMREADBW Reads an image as gray-scale +% I=IMREADBW(FILE) reads the image FILE and converts the result to a +% gray scale image (with DOUBLE storage class anr range normalized +% in [0,1]). + +% AUTORIGHTS +% Copyright (c) 2006 The Regents of the University of California. +% All Rights Reserved. +% +% Created by Andrea Vedaldi +% UCLA Vision Lab - Department of Computer Science +% +% Permission to use, copy, modify, and distribute this software and its +% documentation for educational, research and non-profit purposes, +% without fee, and without a written agreement is hereby granted, +% provided that the above copyright notice, this paragraph and the +% following three paragraphs appear in all copies. +% +% This software program and documentation are copyrighted by The Regents +% of the University of California. The software program and +% documentation are supplied "as is", without any accompanying services +% from The Regents. The Regents does not warrant that the operation of +% the program will be uninterrupted or error-free. The end-user +% understands that the program was developed for research purposes and +% is advised not to rely exclusively on the program for any reason. +% +% This software embodies a method for which the following patent has +% been issued: "Method and apparatus for identifying scale invariant +% features in an image and use of same for locating an object in an +% image," David G. Lowe, US Patent 6,711,293 (March 23, +% 2004). Provisional application filed March 8, 1999. Asignee: The +% University of British Columbia. +% +% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. + +I=im2double(imread(file)) ; + +if(size(I,3) > 1) + I = rgb2gray( I ) ; +end + + diff --git a/SD-VBS/benchmarks/sift/src/matlab/imsmooth.c b/SD-VBS/benchmarks/sift/src/matlab/imsmooth.c new file mode 100644 index 0000000..5d6a660 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/imsmooth.c @@ -0,0 +1,115 @@ +/** file: imsmooth.c + ** author: Andrea Vedaldi + ** description: Smooth an image. + **/ + +#include"mex.h" +#include +#include +#include +#include + +#define greater(a,b) ((a) > (b)) +#define min(a,b) (((a)<(b))?(a):(b)) +#define max(a,b) (((a)>(b))?(a):(b)) + +const double win_factor = 1.5 ; +const int nbins = 36 ; + +void +mexFunction(int nout, mxArray *out[], + int nin, const mxArray *in[]) +{ + int M,N ; + double* I_pt ; + double* J_pt ; + double s ; + enum {I=0,S} ; + enum {J=0} ; + + /* ------------------------------------------------------------------ + ** Check the arguments + ** --------------------------------------------------------------- */ + if (nin != 2) { + mexErrMsgTxt("Exactly two input arguments required."); + } else if (nout > 1) { + mexErrMsgTxt("Too many output arguments."); + } + + if (!mxIsDouble(in[I]) || + !mxIsDouble(in[S]) ) { + mexErrMsgTxt("All arguments must be real.") ; + } + + if(mxGetNumberOfDimensions(in[I]) > 2|| + mxGetNumberOfDimensions(in[S]) > 2) { + mexErrMsgTxt("I must be a two dimensional array and S a scalar.") ; + } + + if(max(mxGetM(in[S]),mxGetN(in[S])) > 1) { + mexErrMsgTxt("S must be a scalar.\n") ; + } + + M = mxGetM(in[I]) ; + N = mxGetN(in[I]) ; + + out[J] = mxCreateDoubleMatrix(M, N, mxREAL) ; + + I_pt = mxGetPr(in[I]) ; + J_pt = mxGetPr(out[J]) ; + s = *mxGetPr(in[S]) ; + + /* ------------------------------------------------------------------ + ** Do the job + ** --------------------------------------------------------------- */ + if(s > 0.01) { + + int W = (int) ceil(4*s) ; + int i ; + int j ; + double* g0 = (double*) mxMalloc( (2*W+1)*sizeof(double) ) ; + double* buffer = (double*) mxMalloc( M*N*sizeof(double) ) ; + double acc=0.0 ; + + for(j = 0 ; j < 2*W+1 ; ++j) { + g0[j] = exp(-0.5 * (j - W)*(j - W)/(s*s)) ; + acc += g0[j] ; + } + for(j = 0 ; j < 2*W+1 ; ++j) { + g0[j] /= acc ; + } + + /* + ** Convolve along the columns + **/ + for(j = 0 ; j < N ; ++j) { + for(i = 0 ; i < M ; ++i) { + double* start = I_pt + max(i-W,0) + j*M ; + double* stop = I_pt + min(i+W,M-1) + j*M + 1 ; + double* g = g0 + max(0, W-i) ; + acc = 0.0 ; + while(stop != start) acc += (*g++) * (*start++) ; + *buffer++ = acc ; + } + } + buffer -= M*N ; + + /* + ** Convolve along the rows + **/ + for(j = 0 ; j < N ; ++j) { + for(i = 0 ; i < M ; ++i) { + double* start = buffer + i + max(j-W,0)*M ; + double* stop = buffer + i + min(j+W,N-1)*M + M ; + double* g = g0 + max(0, W-j) ; + acc = 0.0 ; + while(stop != start) { acc += (*g++) * (*start) ; start+=M ;} + *J_pt++ = acc ; + } + } + mxFree(buffer) ; + mxFree(g0) ; + } else { + memcpy(J_pt, I_pt, sizeof(double)*M*N) ; + } +} diff --git a/SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexa64 b/SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexa64 new file mode 100755 index 0000000..0997262 Binary files /dev/null and b/SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexa64 differ diff --git a/SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexglx b/SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexglx new file mode 100755 index 0000000..a4a7a20 Binary files /dev/null and b/SD-VBS/benchmarks/sift/src/matlab/imsmooth.mexglx differ diff --git a/SD-VBS/benchmarks/sift/src/matlab/imsmooth_.c b/SD-VBS/benchmarks/sift/src/matlab/imsmooth_.c new file mode 100644 index 0000000..e0adeff --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/imsmooth_.c @@ -0,0 +1,103 @@ +#include +#include +#include +#include "imsmooth.h" + +#define greater(a,b) ((a) > (b)) +#define min(a,b) (((a)<(b))?(a):(b)) +#define max(a,b) (((a)>(b))?(a):(b)) + +const double win_factor = 1.5 ; +const int nbins = 36 ; + +float* imsmooth(float* I_pt_, float dsigma) +{ + int M,N ; + float *I_pt ; + float* J_pt_, *J_pt ; + float* out_, *out; + float s ; + enum {I=0,S} ; + enum {J=0} ; + + /* ------------------------------------------------------------------ + ** Check the arguments + ** --------------------------------------------------------------- */ + + M = (int)in_[0]; + N = (int)in_[1]; + + out_ = fMallocHandle(M, N); + J_pt_ = &(out_[0]); + J_pt = &(J_pt_[4]); + s = dsigma; + + + /* ------------------------------------------------------------------ + ** Do the job + ** --------------------------------------------------------------- */ + if(s > 0.01) + { + + int W = (int) ceil(4*s) ; + int i ; + int j ; + float* g0_, *g0, *buffer_, *buffer, acc; + g0_ = fMallocHandle(1, 2*W+1); + g0 = &(g0_[4]); + buffer_ = fMallocHandle(M, N); + acc=0.0 ; + + for(j = 0 ; j < 2*W+1 ; ++j) + { + g0[j] = exp(-0.5 * (j - W)*(j - W)/(s*s)) ; + acc += g0[j] ; + } + + for(j = 0 ; j < 2*W+1 ; ++j) + { + g0[j] /= acc ; + } + + /* + ** Convolve along the columns + **/ + + for(j = 0 ; j < N ; ++j) + { + for(i = 0 ; i < M ; ++i) + { + float* start = I_pt + max(i-W,0) + j*M ; + float* stop = I_pt + min(i+W,M-1) + j*M + 1 ; + float* g = g0 + max(0, W-i) ; + acc = 0.0 ; + while(stop != start) acc += (*g++) * (*start++) ; + *buffer++ = acc ; + } + } + buffer -= M*N ; + + /* + ** Convolve along the rows + **/ + for(j = 0 ; j < N ; ++j) + { + for(i = 0 ; i < M ; ++i) + { + float* start = buffer + i + max(j-W,0)*M ; + float* stop = buffer + i + min(j+W,N-1)*M + M ; + float* g = g0 + max(0, W-j) ; + acc = 0.0 ; + while(stop != start) { acc += (*g++) * (*start) ; start+=M ;} + *J_pt++ = acc ; + } + } + free(buffer) ; + free(g0) ; + } + else + { + memcpy(J_pt, I_pt, sizeof(double)*M*N) ; + } + return J_pt_; +} diff --git a/SD-VBS/benchmarks/sift/src/matlab/mexutils.c b/SD-VBS/benchmarks/sift/src/matlab/mexutils.c new file mode 100644 index 0000000..45843cf --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/mexutils.c @@ -0,0 +1,98 @@ +/* file: mexutils.c +** author: Andrea Vedaldi +** description: Utility functions to write MEX files. +**/ + +#include"mex.h" +#include + +#undef M_PI +#define M_PI 3.14159265358979 + +/** @biref Is real scalar? + ** + ** @return @c true if the array @a A is a real scalar. + **/ +int +uIsRealScalar(const mxArray* A) +{ + return + mxIsDouble(A) && + !mxIsComplex(A) && + mxGetNumberOfDimensions(A) == 2 && + mxGetM(A) == 1 && + mxGetN(A) == 1 ; +} + +/** @brief Is real matrix? + ** + ** The function checks wether the argument @a A is a real matrix. In + ** addition, if @a M >= 0, it checks wether the number of rows is + ** equal to @a M and, if @a N >= 0, if the number of columns is equal + ** to @a N. + ** + ** @param M number of rows. + ** @param N number of columns. + ** @return @c true if the array is a real matrix with the specified format. + **/ +int +uIsRealMatrix(const mxArray* A, int M, int N) +{ + return + mxIsDouble(A) && + !mxIsComplex(A) && + mxGetNumberOfDimensions(A) == 2 && + ((M>=0)?(mxGetM(A) == M):1) && + ((N>=0)?(mxGetN(A) == N):1) ; +} + +/** @brief Is real vector? + ** + ** The function checks wether the argument @a V is a real vector. By + ** definiton, a matrix is a vector if one of its dimension is one. + ** In addition, if @a D >= 0, it checks wether the dimension of the + ** vecotr is equal to @a D. + ** + ** @param D lenght of the vector. + ** @return @c true if the array is a real vector of the specified dimension. + **/ +int +uIsRealVector(const mxArray* V, int D) +{ + int M = mxGetM(V) ; + int N = mxGetN(V) ; + int is_vector = (N == 1) || (M == 1) ; + + return + mxIsDouble(V) && + !mxIsComplex(V) && + mxGetNumberOfDimensions(V) == 2 && + is_vector && + ( D < 0 || N == D || M == D) ; +} + + +/** @brief Is a string? + ** + ** The function checks wether the array @a S is a string. If + ** @a L is non-negative, it also check wether the strign has + ** length @a L. + ** + ** @return @a c true if S is a string of the specified length. + **/ +int +uIsString(const mxArray* S, int L) +{ + int M = mxGetM(S) ; + int N = mxGetN(S) ; + + return + mxIsChar(S) && + M == 1 && + (L < 0 || N == L) ; +} + +/** + ** + **/ + diff --git a/SD-VBS/benchmarks/sift/src/matlab/plotmatches.m b/SD-VBS/benchmarks/sift/src/matlab/plotmatches.m new file mode 100644 index 0000000..b6324e9 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/plotmatches.m @@ -0,0 +1,124 @@ +function h=plotmatches(I1,I2,P1,P2,matches,varargin) +% PLOTMATCHES Plot keypoint matches +% PLOTMATCHES(I1,I2,P1,P2,MATCHES) plots the two images I1 and I2 +% and lines connecting the frames (keypoints) P1 and P2 as specified +% by MATCHES. +% +% P1 and P2 specify two sets of frames, one per column. The first +% two elements of each column specify the X,Y coordinates of the +% corresponding frame. Any other element is ignored. +% +% MATCHES specifies a set of matches, one per column. The two +% elementes of each column are two indexes in the sets P1 and P2 +% respectively. +% +% The images I1 and I2 might be either both grayscale or both color +% and must have DOUBLE storage class. If they are color the range +% must be normalized in [0,1]. +% +% The function accepts the following option-value pairs: +% +% 'Stacking' ['h'] +% Stacking of images: horizontal [h], vertical [v], diagonal +% [h], overlap ['o'] +% +% See also PLOTSIFTDESCRIPTOR(), PLOTSIFTFRAME(), PLOTSS(). + +% AUTORIGHTS +% Copyright (c) 2006 The Regents of the University of California. +% All Rights Reserved. +% +% Created by Andrea Vedaldi +% UCLA Vision Lab - Department of Computer Science +% +% Permission to use, copy, modify, and distribute this software and its +% documentation for educational, research and non-profit purposes, +% without fee, and without a written agreement is hereby granted, +% provided that the above copyright notice, this paragraph and the +% following three paragraphs appear in all copies. +% +% This software program and documentation are copyrighted by The Regents +% of the University of California. The software program and +% documentation are supplied "as is", without any accompanying services +% from The Regents. The Regents does not warrant that the operation of +% the program will be uninterrupted or error-free. The end-user +% understands that the program was developed for research purposes and +% is advised not to rely exclusively on the program for any reason. +% +% This software embodies a method for which the following patent has +% been issued: "Method and apparatus for identifying scale invariant +% features in an image and use of same for locating an object in an +% image," David G. Lowe, US Patent 6,711,293 (March 23, +% 2004). Provisional application filed March 8, 1999. Asignee: The +% University of British Columbia. +% +% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. + +% -------------------------------------------------------------------- +% Check the arguments +% -------------------------------------------------------------------- + +stack='h' ; + +for k=1:2:length(varargin) + switch varargin{k} + case 'Stacking' + stack=varargin{k+1} ; + end +end + +% -------------------------------------------------------------------- +% Do the job +% -------------------------------------------------------------------- + +[M1,N1,K1]=size(I1) ; +[M2,N2,K2]=size(I2) ; + +switch stack + case 'h' + N3=N1+N2 ; + M3=max(M1,M2) ; + oj=N1 ; + oi=0 ; + case 'v' + M3=M1+M2 ; + N3=max(N1,N2) ; + oj=0 ; + oi=M1 ; + case 'd' + M3=M1+M2 ; + N3=N1+N2 ; + oj=N1 ; + oi=M1 ; + case 'o' + M3=max(M1,M2) ; + N3=max(N1,N2) ; + oj=0; + oi=0; +end + +I=zeros(M3,N3,K1) ; +I(1:M1,1:N1,:) = I1 ; +I(oi+(1:M2),oj+(1:N2),:) = I2 ; + +axes('Position', [0 0 1 1]) ; +imagesc(I) ; colormap gray ; hold on ; axis image ; axis off ; +drawnow ; + +K = size(matches, 2) ; +nans = NaN * ones(1,K) ; + +x = [ P1(1,matches(1,:)) ; P2(1,matches(2,:))+oj ; nans ] ; +y = [ P1(2,matches(1,:)) ; P2(2,matches(2,:))+oi ; nans ] ; + +h = line(x(:)', y(:)') ; +set(h,'Marker','.','Color','g') ; diff --git a/SD-VBS/benchmarks/sift/src/matlab/plotsiftdescriptor.m b/SD-VBS/benchmarks/sift/src/matlab/plotsiftdescriptor.m new file mode 100644 index 0000000..2f4af50 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/plotsiftdescriptor.m @@ -0,0 +1,169 @@ +function h=plotsiftdescriptor(d,f) +% PLOTSIFTDESCRIPTOR Plot SIFT descriptor +% PLOTSIFTDESCRIPTOR(D) plots the SIFT descriptors D, stored as +% columns of the matrix D. D has the same format used by SIFT(). +% +% PLOTSIFTDESCRIPTOR(D,F) plots the SIFT descriptors warped to the +% SIFT frames F, specified as columns of the matrix F. F has +% the same format used by SIFT(). +% +% H=PLOTSIFTDESCRIPTOR(...) returns the handle H to the line drawing +% representing the descriptors. +% +% REMARK. Currently the function supports only descriptors with 4x4 +% spatial bins and 8 orientation bins (Lowe's default.) +% +% See also PLOTSIFTFRAME(), PLOTMATCHES(), PLOTSS(). + +% AUTORIGHTS +% Copyright (c) 2006 The Regents of the University of California. +% All Rights Reserved. +% +% Created by Andrea Vedaldi +% UCLA Vision Lab - Department of Computer Science +% +% Permission to use, copy, modify, and distribute this software and its +% documentation for educational, research and non-profit purposes, +% without fee, and without a written agreement is hereby granted, +% provided that the above copyright notice, this paragraph and the +% following three paragraphs appear in all copies. +% +% This software program and documentation are copyrighted by The Regents +% of the University of California. The software program and +% documentation are supplied "as is", without any accompanying services +% from The Regents. The Regents does not warrant that the operation of +% the program will be uninterrupted or error-free. The end-user +% understands that the program was developed for research purposes and +% is advised not to rely exclusively on the program for any reason. +% +% This software embodies a method for which the following patent has +% been issued: "Method and apparatus for identifying scale invariant +% features in an image and use of same for locating an object in an +% image," David G. Lowe, US Patent 6,711,293 (March 23, +% 2004). Provisional application filed March 8, 1999. Asignee: The +% University of British Columbia. +% +% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. + +lowe_compatible = 1 ; + +% -------------------------------------------------------------------- +% Check the arguments +% -------------------------------------------------------------------- + +if(size(d,1) ~= 128) + error('D should be a 128xK matrix (only standard descriptors accepted)') ; +end + +if nargin > 1 + if(size(f,1) ~= 4) + error('F should be a 4xK matrix'); + end + + if(size(f,2) ~= size(f,2)) + error('D and F must have the same number of columns') ; + end +end + +% Descriptors are often non-double numeric arrays +d = double(d) ; +K = size(d,2) ; +if nargin < 2 + f = repmat([0;0;1;0],1,K) ; +end + +maginf = 3.0 ; +NBP=4 ; +NBO=8 ; + +% -------------------------------------------------------------------- +% Do the job +% -------------------------------------------------------------------- + +xall=[] ; +yall=[] ; + +for k=1:K + SBP = maginf * f(3,k) ; + th=f(4,k) ; + c=cos(th) ; + s=sin(th) ; + + [x,y] = render_descr(d(:,k)) ; + xall = [xall SBP*(c*x-s*y)+f(1,k)] ; + yall = [yall SBP*(s*x+c*y)+f(2,k)] ; +end + +h=line(xall,yall) ; + +% -------------------------------------------------------------------- +% Helper functions +% -------------------------------------------------------------------- + +% Renders a single descriptor +function [x,y] = render_descr( d ) + +lowe_compatible=1; +NBP=4 ; +NBO=8 ; + +[x,y] = meshgrid(-NBP/2:NBP/2,-NBP/2:NBP/2) ; + +% Rescale d so that the biggest peak fits inside the bin diagram +d = 0.4 * d / max(d(:)) ; + +% We have NBP*NBP bins to plot. Here are the centers: +xc = x(1:end-1,1:end-1) + 0.5 ; +yc = y(1:end-1,1:end-1) + 0.5 ; + +% We swap the order of the bin diagrams because they are stored row +% major into the descriptor (Lowe's convention that we follow.) +xc = xc' ; +yc = yc' ; + +% Each bin contains a star with eight tips +xc = repmat(xc(:)',NBO,1) ; +yc = repmat(yc(:)',NBO,1) ; + +% Do the stars +th=linspace(0,2*pi,NBO+1) ; +th=th(1:end-1) ; +if lowe_compatible + xd = repmat(cos(-th), 1, NBP*NBP ) ; + yd = repmat(sin(-th), 1, NBP*NBP ) ; +else + xd = repmat(cos(th), 1, NBP*NBP ) ; + yd = repmat(sin(th), 1, NBP*NBP ) ; +end +xd = xd .* d(:)' ; +yd = yd .* d(:)' ; + +% Re-arrange in sequential order the lines to draw +nans = NaN * ones(1,NBP^2*NBO) ; +x1 = xc(:)' ; +y1 = yc(:)' ; +x2 = x1 + xd ; +y2 = y1 + yd ; +xstars = [x1;x2;nans] ; +ystars = [y1;y2;nans] ; + +% Horizontal lines of the grid +nans = NaN * ones(1,NBP+1); +xh = [x(:,1)' ; x(:,end)' ; nans] ; +yh = [y(:,1)' ; y(:,end)' ; nans] ; + +% Verical lines of the grid +xv = [x(1,:) ; x(end,:) ; nans] ; +yv = [y(1,:) ; y(end,:) ; nans] ; + +x=[xstars(:)' xh(:)' xv(:)'] ; +y=[ystars(:)' yh(:)' yv(:)'] ; diff --git a/SD-VBS/benchmarks/sift/src/matlab/plotsiftframe.m b/SD-VBS/benchmarks/sift/src/matlab/plotsiftframe.m new file mode 100644 index 0000000..984dfc9 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/plotsiftframe.m @@ -0,0 +1,119 @@ +function h=plotsiftframe(frames,labels) +% PLOTSIFTFRAME Plot SIFT frame +% H=PLOTSIFTFRAME(FRAMES) plots the SIFT frames FRAMES and returns +% and handle H to the resulting line set. FRAMES has the same format +% used by SIFT(). +% +% PLOTSIFTFRAME(FRAMES, LABELS) displays nearby the frame centers +% the indexes specified by the vector LABELS. This operation is slow +% for large sets of frames. +% +% A SIFT frame is denoted by a circle, representing its support, and +% one of its radii, representing its orientation. The support is a +% disk with radius equal to six times the scale SIGMA of the +% frame. If the standard parameters are used for the detector, this +% corresponds to four times the standard deviation of the Gaussian +% window that has been uses to estimate the orientation, which is in +% fact equal to 1.5 times the scale SIGMA. +% +% This function is considerably more efficient if called once on a +% whole set of frames as opposed to multiple times, one for each +% frame. +% +% See also PLOTMATCHES(), PLOTSIFTDESCRIPTOR(), PLOTSS(). + +% AUTORIGHTS +% Copyright (c) 2006 The Regents of the University of California. +% All Rights Reserved. +% +% Created by Andrea Vedaldi +% UCLA Vision Lab - Department of Computer Science +% +% Permission to use, copy, modify, and distribute this software and its +% documentation for educational, research and non-profit purposes, +% without fee, and without a written agreement is hereby granted, +% provided that the above copyright notice, this paragraph and the +% following three paragraphs appear in all copies. +% +% This software program and documentation are copyrighted by The Regents +% of the University of California. The software program and +% documentation are supplied "as is", without any accompanying services +% from The Regents. The Regents does not warrant that the operation of +% the program will be uninterrupted or error-free. The end-user +% understands that the program was developed for research purposes and +% is advised not to rely exclusively on the program for any reason. +% +% This software embodies a method for which the following patent has +% been issued: "Method and apparatus for identifying scale invariant +% features in an image and use of same for locating an object in an +% image," David G. Lowe, US Patent 6,711,293 (March 23, +% 2004). Provisional application filed March 8, 1999. Asignee: The +% University of British Columbia. +% +% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. + +% -------------------------------------------------------------------- +% Check the arguments +% -------------------------------------------------------------------- + +if size(frames,1) ~= 4 + error('FRAMES should be a 4xK matrix') ; +end + +K = size(frames,2) ; + +if nargin > 1 + putlabels = 1 ; +end + +% -------------------------------------------------------------------- +% Do the work +% -------------------------------------------------------------------- + +hold on ; +K=size(frames,2) ; +thr=linspace(0,2*pi,40) ; + +allx = nan*ones(1, 40*K+(K-1)) ; +ally = nan*ones(1, 40*K+(K-1)) ; + +allxf = nan*ones(1, 3*K) ; +allyf = nan*ones(1, 3*K) ; + +putlabel=0 ; + +for k=1:K + xc=frames(1,k) ; + yc=frames(2,k) ; + r=1.5*4*frames(3,k) ; + th=frames(4,k) ; + + x = r*cos(thr) + xc ; + y = r*sin(thr) + yc ; + + allx((k-1)*(41) + (1:40)) = x ; + ally((k-1)*(41) + (1:40)) = y ; + + allxf((k-1)*3 + (1:2)) = [xc xc+r*cos(th)] ; + allyf((k-1)*3 + (1:2)) = [yc yc+r*sin(th)] ; + + if putlabel + x=xc+r ; + y=yc ; + h=text(x+2,y,sprintf('%d',labels(k))) ; + set(h,'Color',[1 0 0]) ; + plot(x,y,'r.') ; + end + +end + +h=line([allx nan allxf], [ally nan allyf], 'Color','g','LineWidth',3) ; diff --git a/SD-VBS/benchmarks/sift/src/matlab/plotss.m b/SD-VBS/benchmarks/sift/src/matlab/plotss.m new file mode 100644 index 0000000..17868b6 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/plotss.m @@ -0,0 +1,68 @@ +function plotss(ss,field) +% PLOTSS Plot scale space +% PLOTSS(SS) plots all octaves of the scale space SS. +% +% See also GAUSSIANSS(), DIFFSS(). + +% AUTORIGHTS +% Copyright (c) 2006 The Regents of the University of California. +% All Rights Reserved. +% +% Created by Andrea Vedaldi +% UCLA Vision Lab - Department of Computer Science +% +% Permission to use, copy, modify, and distribute this software and its +% documentation for educational, research and non-profit purposes, +% without fee, and without a written agreement is hereby granted, +% provided that the above copyright notice, this paragraph and the +% following three paragraphs appear in all copies. +% +% This software program and documentation are copyrighted by The Regents +% of the University of California. The software program and +% documentation are supplied "as is", without any accompanying services +% from The Regents. The Regents does not warrant that the operation of +% the program will be uninterrupted or error-free. The end-user +% understands that the program was developed for research purposes and +% is advised not to rely exclusively on the program for any reason. +% +% This software embodies a method for which the following patent has +% been issued: "Method and apparatus for identifying scale invariant +% features in an image and use of same for locating an object in an +% image," David G. Lowe, US Patent 6,711,293 (March 23, +% 2004). Provisional application filed March 8, 1999. Asignee: The +% University of British Columbia. +% +% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. + +if nargin > 2 + error('Too many arguments.') ; +end + +omin = ss.omin ; +smin = ss.smin ; +nlevels = ss.smax-ss.smin+1 ; + +for oi=1:ss.O + for si=1:nlevels + tightsubplot(nlevels, ss.O, nlevels*(oi-1)+si) ; + s = si-1 + smin ; + o = oi-1 + omin ; + sigma = ss.sigma0 * 2^(s/ss.S + o) ; + F=squeeze(ss.octave{oi}(:,:,si)) ; + [M,N]=size(F) ; + imagesc(squeeze(ss.octave{oi}(:,:,si))) ; axis image ; axis off ; + h=text(M/10,N/20,sprintf('(o,s)=(%d,%d), sigma=%f',o,s,sigma)) ; + set(h,'BackgroundColor','w','Color','k') ; + end +end + + diff --git a/SD-VBS/benchmarks/sift/src/matlab/script_run_profile.m b/SD-VBS/benchmarks/sift/src/matlab/script_run_profile.m new file mode 100644 index 0000000..c157507 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/script_run_profile.m @@ -0,0 +1,23 @@ +function script_run_profile(dataDir, resultDir, type, common, tooldir) + +path(path,common); +sift_compile; + +I1=readImage([dataDir, '/1.bmp']) ; + +rows = size(I1,1); +cols = size(I1,2); +fprintf(1,'Input size\t\t- (%dx%d)\n', rows, cols); + +I1=I1-min(I1(:)) ; +I1=I1/max(I1(:)) ; + +%% Timing +start = photonStartTiming; +frames1 = sift( I1) ; +stop = photonEndTiming; +elapsed = photonReportTiming(start, stop); +photonPrintTiming(elapsed); + +fWriteMatrix(frames1, dataDir); + diff --git a/SD-VBS/benchmarks/sift/src/matlab/sift.m b/SD-VBS/benchmarks/sift/src/matlab/sift.m new file mode 100644 index 0000000..ee79371 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/sift.m @@ -0,0 +1,110 @@ +% Input: +% I the input image, with pixel values normalize to lie betwen [0,1]. +% +% Output: +% features a structure which contains the following fields: +% pos an n*2 matrix containing the (x,y) coordinates of the keypoints +% stored in rows. +% scale an n*3 matrix with rows describing the scale of each keypoint +% (i.e., first column specifies the octave, second column specifies the interval, and +% third column specifies sigma). +% orient a n*1 vector containing the orientations of the keypoints [-pi,pi). +% desc an n*128 matrix with rows containing the feature descriptors +% corresponding to the keypoints. + +function [frames]=sift(I) + +[M,N,C] = size(I) ; + +% Lowe's choices +S=3 ; +omin=-1 ; +O=floor(log2(min(M,N)))-omin-4 ; % Up to 16x16 images +sigma0=1.6*2^(1/S) ; +sigman=0.5 ; +thresh = 0.04 / S / 2 ; +r = 10 ; + +NBP = 4 ; +NBO = 8 ; +magnif = 3.0 ; + +% Parese input +compute_descriptor = 0 ; +discard_boundary_points = 1 ; +verb = 0 ; + +smin = -1; +smax = S+1; +intervals = smax - smin + 1; + + +% -------------------------------------------------------------------- +% Parameters +% -------------------------------------------------------------------- + +oframes = [] ; +frames = [] ; +descriptors = [] ; + +% -------------------------------------------------------------------- +% SIFT Detector and Descriptor +% -------------------------------------------------------------------- + +% Compute scale spaces +%if verb > 0, fprintf('SIFT: computing scale space...') ; end + +gss = gaussianss(I,sigman,O,S,omin,-1,S+1,sigma0) ; +dogss = diffss(gss) ; +frames = []; + +%% To maintain consistency with C code. Once C code is ready, this will be uncommented. +for o=1:O %for o=1:O + + % Local maxima of the DOG octave + % The 80% tricks discards early very weak points before refinement. + + idx = siftlocalmax( dogss.octave{o}, 0.8*thresh ) ; + idx = [idx , siftlocalmax( - dogss.octave{o}, 0.8*thresh)] ; + + K=length(idx) ; + [i,j,s] = ind2sub( size( dogss.octave{o} ), idx ) ; + + y=i-1 ; + x=j-1 ; + s=s-1+dogss.smin ; + oframes = [x(:)';y(:)';s(:)'] ; + +% fWriteMatrix(oframes, '../data/sim'); + + % Remove points too close to the boundary + if discard_boundary_points + oframes = filter_boundary_points(size(dogss.octave{o}), oframes ) ; + end + + % Refine the location, threshold strength and remove points on edges + oframes = siftrefinemx(... + oframes, ... + dogss.octave{o}, ... + dogss.smin, ... + thresh, ... + r) ; + + frames = [frames, oframes]; + +end +end + +%% -------------------------------------------------------------------- +%% Helpers +%% -------------------------------------------------------------------- +function oframes=filter_boundary_points(sz, oframes) + +sel=find(... + oframes(1,:) > 3 & ... + oframes(1,:) < sz(2)-3 & ... + oframes(2,:) > 3 & ... + oframes(2,:) < sz(1)-3 ) ; + +oframes=oframes(:,sel) ; +end diff --git a/SD-VBS/benchmarks/sift/src/matlab/sift_compile.m b/SD-VBS/benchmarks/sift/src/matlab/sift_compile.m new file mode 100644 index 0000000..78c0034 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/sift_compile.m @@ -0,0 +1,60 @@ +function sift_compile(type) +% SIFT_COMPILE Compile MEX files +% Compiling under Windows requires at least Visual C 6 or LCC. You +% might try other compilers, but most likely you will need to edit +% this file. + +siftroot = fileparts(which('siftcompile')) ; +opts = { '-O', '-I.' } ; +%opts = { opts{:}, '-v' } ; + +if nargin < 1 + type = 'visualc' ; +end + +switch computer + case 'PCWIN' + warning('NOTE: compiling has been tested only with Visual C 6-7 and LCC') ; + switch type + case 'visualc' + lib{1}=[matlabroot '\extern\lib\win32\microsoft\libmwlapack.lib'] ; + lib{2}=[matlabroot '\extern\lib\win32\microsoft\msvc60\libmwlapack.lib']; + lib{3}=[matlabroot '\extern\lib\win32\microsoft\msvc71\libmwlapack.lib']; + case 'lcc' + lib{1}=[matlabroot '\extern\lib\win32\lcc\libmwlapack.lib'] ; + end + found=0; + for k=1:length(lib) + fprintf('Trying LAPACK lib ''%s''\n',lib{k}) ; + found=exist(lib{k}) ; + if found ~= 0 + break ; + end + end + if found == 0 + error('Could not find LAPACK library. Please edit this M-file to fix the issue.'); + end + opts = {opts{:}, '-DWINDOWS'} ; + opts = {opts{:}, lib{k}} ; + + case 'MAC' + opts = {opts{:}, '-DMACOSX'} ; + opts = {opts{:}, 'CFLAGS=\$CFLAGS -faltivec'} ; + + case 'GLNX86' + opts = {opts{:}, '-DLINUX' } ; + + otherwise + error(['Unsupported architecture ', computer, '. Please edit this M-mfile to fix the issue.']) ; +end + +mex('imsmooth.c',opts{:}) ; +mex('siftlocalmax.c',opts{:}) ; +mex('siftrefinemx.c',opts{:}) ; +mex('siftormx.c',opts{:}) ; +mex('siftdescriptor.c',opts{:}) ; +mex('siftmatch.c',opts{:}) ; + + + + diff --git a/SD-VBS/benchmarks/sift/src/matlab/sift_demo2.m b/SD-VBS/benchmarks/sift/src/matlab/sift_demo2.m new file mode 100644 index 0000000..a6f0c66 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/sift_demo2.m @@ -0,0 +1,110 @@ +cd % SIFT_DEMO2 Demonstrate SIFT code (2) +% This is similar to SIFT_DEMO(). +% +% See also SIFT_DEMO(). + +% AUTORIGHTS +% Copyright (c) 2006 The Regents of the University of California. +% All Rights Reserved. +% +% Created by Andrea Vedaldi +% UCLA Vision Lab - Department of Computer Science +% +% Permission to use, copy, modify, and distribute this software and its +% documentation for educational, research and non-profit purposes, +% without fee, and without a written agreement is hereby granted, +% provided that the above copyright notice, this paragraph and the +% following three paragraphs appear in all copies. +% +% This software program and documentation are copyrighted by The Regents +% of the University of California. The software program and +% documentation are supplied "as is", without any accompanying services +% from The Regents. The Regents does not warrant that the operation of +% the program will be uninterrupted or error-free. The end-user +% understands that the program was developed for research purposes and +% is advised not to rely exclusively on the program for any reason. +% +% This software embodies a method for which the following patent has +% been issued: "Method and apparatus for identifying scale invariant +% features in an image and use of same for locating an object in an +% image," David G. Lowe, US Patent 6,711,293 (March 23, +% 2004). Provisional application filed March 8, 1999. Asignee: The +% University of British Columbia. +% +% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. + +I1=imreadbw('data/landscape-a.jpg') ; % I1=I1(1:2:end,:) ; +I2=imreadbw('data/landscape-b.jpg') ; % I2=I2(1:2:end,:) ; +I1c=double(imread('data/landscape-a.jpg'))/255.0 ; +I2c=double(imread('data/landscape-b.jpg'))/255.0 ; + +I1=imsmooth(I1,.1) ; +I2=imsmooth(I2,.1) ; + +I1=I1-min(I1(:)) ; +I1=I1/max(I1(:)) ; +I2=I2-min(I2(:)) ; +I2=I2/max(I2(:)) ; + +S=3 ; + +fprintf('Computing frames and descriptors.\n') ; +[frames1,descr1,gss1,dogss1] = sift( I1, 'Verbosity', 1, 'Threshold', ... + 0.005, 'NumLevels', S ) ; +[frames2,descr2,gss2,dogss2] = sift( I2, 'Verbosity', 1, 'Threshold', ... + 0.005, 'NumLevels', S ) ; + +figure(11) ; clf ; plotss(dogss1) ; colormap gray ; +figure(12) ; clf ; plotss(dogss2) ; colormap gray ; +drawnow ; + +figure(2) ; clf ; +tightsubplot(1,2,1) ; imagesc(I1) ; colormap gray ; axis image ; +hold on ; +h=plotsiftframe( frames1 ) ; set(h,'LineWidth',2,'Color','g') ; +h=plotsiftframe( frames1 ) ; set(h,'LineWidth',1,'Color','k') ; + +tightsubplot(1,2,2) ; imagesc(I2) ; colormap gray ; axis image ; +hold on ; +h=plotsiftframe( frames2 ) ; set(h,'LineWidth',2,'Color','g') ; +h=plotsiftframe( frames2 ) ; set(h,'LineWidth',1,'Color','k') ; + +fprintf('Computing matches.\n') ; +% By passing to integers we greatly enhance the matching speed (we use +% the scale factor 512 as Lowe's, but it could be greater without +% overflow) +descr1=uint8(512*descr1) ; +descr2=uint8(512*descr2) ; +tic ; +matches=siftmatch( descr1, descr2, 3 ) ; +fprintf('Matched in %.3f s\n', toc) ; + +figure(3) ; clf ; +plotmatches(I1c,I2c,frames1(1:2,:),frames2(1:2,:),matches,... + 'Stacking','v') ; +drawnow ; + +% Movie +figure(4) ; set(gcf,'Position',[10 10 1024 512]) ; +figure(4) ; clf ; +tightsubplot(1,1); +imagesc(I1) ; colormap gray ; axis image ; hold on ; +h=plotsiftframe( frames1 ) ; set(h,'LineWidth',1,'Color','g') ; +h=plot(frames1(1,:),frames1(2,:),'r.') ; +MOV(1)=getframe ; + +figure(4) ; clf ; +tightsubplot(1,1); +imagesc(I2) ; colormap gray ; axis image ; hold on ; +h=plotsiftframe( frames2 ) ; set(h,'LineWidth',1,'Color','g') ; +h=plot(frames2(1,:),frames2(2,:),'r.') ; +MOV(2)=getframe ; diff --git a/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.css b/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.css new file mode 100644 index 0000000..b9bd8e3 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.css @@ -0,0 +1,159 @@ +/** file: default.css + ** author: Andrea Vedaldi + ** description: Default CSS sylesheet for SIFT_GENDOC.PL + **/ + +/* AUTORIGHTS */ + +h1 +{ + color: #3a3a3a ; +} + +pre +{ + font-family: monaco, courier, monospace ; + font-size: 14px ; + color: #3a3a3a ; +} + +body +{ + margin: 10px ; + padding: 1em ; + right: 0px ; + font-family: arial, sans-serif ; + + color: black ; + background-color: #fffffa ; +} + +/* ------------------------------------------------------------------ + * Module + * --------------------------------------------------------------- */ + +div.module +{ + margin-bottom: 1em ; + //background-color: #ffffff ; +} + +div.module h1 +{ + margin: 0px ; + margin-bottom: 0.5em ; + padding: 0px ; + font-size: 1.5em ; + border-bottom: 2px #3a3a3a solid ; +} + +div.module pre +{ + padding-left: 1em ; +} + +div.module div.footer +{ + clear: both ; +} + +div.module div.footer :link +{ + color: white ; + text-decoration: none ; +} + +/* ------------------------------------------------------------------ + * Module index + * --------------------------------------------------------------- */ + +div.module div.index +{ + font-family: sans-serif ; + font-size: 0.8em ; + + width: 15em ; + float: right ; + border: 1px #ccc0c0 solid ; + background-color: #fcf0f0 ; +} + +div.module div.index h1 +{ + font-size: 1em ; + font-weight: bold ; + text-align: center ; + border: none ; + padding: 0 ; + margin: 0 ; +} + +div.module div.index ul +{ + list-style-type: square ; + list-style-position: inside ; + color: #3a3a3a ; + padding: 0 ; + margin: 0 ; + padding: 0.3em ; +} + +div.module div.index ul li +{ + margin-bottom: 0.1em ; +} + + +/* ------------------------------------------------------------------ + * Mfile + * --------------------------------------------------------------- */ + +div.mfile +{ + background-color: white ; + margin-bottom: 1em ; + border: 1px #aaaaaa solid ; +} + +div.mfile h1 +{ + margin: 0 ; + padding: 0 ; + padding-left: 10px ; + font-size: 1.3em ; + background-color: #f0f0ff ; +} + +div.mfile h1 span.name +{ + font-family: monaco, courier, fixed ; +} + +div.mfile h1 span.brief +{ + padding-left: 10px ; + font-size: 0.9em ; + font-weight: normal ; +} + +div.mfile pre +{ + padding-left: 1em ; +} + +div.mfile div.footer +{ + font-size: 0.8em ; + padding: 0.3em ; +} + +div.mfile div.footer a +{ + text-decoration: none ; + color: #777777 ; +} + +div.mfile div.footer a:hover +{ + text-decoration: underline ; +} diff --git a/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.m b/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.m new file mode 100644 index 0000000..417f9c5 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.m @@ -0,0 +1,50 @@ +function sift_gendoc +% SIFT_GENDOC Generate documentation +% This function extracts the documentation from the MEX files and +% creates an HTML manual. + +% AUTORIGHTS +% Copyright (c) 2006 The Regents of the University of California. +% All Rights Reserved. +% +% Created by Andrea Vedaldi +% UCLA Vision Lab - Department of Computer Science +% +% Permission to use, copy, modify, and distribute this software and its +% documentation for educational, research and non-profit purposes, +% without fee, and without a written agreement is hereby granted, +% provided that the above copyright notice, this paragraph and the +% following three paragraphs appear in all copies. +% +% This software program and documentation are copyrighted by The Regents +% of the University of California. The software program and +% documentation are supplied "as is", without any accompanying services +% from The Regents. The Regents does not warrant that the operation of +% the program will be uninterrupted or error-free. The end-user +% understands that the program was developed for research purposes and +% is advised not to rely exclusively on the program for any reason. +% +% This software embodies a method for which the following patent has +% been issued: "Method and apparatus for identifying scale invariant +% features in an image and use of same for locating an object in an +% image," David G. Lowe, US Patent 6,711,293 (March 23, +% 2004). Provisional application filed March 8, 1999. Asignee: The +% University of British Columbia. +% +% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. + +rootdir = fileparts(which('siftgendoc')) ; +d=pwd; +cd([rootdir]) ; +res = perl([rootdir '/siftgendoc.pl'],[rootdir '/dom/index.html']) ; +cd(d) ; +fprintf(res) ; diff --git a/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.pl b/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.pl new file mode 100755 index 0000000..6bc0b5e --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/sift_gendoc.pl @@ -0,0 +1,296 @@ +#!/usr/bin/perl -w +## file: sift_gendoc.pl +## author: Andrea Vedaldi +## description: Summarize MATLAB M-Files docs + +# AUTORIGHTS +# Copyright (c) 2006 The Regents of the University of California. +# All Rights Reserved. +# +# Created by Andrea Vedaldi +# UCLA Vision Lab - Department of Computer Science +# +# Permission to use, copy, modify, and distribute this software and its +# documentation for educational, research and non-profit purposes, +# without fee, and without a written agreement is hereby granted, +# provided that the above copyright notice, this paragraph and the +# following three paragraphs appear in all copies. +# +# This software program and documentation are copyrighted by The Regents +# of the University of California. The software program and +# documentation are supplied "as is", without any accompanying services +# from The Regents. The Regents does not warrant that the operation of +# the program will be uninterrupted or error-free. The end-user +# understands that the program was developed for research purposes and +# is advised not to rely exclusively on the program for any reason. +# +# This software embodies a method for which the following patent has +# been issued: "Method and apparatus for identifying scale invariant +# features in an image and use of same for locating an object in an +# image," David G. Lowe, US Patent 6,711,293 (March 23, +# 2004). Provisional application filed March 8, 1999. Asignee: The +# University of British Columbia. +# +# IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +# FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +# INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +# ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +# CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +# A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +# BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +# MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. + +# Debugging level +$verb = 0 ; + +# PDF document location +$pdfdoc = 'sift.pdf' ; + +# This is the queue of directories to process. +# It gets filled recursively. +@dir_fifo = ('.') ; + +# This will hold an entry for each m-file found +%mfiles = () ; + +# This will hold an entry for each module found +%modules = () ; + +# #ARGV is the index of the last element, which is 0 +if ($#ARGV == 0) { + open(FOUT,">$ARGV[0]") ; + print STDERR "Writing to file '$ARGV[0]'.\n" if $verb ; +} else { + *FOUT= *STDOUT; + print STDERR "Using standard output.\n" if $verb ; +} + +# Each module is a subdirectory. The subdirectory is used +# as the module ID. + +while ($module_path = shift(@dir_fifo)) { + print STDERR "=> '$module_path'\n" if $verb ; + + # get a first version of the module name + $module_path =~ m/.*\/([^\\]+)$/ ; + $module_name = $1 ; + + # start a new module + $module = { + 'id' => $module_path, + 'path' => $module_path, + 'name' => $module_name, + 'mfiles' => [], + 'description' => "" + } ; + + # ................................................................. + opendir(DIRHANDLE, $module->{'path'}) ; + FILE: foreach (sort readdir(DIRHANDLE)) { + $name = $_ ; + $path = $module->{'path'} . "/" . $_ ; + + # push if a directory + $_ = $path ; + if (-d) { + next FILE if /(\.$|\.\.$)/ ; + push( @dir_fifo, "$_" ) ; + next FILE ; + } + + # parse if .m and not test_ + next FILE unless (/.+\.m$/) ; + next FILE if (/test_.x*/) ; + $name =~ m/(.*).m/ ; + $name = $1 ; + print STDERR " . m-file '$name'\n" if $verb ; + my ($brief,$description) = get_comment($path) ; + + # topic description? + if (/overview/) { + print STDERR " * module description\n" if $verb ; + $module->{'id'} = $name ; + $module->{'name'} = $brief ; + $module->{'description'} = $description ; + next FILE ; + } + + # use names as IDs + $id = $name ; + $id =~ tr/A-Z/a-z/ ; + + # create a new mfile object + $mfile = { + 'id' => $id, + 'path' => $path, + 'name' => $name, + 'brief' => $brief, + 'module' => $module, + 'description' => $description + } ; + + # add a reference to this mfile + # object to the global mfile list + $mfiles{$id} = $mfile ; + + # add a reference to this mfile + # object to the current module + push( @{$module->{'mfiles'}}, $mfile) ; + } + closedir(DIRHANDLE) ; + # ................................................................ + + # add a reference to the current module to the global + # module list + $modules{$module->{'id'}} = $module ; +} + +# .................................................................... +# write documentation +# .................................................................... + +print FOUT < + + + + +EOF + +# sort modules by path +sub criterion { $modules{$a}->{'path'} cmp $modules{$b}->{'path'} ; } + +MODULE: +foreach $id ( sort criterion keys %modules ) { + my $module = $modules{$id} ; + my $rich_description = make_rich($module->{'description'}) ; + + next MODULE if $#{$module->{'mfiles'}} < 0 and length($rich_description) == 0; + + print FOUT < +

$module->{'name'}

+
+

Module contents

+
    +EOF + foreach( @{$module->{'mfiles'}} ) { + print FOUT "
  • " + . "$_->{'name'}
    • \n" ; + } + print FOUT < +
+ 
+      $rich_description
+
+
+
+ +EOF +} + +foreach $id (sort keys %mfiles) { + my $mfile = $mfiles{$id} ; + my $rich_description = make_rich($mfile->{'description'}) ; + + print FOUT < +

+ $mfile->{"name"} + $mfile->{"brief"} +

+ 
+$rich_description
+
+
+ + $mfile->{'module'}->{'name'} + +
+ +EOF +} + +print FOUT "" ; + +# Close file +close FOUT ; + + +# ------------------------------------------------------------------------- +sub get_comment { +# ------------------------------------------------------------------------- + local $_ ; + my $full_name = $_[0] ; + my @comment = () ; + + open IN,$full_name ; + SCAN: + while( ) { + next if /^function/ ; + last SCAN unless ( /^%/ ) ; + push(@comment, substr("$_",1)) ; + } + close IN ; + + my $brief = "" ; + if( $#comment >= 0 && $comment[0] =~ m/^\s*\w+\s+(.*)$/ ) { + $brief = $1 ; + splice (@comment, 0, 1) ; + } + + # from the first line + return ($brief, join("",@comment)) ; +} + +# ------------------------------------------------------------------------- +sub make_rich { +# ------------------------------------------------------------------------- + local $_ = $_[0] ; + s/([A-Z]+[A-Z0-9_]*)\(([^\)]*)\)/${\make_link($1,$2)}/g ; + s/PDF:([A-Z0-9_\-:.]+[A-Z0-9])/${\make_pdf_link($1)}/g ; + return $_ ; +} + +# ------------------------------------------------------------------------- +sub make_link { +# ------------------------------------------------------------------------- + local $_ ; + my $name = $_[0] ; + my $arg = $_[1] ; + my $id = $name ; + + # convert name to lower case and put into $_ + $id =~ tr/A-Z/a-z/ ; + + # get mfile + my $mfile = $mfiles{$id} ; + my $module = $modules{$id} ; + + # return as appropriate + if($mfile) { + return "" . $name . "" . "(" . $arg . ")" ; + } elsif($module) { + return "" . $name . + "" . "(" . $arg . ")" ; + } else { + return $name . "(" . $arg .")" ; + } +} + + +# ------------------------------------------------------------------------- +sub make_pdf_link { +# ------------------------------------------------------------------------- + local $_ ; + my $name = $_[0] ; + my $id = $name ; + + # convert name to lower case and put into $_ + $id =~ tr/A-Z/a-z/ ; + + return "PDF:$1" ; +} diff --git a/SD-VBS/benchmarks/sift/src/matlab/sift_overview.m b/SD-VBS/benchmarks/sift/src/matlab/sift_overview.m new file mode 100644 index 0000000..71557e4 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/sift_overview.m @@ -0,0 +1,33 @@ +% SIFT_OVERVIEW Scale-Invariant Feature Transfrom +% +% This is a MATLAB/C implementation of SIFT detector and descriptor +% [1]. You can: +% +% * Use SIFT() to detect the SIFT frames (keypoints) of a given image +% and compute their descriptors. Then you can use SIFTMATCH() to +% match the descriptors. +% +% * Use PLOTSS(), PLOTSIFTDESCRIPTOR(), PLOTSIFTFRAME(), +% PLOTMATCHES() to visualize the results. +% +% As SIFT is implemented by several reusable M and MEX files, you can +% also run portions of the algorithm, or change them. Specifically, +% you can: +% +% * Use SIFTDESCRIPTOR() to compute the SIFT descriptor from a list +% of frames and a scale space or plain image. +% +% * Use GAUSSIANSS() and DIFFSS() to compute the Gaussian and DOG +% scale spaces. +% +% * Use SIFTLOCALMAX(), SIFTREFINEMX(), SIFTORMX() to manually +% extract the SIFT frames from the DOG scale space. More in +% general, you can use SIFTLOCALMAX() to find maximizers of any +% multi-dimensional arrays. +% +% REFERENCES +% [1] D. G. Lowe, "Distinctive image features from scale-invariant +% keypoints," IJCV, vol. 2, no. 60, pp. 91 110, 2004. +% +% See also PDF:SIFT.INTRODUCTION. + diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.c b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.c new file mode 100644 index 0000000..63f4830 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.c @@ -0,0 +1,524 @@ +/* file: siftdescriptor +** author: Andrea Vedaldi +** description: Compute SIFT descriptors +**/ + +/* AUTORIGHTS +Copyright (c) 2006 The Regents of the University of California. +All Rights Reserved. + +Created by Andrea Vedaldi +UCLA Vision Lab - Department of Computer Science + +Permission to use, copy, modify, and distribute this software and its +documentation for educational, research and non-profit purposes, +without fee, and without a written agreement is hereby granted, +provided that the above copyright notice, this paragraph and the +following three paragraphs appear in all copies. + +This software program and documentation are copyrighted by The Regents +of the University of California. The software program and +documentation are supplied "as is", without any accompanying services +from The Regents. The Regents does not warrant that the operation of +the program will be uninterrupted or error-free. The end-user +understands that the program was developed for research purposes and +is advised not to rely exclusively on the program for any reason. + +This software embodies a method for which the following patent has +been issued: "Method and apparatus for identifying scale invariant +features in an image and use of same for locating an object in an +image," David G. Lowe, US Patent 6,711,293 (March 23, +2004). Provisional application filed March 8, 1999. Asignee: The +University of British Columbia. + +IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +*/ + +/* + REMARKS. The use of strcasecmp makes the function POSIX but not ANSI + compliant. When compling with Altivec, GCC Altivec extensions are + supported. +*/ + +#define LOWE_COMPATIBLE + +#include"mexutils.c" +#include +#include + +#ifdef WINDOWS +#include +#ifndef __cplusplus +#define sqrtf(x) ((float)sqrt((double)(x))) +#define powf(x,y) ((float)pow((double)(x),(double)(y))) +#define fabsf(x) ((float)fabs((double)(x))) +#define sinf(x) ((float)sin((double)(x))) +#define cosf(x) ((float)cos((double)(x))) +#define expf(x) ((float)exp((double)(x))) +#define atan2f(x,y) ((float)atan2((double)(x),(double)(y))) +#endif +#else +#include +#endif + +/* Altivec and Accelerate support. + * Very crude at this time. + */ +#if defined( MACOSX ) && defined( __ALTIVEC__ ) +#include +typedef union +{ + float x[4] ; + vFloat vec ; +} float4 ; +#endif + +#define greater(a,b) a > b +#define min(a,b) (((a)<(b))?(a):(b)) +#define max(a,b) (((a)>(b))?(a):(b)) + + +enum {SCALESPACE, NOSCALESPACE} ; + +enum {PROP_MAGNIF=0, + PROP_NBP, + PROP_NBO, + PROP_UNKNOWN} ; + +char const * properties [4] = + { "Magnif", + "NumSpatialBins", + "NumOrientBins", + 0L + } ; + +/** Fast fmodf for 2*PI + **/ +/*inline*/ +float fast_mod(float th) +{ + while(th < 0) th += 2*M_PI ; + while(th > 2*M_PI) th -= 2*M_PI ; + return th ; +} + +/** Fast floor. Equivalent to (int) floor(x) + **/ +/*inline*/ +int fast_floor(float x) +{ + return (int)( x - ((x>=0)?0:1) ) ; +} + +/** Normalizes in norm L_2 a descriptor. + **/ +void +normalize_histogram(float* L_begin, float* L_end) +{ + float* L_iter ; + float norm=0.0 ; + + for(L_iter = L_begin; L_iter != L_end ; ++L_iter) + norm += (*L_iter) * (*L_iter) ; + + norm = sqrtf(norm) ; + /* mexPrintf("%f\n",norm) ;*/ + + for(L_iter = L_begin; L_iter != L_end ; ++L_iter) + *L_iter /= norm ; +} + +/** @brief MATLAB Driver. + **/ +void +mexFunction(int nout, mxArray *out[], + int nin, const mxArray *in[]) +{ + int M,N,S=0,smin=0,K,num_levels=0 ; + const int* dimensions ; + const double* P_pt ; + const double* G_pt ; + float* descr_pt ; + float* buffer_pt ; + float sigma0 ; + float magnif = 3.0f ; /* Spatial bin extension factor. */ + int NBP = 4 ; /* Number of bins for one spatial direction (even). */ + int NBO = 8 ; /* Number of bins for the ortientation. */ + int mode = NOSCALESPACE ; + int buffer_size=0; + + enum {IN_G=0,IN_P,IN_SIGMA0,IN_S,IN_SMIN} ; + enum {OUT_L=0} ; + + /* ------------------------------------------------------------------ + ** Check the arguments + ** --------------------------------------------------------------- */ + + if (nin < 3) { + mexErrMsgTxt("At least three arguments are required") ; + } else if (nout > 1) { + mexErrMsgTxt("Too many output arguments."); + } + + if( !uIsRealScalar(in[IN_SIGMA0]) ) { + mexErrMsgTxt("SIGMA0 should be a real scalar") ; + } + + if(!mxIsDouble(in[IN_G]) || + mxGetNumberOfDimensions(in[IN_G]) > 3) { + mexErrMsgTxt("G should be a real matrix or 3-D array") ; + } + + sigma0 = (float) *mxGetPr(in[IN_SIGMA0]) ; + + dimensions = mxGetDimensions(in[IN_G]) ; + M = dimensions[0] ; + N = dimensions[1] ; + G_pt = mxGetPr(in[IN_G]) ; + + P_pt = mxGetPr(in[IN_P]) ; + K = mxGetN(in[IN_P]) ; + + if( !uIsRealMatrix(in[IN_P],-1,-1)) { + mexErrMsgTxt("P should be a real matrix") ; + } + + if ( mxGetM(in[IN_P]) == 4) { + /* Standard (scale space) mode */ + mode = SCALESPACE ; + num_levels = dimensions[2] ; + + if(nin < 5) { + mexErrMsgTxt("Five arguments are required in standard mode") ; + } + + if( !uIsRealScalar(in[IN_S]) ) { + mexErrMsgTxt("S should be a real scalar") ; + } + + if( !uIsRealScalar(in[IN_SMIN]) ) { + mexErrMsgTxt("SMIN should be a real scalar") ; + } + + if( !uIsRealMatrix(in[IN_P],4,-1)) { + mexErrMsgTxt("When the e mode P should be a 4xK matrix.") ; + } + + S = (int)(*mxGetPr(in[IN_S])) ; + smin = (int)(*mxGetPr(in[IN_SMIN])) ; + + } else if ( mxGetM(in[IN_P]) == 3 ) { + mode = NOSCALESPACE ; + num_levels = 1 ; + S = 1 ; + smin = 0 ; + } else { + mexErrMsgTxt("P should be either a 3xK or a 4xK matrix.") ; + } + + /* Parse the property-value pairs */ + { + char str [80] ; + int arg = (mode == SCALESPACE) ? IN_SMIN + 1 : IN_SIGMA0 + 1 ; + + while(arg < nin) { + int k ; + + if( !uIsString(in[arg],-1) ) { + mexErrMsgTxt("The first argument in a property-value pair" + " should be a string") ; + } + mxGetString(in[arg], str, 80) ; + +#ifdef WINDOWS + for(k = 0 ; properties[k] && strcmpi(str, properties[k]) ; ++k) ; +#else + for(k = 0 ; properties[k] && strcasecmp(str, properties[k]) ; ++k) ; +#endif + + switch (k) { + case PROP_NBP: + if( !uIsRealScalar(in[arg+1]) ) { + mexErrMsgTxt("'NumSpatialBins' should be real scalar") ; + } + NBP = (int) *mxGetPr(in[arg+1]) ; + if( NBP <= 0 || (NBP & 0x1) ) { + mexErrMsgTxt("'NumSpatialBins' must be positive and even") ; + } + break ; + + case PROP_NBO: + if( !uIsRealScalar(in[arg+1]) ) { + mexErrMsgTxt("'NumOrientBins' should be a real scalar") ; + } + NBO = (int) *mxGetPr(in[arg+1]) ; + if( NBO <= 0 ) { + mexErrMsgTxt("'NumOrientlBins' must be positive") ; + } + break ; + + case PROP_MAGNIF: + if( !uIsRealScalar(in[arg+1]) ) { + mexErrMsgTxt("'Magnif' should be a real scalar") ; + } + magnif = (float) *mxGetPr(in[arg+1]) ; + if( magnif <= 0 ) { + mexErrMsgTxt("'Magnif' must be positive") ; + } + break ; + + case PROP_UNKNOWN: + mexErrMsgTxt("Property unknown.") ; + break ; + } + arg += 2 ; + } + } + + /* ----------------------------------------------------------------- + * Pre-compute gradient and angles + * -------------------------------------------------------------- */ + /* Alloc two buffers and make sure their size is multiple of 128 for + * better alignment (used also by the Altivec code below.) + */ + buffer_size = (M*N*num_levels + 0x7f) & (~ 0x7f) ; + buffer_pt = (float*) mxMalloc( sizeof(float) * 2 * buffer_size ) ; + descr_pt = (float*) mxCalloc( NBP*NBP*NBO*K, sizeof(float) ) ; + + { + /* Offsets to move in the scale space. */ + const int yo = 1 ; + const int xo = M ; + const int so = M*N ; + int x,y,s ; + +#define at(x,y) (*(pt + (x)*xo + (y)*yo)) + + /* Compute the gradient */ + for(s = 0 ; s < num_levels ; ++s) { + const double* pt = G_pt + s*so ; + for(x = 1 ; x < N-1 ; ++x) { + for(y = 1 ; y < M-1 ; ++y) { + float Dx = 0.5 * ( at(x+1,y) - at(x-1,y) ) ; + float Dy = 0.5 * ( at(x,y+1) - at(x,y-1) ) ; + buffer_pt[(x*xo+y*yo+s*so) + 0 ] = Dx ; + buffer_pt[(x*xo+y*yo+s*so) + buffer_size] = Dy ; + } + } + } + + /* Compute angles and modules */ + { + float* pt = buffer_pt ; + int j = 0 ; + while (j < N*M*num_levels) { + +#if defined( MACOSX ) && defined( __ALTIVEC__ ) + if( ((unsigned int)pt & 0x7) == 0 && j+3 < N*M*num_levels ) { + /* If aligned to 128 bit and there are at least 4 pixels left */ + float4 a, b, c, d ; + a.vec = vec_ld(0,(vector float*)(pt )) ; + b.vec = vec_ld(0,(vector float*)(pt + buffer_size)) ; + c.vec = vatan2f(b.vec,a.vec) ; + a.x[0] = a.x[0]*a.x[0]+b.x[0]*b.x[0] ; + a.x[1] = a.x[1]*a.x[1]+b.x[1]*b.x[1] ; + a.x[2] = a.x[2]*a.x[2]+b.x[2]*b.x[2] ; + a.x[3] = a.x[3]*a.x[3]+b.x[3]*b.x[3] ; + d.vec = vsqrtf(a.vec) ; + vec_st(c.vec,0,(vector float*)(pt + buffer_size)) ; + vec_st(d.vec,0,(vector float*)(pt )) ; + j += 4 ; + pt += 4 ; + } else { +#endif + float Dx = *(pt ) ; + float Dy = *(pt + buffer_size) ; + *(pt ) = sqrtf(Dx*Dx + Dy*Dy) ; + *(pt + buffer_size) = atan2f(Dy, Dx) ; + j += 1 ; + pt += 1 ; +#if defined( MACOSX ) && defined( __ALTIVEC__ ) + } +#endif + + } + } + } + + /* ----------------------------------------------------------------- + * Do the job + * -------------------------------------------------------------- */ + if(K > 0) { + int p ; + + /* Offsets to move in the buffer */ + const int yo = 1 ; + const int xo = M ; + const int so = M*N ; + + /* Offsets to move in the descriptor. */ + /* Use Lowe's convention. */ + const int binto = 1 ; + const int binyo = NBO * NBP ; + const int binxo = NBO ; + const int bino = NBO * NBP * NBP ; + + for(p = 0 ; p < K ; ++p, descr_pt += bino) { + /* The SIFT descriptor is a three dimensional histogram of the position + * and orientation of the gradient. There are NBP bins for each spatial + * dimesions and NBO bins for the orientation dimesion, for a total of + * NBP x NBP x NBO bins. + * + * The support of each spatial bin has an extension of SBP = 3sigma + * pixels, where sigma is the scale of the keypoint. Thus all the bins + * together have a support SBP x NBP pixels wide . Since weighting and + * interpolation of pixel is used, another half bin is needed at both + * ends of the extension. Therefore, we need a square window of SBP x + * (NBP + 1) pixels. Finally, since the patch can be arbitrarly rotated, + * we need to consider a window 2W += sqrt(2) x SBP x (NBP + 1) pixels + * wide. + */ + const float x = (float) *P_pt++ ; + const float y = (float) *P_pt++ ; + const float s = (float) (mode == SCALESPACE) ? (*P_pt++) : 0.0 ; + const float theta0 = (float) *P_pt++ ; + + const float st0 = sinf(theta0) ; + const float ct0 = cosf(theta0) ; + const int xi = (int) floor(x+0.5) ; /* Round-off */ + const int yi = (int) floor(y+0.5) ; + const int si = (int) floor(s+0.5) - smin ; + const float sigma = sigma0 * powf(2, s / S) ; + const float SBP = magnif * sigma ; + const int W = (int) floor( sqrt(2.0) * SBP * (NBP + 1) / 2.0 + 0.5) ; + int bin ; + int dxi ; + int dyi ; + const float* pt ; + float* dpt ; + + /* Check that keypoints are within bounds . */ + + if(xi < 0 || + xi > N-1 || + yi < 0 || + yi > M-1 || + ((mode==SCALESPACE) && + (si < 0 || + si > dimensions[2]-1) ) ) + continue ; + + /* Center the scale space and the descriptor on the current keypoint. + * Note that dpt is pointing to the bin of center (SBP/2,SBP/2,0). + */ + pt = buffer_pt + xi*xo + yi*yo + si*so ; + dpt = descr_pt + (NBP/2) * binyo + (NBP/2) * binxo ; + +#define atd(dbinx,dbiny,dbint) (*(dpt + (dbint)*binto + (dbiny)*binyo + (dbinx)*binxo)) + + /* + * Process each pixel in the window and in the (1,1)-(M-1,N-1) + * rectangle. + */ + for(dxi = max(-W, 1-xi) ; dxi <= min(+W, N-2-xi) ; ++dxi) { + for(dyi = max(-W, 1-yi) ; dyi <= min(+W, M-2-yi) ; ++dyi) { + + /* Compute the gradient. */ + float mod = *(pt + dxi*xo + dyi*yo + 0 ) ; + float angle = *(pt + dxi*xo + dyi*yo + buffer_size ) ; +#ifdef LOWE_COMPATIBLE + float theta = fast_mod(-angle + theta0) ; +#else + float theta = fast_mod(angle - theta0) ; +#endif + /* Get the fractional displacement. */ + float dx = ((float)(xi+dxi)) - x; + float dy = ((float)(yi+dyi)) - y; + + /* Get the displacement normalized w.r.t. the keypoint orientation + * and extension. */ + float nx = ( ct0 * dx + st0 * dy) / SBP ; + float ny = (-st0 * dx + ct0 * dy) / SBP ; + float nt = NBO * theta / (2*M_PI) ; + + /* Get the gaussian weight of the sample. The gaussian window + * has a standard deviation of NBP/2. Note that dx and dy are in + * the normalized frame, so that -NBP/2 <= dx <= NBP/2. */ + const float wsigma = NBP/2 ; + float win = expf(-(nx*nx + ny*ny)/(2.0 * wsigma * wsigma)) ; + + /* The sample will be distributed in 8 adijacient bins. + * Now we get the ``lower-left'' bin. */ + int binx = fast_floor( nx - 0.5 ) ; + int biny = fast_floor( ny - 0.5 ) ; + int bint = fast_floor( nt ) ; + float rbinx = nx - (binx+0.5) ; + float rbiny = ny - (biny+0.5) ; + float rbint = nt - bint ; + int dbinx ; + int dbiny ; + int dbint ; + + /* Distribute the current sample into the 8 adijacient bins. */ + for(dbinx = 0 ; dbinx < 2 ; ++dbinx) { + for(dbiny = 0 ; dbiny < 2 ; ++dbiny) { + for(dbint = 0 ; dbint < 2 ; ++dbint) { + + if( binx+dbinx >= -(NBP/2) && + binx+dbinx < (NBP/2) && + biny+dbiny >= -(NBP/2) && + biny+dbiny < (NBP/2) ) { + float weight = win + * mod + * fabsf(1 - dbinx - rbinx) + * fabsf(1 - dbiny - rbiny) + * fabsf(1 - dbint - rbint) ; + + atd(binx+dbinx, biny+dbiny, (bint+dbint) % NBO) += weight ; + } + } + } + } + } + } + + { + /* Normalize the histogram to L2 unit length. */ + normalize_histogram(descr_pt, descr_pt + NBO*NBP*NBP) ; + + /* Truncate at 0.2. */ + for(bin = 0; bin < NBO*NBP*NBP ; ++bin) { + if (descr_pt[bin] > 0.2) descr_pt[bin] = 0.2; + } + + /* Normalize again. */ + normalize_histogram(descr_pt, descr_pt + NBO*NBP*NBP) ; + } + } + } + + /* Restore pointer to the beginning of the descriptors. */ + descr_pt -= NBO*NBP*NBP*K ; + + { + int k ; + double* L_pt ; + out[OUT_L] = mxCreateDoubleMatrix(NBP*NBP*NBO, K, mxREAL) ; + L_pt = mxGetPr(out[OUT_L]) ; + for(k = 0 ; k < NBP*NBP*NBO*K ; ++k) { + L_pt[k] = descr_pt[k] ; + } + } + + mxFree(descr_pt) ; + mxFree(buffer_pt) ; +} diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.m b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.m new file mode 100644 index 0000000..7764897 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.m @@ -0,0 +1,79 @@ +% SIFTDESCRIPTOR Compute SIFT descriptors +% DESCR = SIFTDESCRIPTOR(G, P, SIGMA0, S, MINS) computes the SIFT +% descriptors DESCR of the SIFT frames P defined on the octave G of +% a Gaussian scale space. SIGMA0, S and MINS are the the parameter +% of the scale space as explained in PDF:SIFT.USER.SS. P has one +% column per frame, specifiying the center X1,X2, the scal index s +% and the orientation THETA of the frame. +% +% As the function operates on a single octave, in order to process +% features spanning several octaves, one must group them and call +% SIFTDESCRIPTOR() once per each octave. +% +% DESCR = SIFTDESCRIPTOR(I, P, SIGMA0) operates on a plain image +% I. The image I is assumed to be pre-smoothed at scale SIGMA0 and P +% is a matrix with a column per frame, specifying the center X1,X2 +% and the orientation THETA (but NOT the scale). +% +% REMARK. The scale parameter s in P is the scale index, NOT the +% scale coordinate (see the PDF doc. for a discussion). +% +% Other parameters can be specfied as option-value paris. These +% are +% +% 'Magnif' [3.0] +% Frame magnification factor. Each spatial bin of the SIFT +% histogram has an exentsion equal to magnif * sigma, where +% magnif is the frame magnification factor and sigma is the scale +% of the frame. +% +% 'NumSpatialBins' [4] +% This parameter specifies the number of spatial bins in each +% spatial direction x1 and x2. It must be a positive and even +% number. +% +% 'NumOrientBins' [8] +% This parameter specifies the number of orietnation bins. It +% must be a positive number. +% +% See also SIFT(), GAUSSIANSS(), DIFFSS(), SIFTLOCALMAX(), +% PDF:SIFT.USER.DESCRIPTOR. + +% AUTORIGHTS +% Copyright (c) 2006 The Regents of the University of California. +% All Rights Reserved. +% +% Created by Andrea Vedaldi +% UCLA Vision Lab - Department of Computer Science +% +% Permission to use, copy, modify, and distribute this software and its +% documentation for educational, research and non-profit purposes, +% without fee, and without a written agreement is hereby granted, +% provided that the above copyright notice, this paragraph and the +% following three paragraphs appear in all copies. +% +% This software program and documentation are copyrighted by The Regents +% of the University of California. The software program and +% documentation are supplied "as is", without any accompanying services +% from The Regents. The Regents does not warrant that the operation of +% the program will be uninterrupted or error-free. The end-user +% understands that the program was developed for research purposes and +% is advised not to rely exclusively on the program for any reason. +% +% This software embodies a method for which the following patent has +% been issued: "Method and apparatus for identifying scale invariant +% features in an image and use of same for locating an object in an +% image," David G. Lowe, US Patent 6,711,293 (March 23, +% 2004). Provisional application filed March 8, 1999. Asignee: The +% University of British Columbia. +% +% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.mexa64 b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.mexa64 new file mode 100755 index 0000000..f094f56 Binary files /dev/null and b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.mexa64 differ diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.mexglx b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.mexglx new file mode 100755 index 0000000..4eb94d8 Binary files /dev/null and b/SD-VBS/benchmarks/sift/src/matlab/siftdescriptor.mexglx differ diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.c b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.c new file mode 100644 index 0000000..3646692 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.c @@ -0,0 +1,291 @@ +/* file: siftlocalmax.c +** author: Andrea Vedaldi +** description: Find local maximizer of multi-dimensional array. +**/ + +/* AUTORIGHTS +Copyright (c) 2006 The Regents of the University of California. +All Rights Reserved. + +Created by Andrea Vedaldi +UCLA Vision Lab - Department of Computer Science + +Permission to use, copy, modify, and distribute this software and its +documentation for educational, research and non-profit purposes, +without fee, and without a written agreement is hereby granted, +provided that the above copyright notice, this paragraph and the +following three paragraphs appear in all copies. + +This software program and documentation are copyrighted by The Regents +of the University of California. The software program and +documentation are supplied "as is", without any accompanying services +from The Regents. The Regents does not warrant that the operation of +the program will be uninterrupted or error-free. The end-user +understands that the program was developed for research purposes and +is advised not to rely exclusively on the program for any reason. + +This software embodies a method for which the following patent has +been issued: "Method and apparatus for identifying scale invariant +features in an image and use of same for locating an object in an +image," David G. Lowe, US Patent 6,711,293 (March 23, +2004). Provisional application filed March 8, 1999. Asignee: The +University of British Columbia. + +IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +*/ + +#include"mex.h" +#include +#include + +/** Matlab driver. + **/ +#define greater(a,b) ((a) > (b)+threshold) + +void +mexFunction(int nout, mxArray *out[], + int nin, const mxArray *in[]) +{ + int M, N ; + const double* F_pt ; + int ndims ; + int pdims = -1 ; + int* offsets ; + int* midx ; + int* neighbors ; + int nneighbors ; + int* dims ; + enum {F=0,THRESHOLD,P} ; + enum {MAXIMA=0} ; + double threshold = - mxGetInf() ; + + + /* ------------------------------------------------------------------ + * Check the arguments + * --------------------------------------------------------------- */ + if(nin > 1) { + if(!uIsRealScalar(in[THRESHOLD])) { + mexErrMsgTxt("THRESHOLD must be a real scalar.") ; + } + threshold = *mxGetPr(in[THRESHOLD]) ; + } + + ndims = mxGetNumberOfDimensions(in[F]) ; + +if(ndims !=3) + printf("Error Error: NDIMS not equal to 3. Not handled by C version\n"); + + { + /* We need to make a copy because in one special case (see below) + we need to adjust dims[]. + */ + int d ; + const int* const_dims = (int*) mxGetDimensions(in[F]) ; + dims = mxMalloc(sizeof(int)*ndims) ; + for(d=0 ; d < ndims ; ++d) + { +/* printf("Const Dimensions = %d\n", const_dims[d]); +*/ + dims[d] = const_dims[d] ; + } + } + M = dims[0] ; + N = dims[1] ; + F_pt = mxGetPr(in[F]) ; + + /* + If there are only two dimensions and if one is singleton, then + assume that a vector has been provided as input (and treat this + as a COLUMN matrix with p=1). We do this because Matlab does not + distinguish between vectors and 1xN or Mx1 matrices and because + the cases 1xN and Mx1 are trivial (the result is alway empty). + */ + if((ndims == 2) && (pdims < 0) && (M == 1 || N == 1)) { + + printf("ERROR ERROR: Entered a different loop here. Not handled by C version"); + + pdims = 1 ; + M = (M>N)?M:N ; + N = 1 ; + dims[0]=M ; + dims[1]=N ; + } + + /* search the local maxima along the first p dimensions only */ + if(pdims < 0) + { + pdims = ndims ; + } + + if(pdims > ndims) { + mxFree(dims) ; + mexErrMsgTxt("P must not be greater than the number of dimensions") ; + } + + /* ------------------------------------------------------------------ + * Do the job + * --------------------------------------------------------------- */ + { + int maxima_size = M*N ; + int* maxima_start = (int*) mxMalloc(sizeof(int) * maxima_size) ; + int* maxima_iterator = maxima_start ; + int* maxima_end = maxima_start + maxima_size ; + int i,j,h,o ; + const double* pt = F_pt ; + + /* Compute the offsets between dimensions. */ + offsets = (int*) mxMalloc(sizeof(int) * ndims) ; + offsets[0] = 1 ; + + for(h = 1 ; h < ndims ; ++h) + { + offsets[h] = offsets[h-1]*dims[h-1] ; +/* printf("%d:%d\t%d\n", h, offsets[h], dims[h-1]); +*/ + } + + /* Multi-index. */ + midx = (int*) mxMalloc(sizeof(int) * ndims) ; + for(h = 0 ; h < ndims ; ++h) + midx[h] = 1 ; + + /* Neighbors. */ + nneighbors = 1 ; + o=0 ; + for(h = 0 ; h < pdims ; ++h) { + nneighbors *= 3 ; + midx[h] = -1 ; + o -= offsets[h] ; + } + nneighbors -= 1 ; + neighbors = (int*) mxMalloc(sizeof(int) * nneighbors) ; + + /* Precompute offsets from offset(-1,...,-1,0,...0) to + * offset(+1,...,+1,0,...,0). */ + i = 0 ; + + while(true) { + if(o != 0) + neighbors[i++] = o ; + h = 0 ; + while( o += offsets[h], (++midx[h]) > 1 ) { + o -= 3*offsets[h] ; + midx[h] = -1 ; + if(++h >= pdims) + goto stop ; + } + } + stop: ; + + /* Starts at the corner (1,1,...,1,0,0,...0) */ + for(h = 0 ; h < pdims ; ++h) { + midx[h] = 1 ; + pt += offsets[h] ; +/* printf("%d:%x\t%d\n", h, pt, offsets[h]); +*/ + } + + for(h = pdims ; h < ndims ; ++h) { + midx[h] = 0 ; + } + + /* --------------------------------------------------------------- + * Loop + * ------------------------------------------------------------ */ + + /* + If any dimension in the first P is less than 3 elements wide + then just return the empty matrix (if we proceed without doing + anything we break the carry reporting algorithm below). + */ + for(h=0 ; h < pdims ; ++h) + if(dims[h] < 3) goto end ; + + while(true) { + + /* Propagate carry along multi index midx */ + h = 0 ; + while((midx[h]) >= dims[h] - 1) { +/* pt += 2*offsets[h] ; skip first and last el. */ + midx[h] = 1 ; + if(++h >= pdims) + goto next_layer ; + ++midx[h] ; + } + + /* Scan neighbors */ + { + double v = *pt ; + bool is_greater = (v >= threshold) ; + /* printf("%f\n", v); + */ + i = 0 ; + while(is_greater && i < nneighbors) + { + is_greater &= v > *(pt + neighbors[i++]) ; + } + + /* Add the local maximum */ + if(is_greater) { + /* Need more space? */ + if(maxima_iterator == maxima_end) { + maxima_size += M*N ; + maxima_start = (int*) mxRealloc(maxima_start, + maxima_size*sizeof(int)) ; + maxima_end = maxima_start + maxima_size ; + maxima_iterator = maxima_end - M*N ; + } + + *maxima_iterator++ = pt - F_pt + 1 ; + } + + /* Go to next element */ + pt += 1 ; + ++midx[0] ; + continue ; + + next_layer: ; + if( h >= ndims ) + goto end ; + + while((++midx[h]) >= dims[h]) { + midx[h] = 0 ; + if(++h >= ndims) + goto end ; + } + } + } + end:; + /* Return. */ + { + double* M_pt ; +/* printf("%x\t%x\t%d\n", maxima_iterator, maxima_start, maxima_iterator-maxima_start); +*/ + out[MAXIMA] = mxCreateDoubleMatrix + (1, maxima_iterator-maxima_start, mxREAL) ; + maxima_end = maxima_iterator ; + maxima_iterator = maxima_start ; + M_pt = mxGetPr(out[MAXIMA]) ; + while(maxima_iterator != maxima_end) + { + *M_pt++ = *maxima_iterator++ ; + } + } + + /* Release space. */ + mxFree(offsets) ; + mxFree(neighbors) ; + mxFree(midx) ; + mxFree(maxima_start) ; + } + mxFree(dims) ; +} diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.m b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.m new file mode 100644 index 0000000..cee30ce --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.m @@ -0,0 +1,33 @@ +% SIFTLOCALMAX Find local maximizers +% SEL=SIFTLOCALMAX(F) returns the indexes of the local maximizers of +% the Q-dimensional array F. +% +% A local maximizer is an element whose value is greater than the +% value of all its neighbors. The neighbors of an element i1...iQ +% are the subscripts j1...jQ such that iq-1 <= jq <= iq (excluding +% i1...iQ itself). For example, if Q=1 the neighbors of an element +% are its predecessor and successor in the linear order; if Q=2, its +% neighbors are the elements immediately to its north, south, west, +% est, north-west, north-est, south-west and south-est +% (8-neighborhood). +% +% Points on the boundary of F are ignored (and never selected as +% local maximizers). +% +% SEL=SIFTLOCALMAX(F,THRESH) accepts an element as a mazimizer only +% if it is at least THRES greater than all its neighbors. +% +% SEL=SIFTLOCALMAX(F,THRESH,P) look for neighbors only in the first +% P dimensions of the Q-dimensional array F. This is useful to +% process F in ``slices''. +% +% REMARK. Matrices (2-array) with a singleton dimension are +% interpreted as vectors (1-array). So for example SIFTLOCALMAX([0 1 +% 0]) and SIFTLOCALMAX([0 1 0]') both return 2 as an aswer. However, +% if [0 1 0] is to be interpreted as a 1x2 matrix, then the correct +% answer is the empty set, as all elements are on the boundary. +% Unfortunately MATLAB does not distinguish between vectors and +% 2-matrices with a singleton dimension. To forece the +% interpretation of all matrices as 2-arrays, use +% SIFTLOCALMAX(F,TRESH,2) (but note that in this case the result is +% always empty!). diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.mexa64 b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.mexa64 new file mode 100755 index 0000000..2b8e264 Binary files /dev/null and b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.mexa64 differ diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.mexglx b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.mexglx new file mode 100755 index 0000000..613a74e Binary files /dev/null and b/SD-VBS/benchmarks/sift/src/matlab/siftlocalmax.mexglx differ diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftmatch.c b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.c new file mode 100644 index 0000000..1150176 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.c @@ -0,0 +1,250 @@ +/* file: siftmatch.c +** author: Andrea Vedaldi +** description: SIFT descriptor matching. +**/ + +/* AUTORIGHTS +Copyright (c) 2006 The Regents of the University of California. +All Rights Reserved. + +Created by Andrea Vedaldi +UCLA Vision Lab - Department of Computer Science + +Permission to use, copy, modify, and distribute this software and its +documentation for educational, research and non-profit purposes, +without fee, and without a written agreement is hereby granted, +provided that the above copyright notice, this paragraph and the +following three paragraphs appear in all copies. + +This software program and documentation are copyrighted by The Regents +of the University of California. The software program and +documentation are supplied "as is", without any accompanying services +from The Regents. The Regents does not warrant that the operation of +the program will be uninterrupted or error-free. The end-user +understands that the program was developed for research purposes and +is advised not to rely exclusively on the program for any reason. + +This software embodies a method for which the following patent has +been issued: "Method and apparatus for identifying scale invariant +features in an image and use of same for locating an object in an +image," David G. Lowe, US Patent 6,711,293 (March 23, +2004). Provisional application filed March 8, 1999. Asignee: The +University of British Columbia. + +IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +*/ + +#include"mexutils.c" + +#include +#include +#include + +#define greater(a,b) ((a) > (b)) +#define min(a,b) (((a)<(b))?(a):(b)) +#define max(a,b) (((a)>(b))?(a):(b)) + +#define TYPEOF_mxDOUBLE_CLASS double +#define TYPEOF_mxSINGLE_CLASS float +#define TYPEOF_mxINT8_CLASS char +#define TYPEOF_mxUINT8_CLASS unsigned char + +#define PROMOTE_mxDOUBLE_CLASS double +#define PROMOTE_mxSINGLE_CLASS float +#define PROMOTE_mxINT8_CLASS int +#define PROMOTE_mxUINT8_CLASS int + +#define MAXVAL_mxDOUBLE_CLASS mxGetInf() +#define MAXVAL_mxSINGLE_CLASS ((float)mxGetInf()) +#define MAXVAL_mxINT8_CLASS 0x7fffffff +#define MAXVAL_mxUINT8_CLASS 0x7fffffff + +typedef struct +{ + int k1 ; + int k2 ; + double score ; +} Pair ; + +/* + * This macro defines the matching function for abstract type; that + * is, it is a sort of C++ template. This is also a good illustration + * of why C++ is preferable for templates :-) + */ +#define _COMPARE_TEMPLATE(MXC) \ + Pair* compare_##MXC (Pair* pairs_iterator, \ + const TYPEOF_##MXC * L1_pt, \ + const TYPEOF_##MXC * L2_pt, \ + int K1, int K2, int ND, float thresh) \ + { \ + int k1, k2 ; \ + const PROMOTE_##MXC maxval = MAXVAL_##MXC ; \ + for(k1 = 0 ; k1 < K1 ; ++k1, L1_pt += ND ) { \ + \ + PROMOTE_##MXC best = maxval ; \ + PROMOTE_##MXC second_best = maxval ; \ + int bestk = -1 ; \ + \ + /* For each point P2[k2] in the second image... */ \ + for(k2 = 0 ; k2 < K2 ; ++k2, L2_pt += ND) { \ + \ + int bin ; \ + PROMOTE_##MXC acc = 0 ; \ + for(bin = 0 ; bin < ND ; ++bin) { \ + PROMOTE_##MXC delta = \ + ((PROMOTE_##MXC) L1_pt[bin]) - \ + ((PROMOTE_##MXC) L2_pt[bin]) ; \ + acc += delta*delta ; \ + } \ + \ + /* Filter the best and second best matching point. */ \ + if(acc < best) { \ + second_best = best ; \ + best = acc ; \ + bestk = k2 ; \ + } else if(acc < second_best) { \ + second_best = acc ; \ + } \ + } \ + \ + L2_pt -= ND*K2 ; \ + \ + /* Lowe's method: accept the match only if unique. */ \ + if(thresh * (float) best <= (float) second_best && \ + bestk != -1) { \ + pairs_iterator->k1 = k1 ; \ + pairs_iterator->k2 = bestk ; \ + pairs_iterator->score = best ; \ + pairs_iterator++ ; \ + } \ + } \ + \ + return pairs_iterator ; \ + } \ + +_COMPARE_TEMPLATE( mxDOUBLE_CLASS ) +_COMPARE_TEMPLATE( mxSINGLE_CLASS ) +_COMPARE_TEMPLATE( mxINT8_CLASS ) +_COMPARE_TEMPLATE( mxUINT8_CLASS ) + +void +mexFunction(int nout, mxArray *out[], + int nin, const mxArray *in[]) +{ + int K1, K2, ND ; + void* L1_pt ; + void* L2_pt ; + double thresh = 1.5 ; + mxClassID data_class ; + enum {L1=0,L2,THRESH} ; + enum {MATCHES=0,D} ; + + /* ------------------------------------------------------------------ + ** Check the arguments + ** --------------------------------------------------------------- */ + if (nin < 2) { + mexErrMsgTxt("At least two input arguments required"); + } else if (nout > 2) { + mexErrMsgTxt("Too many output arguments"); + } + + if(!mxIsNumeric(in[L1]) || + !mxIsNumeric(in[L2]) || + mxGetNumberOfDimensions(in[L1]) > 2 || + mxGetNumberOfDimensions(in[L2]) > 2) { + mexErrMsgTxt("L1 and L2 must be two dimensional numeric arrays") ; + } + + K1 = mxGetN(in[L1]) ; + K2 = mxGetN(in[L2]) ; + ND = mxGetM(in[L1]) ; + + if(mxGetM(in[L2]) != ND) { + mexErrMsgTxt("L1 and L2 must have the same number of rows") ; + } + + data_class = mxGetClassID(in[L1]) ; + if(mxGetClassID(in[L2]) != data_class) { + mexErrMsgTxt("L1 and L2 must be of the same class") ; + } + + L1_pt = mxGetData(in[L1]) ; + L2_pt = mxGetData(in[L2]) ; + + if(nin == 3) { + if(!uIsRealScalar(in[THRESH])) { + mexErrMsgTxt("THRESH should be a real scalar") ; + } + thresh = *mxGetPr(in[THRESH]) ; + } else if(nin > 3) { + mexErrMsgTxt("At most three arguments are allowed") ; + } + + /* ------------------------------------------------------------------ + ** Do the job + ** --------------------------------------------------------------- */ + { + Pair* pairs_begin = (Pair*) mxMalloc(sizeof(Pair) * (K1+K2)) ; + Pair* pairs_iterator = pairs_begin ; + + +#define _DISPATCH_COMPARE( MXC ) \ + case MXC : \ + pairs_iterator = compare_##MXC(pairs_iterator, \ + (const TYPEOF_##MXC*) L1_pt, \ + (const TYPEOF_##MXC*) L2_pt, \ + K1,K2,ND,thresh) ; \ + break ; \ + + switch (data_class) { + _DISPATCH_COMPARE( mxDOUBLE_CLASS ) ; + _DISPATCH_COMPARE( mxSINGLE_CLASS ) ; + _DISPATCH_COMPARE( mxINT8_CLASS ) ; + _DISPATCH_COMPARE( mxUINT8_CLASS ) ; + default : + mexErrMsgTxt("Unsupported numeric class") ; + break ; + } + + /* --------------------------------------------------------------- + * Finalize + * ------------------------------------------------------------ */ + { + Pair* pairs_end = pairs_iterator ; + double* M_pt ; + double* D_pt = NULL ; + + out[MATCHES] = mxCreateDoubleMatrix + (2, pairs_end-pairs_begin, mxREAL) ; + + M_pt = mxGetPr(out[MATCHES]) ; + + if(nout > 1) { + out[D] = mxCreateDoubleMatrix(1, + pairs_end-pairs_begin, + mxREAL) ; + D_pt = mxGetPr(out[D]) ; + } + + for(pairs_iterator = pairs_begin ; + pairs_iterator < pairs_end ; + ++pairs_iterator) { + *M_pt++ = pairs_iterator->k1 + 1 ; + *M_pt++ = pairs_iterator->k2 + 1 ; + if(nout > 1) { + *D_pt++ = pairs_iterator->score ; + } + } + } + mxFree(pairs_begin) ; + } +} diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftmatch.m b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.m new file mode 100644 index 0000000..d4398c1 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.m @@ -0,0 +1,60 @@ +% SIFTMATCH Match SIFT features +% MATCHES=SIFTMATCH(DESCR1, DESCR2) matches the two sets of SIFT +% descriptors DESCR1 and DESCR2. +% +% The function uses the same algorithm suggested by D. Lowe [1] to +% reject matches that are too ambiguous. +% +% SIFTMATCH(DESCR1, DESCR2, THRESH) uses [1] with the specified +% threshold THRESH. A descriptor D1 is matched to a descriptor D2 +% only if the distance d(D1,D2) multiplied by THRESH is not greather +% than the distance of D1 to all other descriptors. The default +% value of THRESH is 1.5. +% +% The storage class of the descriptors can be either DOUBLE, FLOAT, +% INT8 or UINT8. Usually interger classes are faster. +% +% [1] D. G. Lowe, +% `Distinctive image features from scale-invariant keypoints,' +% IJCV, vol. 2, no. 60, pp. 91–110, 2004. +% +% See also SIFT(), SIFTDESCRIPTOR(). + +% AUTORIGHTS +% Copyright (c) 2006 The Regents of the University of California. +% All Rights Reserved. +% +% Created by Andrea Vedaldi +% UCLA Vision Lab - Department of Computer Science +% +% Permission to use, copy, modify, and distribute this software and its +% documentation for educational, research and non-profit purposes, +% without fee, and without a written agreement is hereby granted, +% provided that the above copyright notice, this paragraph and the +% following three paragraphs appear in all copies. +% +% This software program and documentation are copyrighted by The Regents +% of the University of California. The software program and +% documentation are supplied "as is", without any accompanying services +% from The Regents. The Regents does not warrant that the operation of +% the program will be uninterrupted or error-free. The end-user +% understands that the program was developed for research purposes and +% is advised not to rely exclusively on the program for any reason. +% +% This software embodies a method for which the following patent has +% been issued: "Method and apparatus for identifying scale invariant +% features in an image and use of same for locating an object in an +% image," David G. Lowe, US Patent 6,711,293 (March 23, +% 2004). Provisional application filed March 8, 1999. Asignee: The +% University of British Columbia. +% +% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexa64 b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexa64 new file mode 100755 index 0000000..e58eaf5 Binary files /dev/null and b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexa64 differ diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexglx b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexglx new file mode 100755 index 0000000..d5e32f7 Binary files /dev/null and b/SD-VBS/benchmarks/sift/src/matlab/siftmatch.mexglx differ diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftormx.c b/SD-VBS/benchmarks/sift/src/matlab/siftormx.c new file mode 100644 index 0000000..40380a2 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/siftormx.c @@ -0,0 +1,251 @@ +/* file: siftormx.c +** author: Andrea Vedaldi +** description: Computes peaks of orientation histogram. +**/ + +/* AUTORIGHTS +Copyright (c) 2006 The Regents of the University of California. +All Rights Reserved. + +Created by Andrea Vedaldi +UCLA Vision Lab - Department of Computer Science + +Permission to use, copy, modify, and distribute this software and its +documentation for educational, research and non-profit purposes, +without fee, and without a written agreement is hereby granted, +provided that the above copyright notice, this paragraph and the +following three paragraphs appear in all copies. + +This software program and documentation are copyrighted by The Regents +of the University of California. The software program and +documentation are supplied "as is", without any accompanying services +from The Regents. The Regents does not warrant that the operation of +the program will be uninterrupted or error-free. The end-user +understands that the program was developed for research purposes and +is advised not to rely exclusively on the program for any reason. + +This software embodies a method for which the following patent has +been issued: "Method and apparatus for identifying scale invariant +features in an image and use of same for locating an object in an +image," David G. Lowe, US Patent 6,711,293 (March 23, +2004). Provisional application filed March 8, 1999. Asignee: The +University of British Columbia. + +IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +*/ + +#include"mex.h" +#include +#include +#include + +#include + +#define greater(a,b) a > b +#define min(a,b) (((a)<(b))?(a):(b)) +#define max(a,b) (((a)>(b))?(a):(b)) + +const double win_factor = 1.5 ; +#define NBINS 36 + +void +mexFunction(int nout, mxArray *out[], + int nin, const mxArray *in[]) +{ + int M,N,S,smin,K ; + const int* dimensions ; + const double* P_pt ; + const double* G_pt ; + double* TH_pt ; + double sigma0 ; + double H_pt [ NBINS ] ; + + enum {IN_P=0,IN_G,IN_S,IN_SMIN,IN_SIGMA0} ; + enum {OUT_Q=0} ; + + /* ----------------------------------------------------------------- + ** Check the arguments + ** -------------------------------------------------------------- */ + if (nin != 5) { + mexErrMsgTxt("Exactly five input arguments required."); + } else if (nout > 1) { + mexErrMsgTxt("Too many output arguments."); + } + + if( !uIsRealScalar(in[IN_S]) ) { + mexErrMsgTxt("S should be a real scalar") ; + } + + if( !uIsRealScalar(in[IN_SMIN]) ) { + mexErrMsgTxt("SMIN should be a real scalar") ; + } + + if( !uIsRealScalar(in[IN_SIGMA0]) ) { + mexErrMsgTxt("SIGMA0 should be a real scalar") ; + } + + if( !uIsRealMatrix(in[IN_P],3,-1)) { + mexErrMsgTxt("P should be a 3xK real matrix") ; + } + + if(mxGetNumberOfDimensions(in[IN_G]) != 3) { + mexErrMsgTxt("SSO must be a three dimensional array") ; + } + + dimensions = mxGetDimensions(in[IN_G]) ; + M = dimensions[0] ; + N = dimensions[1] ; + S = (int)(*mxGetPr(in[IN_S])) ; + smin = (int)(*mxGetPr(in[IN_SMIN])) ; + sigma0 = *mxGetPr(in[IN_SIGMA0]) ; + + K = mxGetN(in[IN_P]) ; + P_pt = mxGetPr(in[IN_P]) ; + G_pt = mxGetPr(in[IN_G]) ; + + + /* If the input array is empty, then output an empty array as well. */ + if(K == 0) { + out[OUT_Q] = mxCreateDoubleMatrix(4,0,mxREAL) ; + return ; + } + + /* ------------------------------------------------------------------ + * Do the job + * --------------------------------------------------------------- */ + { + int p ; + const int yo = 1 ; + const int xo = M ; + const int so = M*N ; + + int buffer_size = K*4 ; + double* buffer_start = (double*) mxMalloc( buffer_size *sizeof(double)) ; + double* buffer_iterator = buffer_start ; + double* buffer_end = buffer_start + buffer_size ; + + for(p = 0 ; p < K ; ++p, TH_pt += 2) { + const double x = *P_pt++ ; + const double y = *P_pt++ ; + const double s = *P_pt++ ; + int xi = ((int) (x+0.5)) ; /* Round them off. */ + int yi = ((int) (y+0.5)) ; + int si = ((int) (s+0.5)) - smin ; + int xs ; + int ys ; + double sigmaw = win_factor * sigma0 * pow(2, ((double)s) / S) ; + int W = (int) ceil(3.0 * sigmaw) ; + int bin ; + const double* pt ; + + /* Make sure that the rounded off keypoint index is within bound. + */ + if(xi < 0 || + xi > N-1 || + yi < 0 || + yi > M-1 || + si < 0 || + si > dimensions[2]-1 ) { + mexPrintf("Dropping %d: W %d x %d y %d si [%d,%d,%d,%d]\n",p,W,xi,yi,si,M,N,dimensions[2]) ; + continue ; + } + + /* Clear histogram buffer. */ + { + int i ; + for(i = 0 ; i < NBINS ; ++i) + H_pt[i] = 0 ; + } + + pt = G_pt + xi*xo + yi*yo + si*so ; + +#define at(dx,dy) (*(pt + (dx)*xo + (dy)*yo)) + + for(xs = max(-W, 1-xi) ; xs <= min(+W, N -2 -xi) ; ++xs) { + for(ys = max(-W, 1-yi) ; ys <= min(+W, M -2 -yi) ; ++ys) { + double Dx = 0.5 * ( at(xs+1,ys) - at(xs-1,ys) ) ; + double Dy = 0.5 * ( at(xs,ys+1) - at(xs,ys-1) ) ; + double dx = ((double)(xi+xs)) - x; + double dy = ((double)(yi+ys)) - y; + + if(dx*dx + dy*dy > W*W+0.5) continue ; + + { + double win = exp( - (dx*dx + dy*dy)/(2*sigmaw*sigmaw) ) ; + double mod = sqrt(Dx*Dx + Dy*Dy) ; + double theta = fmod(atan2(Dy, Dx) + 2*M_PI, 2*M_PI) ; + bin = (int) floor( NBINS * theta / (2*M_PI) ) ; + H_pt[bin] += mod*win ; + } + } + } + + /* Smooth histogram */ + { + int iter, i ; + for (iter = 0; iter < 6; iter++) { + double prev; + prev = H_pt[NBINS-1]; + for (i = 0; i < NBINS; i++) { + float newh = (prev + H_pt[i] + H_pt[(i+1) % NBINS]) / 3.0; + prev = H_pt[i] ; + H_pt[i] = newh ; + } + } + } + + /* Find most strong peaks. */ + { + int i ; + double maxh = H_pt[0] ; + for(i = 1 ; i < NBINS ; ++i) + maxh = max(maxh, H_pt[i]) ; + + for(i = 0 ; i < NBINS ; ++i) { + double h0 = H_pt[i] ; + double hm = H_pt[(i-1+NBINS) % NBINS] ; + double hp = H_pt[(i+1+NBINS) % NBINS] ; + + if( h0 > 0.8*maxh && h0 > hm && h0 > hp ) { + + double di = -0.5 * (hp-hm) / (hp+hm-2*h0) ; /*di=0;*/ + double th = 2*M_PI*(i+di+0.5)/NBINS ; + + if( buffer_iterator == buffer_end ) { + int offset = buffer_iterator - buffer_start ; + buffer_size += 4*max(1, K/16) ; + buffer_start = (double*) mxRealloc(buffer_start, + buffer_size*sizeof(double)) ; + buffer_end = buffer_start + buffer_size ; + buffer_iterator = buffer_start + offset ; + } + + *buffer_iterator++ = x ; + *buffer_iterator++ = y ; + *buffer_iterator++ = s ; + *buffer_iterator++ = th ; + } + } /* Scan histogram */ + } /* Find peaks */ + } + + /* Save back the result. */ + { + double* result ; + int NL = (buffer_iterator - buffer_start)/4 ; + out[OUT_Q] = mxCreateDoubleMatrix(4, NL, mxREAL) ; + result = mxGetPr(out[OUT_Q]); + memcpy(result, buffer_start, sizeof(double) * 4 * NL) ; + } + mxFree(buffer_start) ; + } +} diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftormx.mexa64 b/SD-VBS/benchmarks/sift/src/matlab/siftormx.mexa64 new file mode 100755 index 0000000..ae1f787 Binary files /dev/null and b/SD-VBS/benchmarks/sift/src/matlab/siftormx.mexa64 differ diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftormx.mexglx b/SD-VBS/benchmarks/sift/src/matlab/siftormx.mexglx new file mode 100755 index 0000000..de90e4b Binary files /dev/null and b/SD-VBS/benchmarks/sift/src/matlab/siftormx.mexglx differ diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftread.m b/SD-VBS/benchmarks/sift/src/matlab/siftread.m new file mode 100644 index 0000000..76945c8 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/siftread.m @@ -0,0 +1,101 @@ +function [frames,descriptors] = siftread(file) +% SIFTREAD Read Lowe's SIFT implementation data files +% [FRAMES, DESCRIPTORS] = READSIFT(FILE) reads the frames and the +% descriptors from the specified file. The function reads files +% produced by Lowe's SIFT implementation. +% +% FRAMES and DESCRIPTORS have the same format used by SIFT(). +% +% REMARK. Lowe's and our implementations use a silightly different +% convention to store the orientation of the frame. When the file +% is read, the orientation is changed to match our convention. +% +% See also SIFT(). + +% AUTORIGHTS +% Copyright (c) 2006 The Regents of the University of California. +% All Rights Reserved. +% +% Created by Andrea Vedaldi +% UCLA Vision Lab - Department of Computer Science +% +% Permission to use, copy, modify, and distribute this software and its +% documentation for educational, research and non-profit purposes, +% without fee, and without a written agreement is hereby granted, +% provided that the above copyright notice, this paragraph and the +% following three paragraphs appear in all copies. +% +% This software program and documentation are copyrighted by The Regents +% of the University of California. The software program and +% documentation are supplied "as is", without any accompanying services +% from The Regents. The Regents does not warrant that the operation of +% the program will be uninterrupted or error-free. The end-user +% understands that the program was developed for research purposes and +% is advised not to rely exclusively on the program for any reason. +% +% This software embodies a method for which the following patent has +% been issued: "Method and apparatus for identifying scale invariant +% features in an image and use of same for locating an object in an +% image," David G. Lowe, US Patent 6,711,293 (March 23, +% 2004). Provisional application filed March 8, 1999. Asignee: The +% University of British Columbia. +% +% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. + +verbosity=0 ; + +g = fopen(file, 'r'); +if g == -1 + error(['Could not open file ''', file, '''.']) ; +end +[header, count] = fscanf(g, '%d', [1 2]) ; +if count ~= 2 + error('Invalid keypoint file header.'); +end +K = header(1) ; +DL = header(2) ; + +if(verbosity > 0) + fprintf('%d keypoints, %d descriptor length.\n', K, DL) ; +end + +%creates two output matrices +P = zeros(4,K) ; +L = zeros(DL,K) ; + +%parse tmp.key +for k = 1:K + + % Record format: i,j,s,th + [record, count] = fscanf(g, '%f', [1 4]) ; + if count ~= 4 + error(... + sprintf('Invalid keypoint file (parsing keypoint %d)',k) ); + end + P(:,k) = record(:) ; + + % Record format: descriptor + [record, count] = fscanf(g, '%d', [1 DL]) ; + if count ~= DL + error(... + sprintf('Invalid keypoint file (parsing keypoint %d)',k) ); + end + L(:,k) = record(:) ; + +end +fclose(g) ; + +L=double(L) ; +P(1:2,:)=flipud(P(1:2,:)) ; % i,j -> x,y + +frames=[ P(1:2,:) ; P(3,:) ; -P(4,:) ] ; +descriptors = L ; diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.c b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.c new file mode 100644 index 0000000..a1ef8c7 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.c @@ -0,0 +1,342 @@ +/* file: siftrefinemx.c + +Benchmark - sift +Data set - sqcif + + < M A T L A B > + Copyright 1984-2006 The MathWorks, Inc. + Version 7.3.0.298 (R2006b) + August 03, 2006 + + + To get started, type one of these: helpwin, helpdesk, or demo. + For product information, visit www.mathworks.com. + +Warning: Function /h/g2/kvs/checkParallel/sdvbs-svn/common/matlab/randn.m has the same name as a MATLAB builtin. We suggest you rename the function to avoid a potential name conflict. +> In path at 113 + In script_run_profile at 3 +Warning: You are using gcc version "4.1.1". The earliest gcc version supported +with mex is "3.4.0". The latest version tested for use with mex is "3.4.5". +To download a different version of gcc, visit http://gcc.gnu.org +Warning: You are using gcc version "4.1.1". The earliest gcc version supported +with mex is "3.4.0". The latest version tested for use with mex is "3.4.5". +To download a different version of gcc, visit http://gcc.gnu.org +Warning: You are using gcc version "4.1.1". The earliest gcc version supported +with mex is "3.4.0". The latest version tested for use with mex is "3.4.5". +To download a different version of gcc, visit http://gcc.gnu.org +Warning: You are using gcc version "4.1.1". The earliest gcc version supported +with mex is "3.4.0". The latest version tested for use with mex is "3.4.5". +To download a different version of gcc, visit http://gcc.gnu.org +Warning: You are using gcc version "4.1.1". The earliest gcc version supported +with mex is "3.4.0". The latest version tested for use with mex is "3.4.5". +To download a different version of gcc, visit http://gcc.gnu.org +Warning: You are using gcc version "4.1.1". The earliest gcc version supported +with mex is "3.4.0". The latest version tested for use with mex is "3.4.5". +To download a different version of gcc, visit http://gcc.gnu.org +Input size - (96x96) +** author: Andrea Vedaldi +** description: Subpixel localization, thresholding and off-edge test. +**/ + +/* AUTORIGHTS +Copyright (c) 2006 The Regents of the University of California. +All Rights Reserved. + +Created by Andrea Vedaldi +UCLA Vision Lab - Department of Computer Science + +Permission to use, copy, modify, and distribute this software and its +documentation for educational, research and non-profit purposes, +without fee, and without a written agreement is hereby granted, +provided that the above copyright notice, this paragraph and the +following three paragraphs appear in all copies. + +This software program and documentation are copyrighted by The Regents +of the University of California. The software program and +documentation are supplied "as is", without any accompanying services +from The Regents. The Regents does not warrant that the operation of +the program will be uninterrupted or error-free. The end-user +understands that the program was developed for research purposes and +is advised not to rely exclusively on the program for any reason. + +This software embodies a method for which the following patent has +been issued: "Method and apparatus for identifying scale invariant +features in an image and use of same for locating an object in an +image," David G. Lowe, US Patent 6,711,293 (March 23, +2004). Provisional application filed March 8, 1999. Asignee: The +University of British Columbia. + +IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +*/ + +#include"mex.h" + +#include + +#include +#include + +/* Prototype of DGESV LAPACK function for the solution of a linear system. */ +#ifdef WINDOWS +#define DGESV dgesv +#undef min +#undef max +#else +#define DGESV dgesv_ +#endif + +#ifdef WINDOWS +#ifdef __cplusplus__ +extern "C" { + extern int DGESV(int *n, int *nrhs, double *a, int *lda, + int *ipiv, double *b, int *ldb, int *info) ; +} +#else + extern int DGESV(int *n, int *nrhs, double *a, int *lda, + int *ipiv, double *b, int *ldb, int *info) ; +#define sqrtf(x) ((float)sqrt((double)x) +#define powf(x) ((float)pow((double)x) +#define fabsf(x) ((float)fabs((double)x) +#endif +#else +extern int DGESV(int *n, int *nrhs, double *a, int *lda, + int *ipiv, double *b, int *ldb, int *info) ; +#endif + +#define greater(a,b) ((a) > (b)) +#define min(a,b) (((a)<(b))?(a):(b)) +#define max(a,b) (((a)>(b))?(a):(b)) + +const int max_iter = 5 ; + +void +mexFunction(int nout, mxArray *out[], + int nin, const mxArray *in[]) +{ + int M,N,S,smin,K ; + const int* dimensions ; + const double* P_pt ; + const double* D_pt ; + double threshold = 0.01 ; /*0.02 ;*/ + double r = 10.0 ; + double* result ; + enum {IN_P=0,IN_D,IN_SMIN,IN_THRESHOLD,IN_R} ; + enum {OUT_Q=0} ; + + /* ----------------------------------------------------------------- + ** Check the arguments + ** -------------------------------------------------------------- */ + if (nin < 3) { + mexErrMsgTxt("At least three input arguments required."); + } else if (nout > 1) { + mexErrMsgTxt("Too many output arguments."); + } + + if( !uIsRealMatrix(in[IN_P],3,-1) ) { + mexErrMsgTxt("P must be a 3xK real matrix") ; + } + + if( !mxIsDouble(in[IN_D]) || mxGetNumberOfDimensions(in[IN_D]) != 3) { + mexErrMsgTxt("G must be a three dimensional real array.") ; + } + + if( !uIsRealScalar(in[IN_SMIN]) ) { + mexErrMsgTxt("SMIN must be a real scalar.") ; + } + + if(nin >= 4) { + if(!uIsRealScalar(in[IN_THRESHOLD])) { + mexErrMsgTxt("THRESHOLD must be a real scalar.") ; + } + threshold = *mxGetPr(in[IN_THRESHOLD]) ; + } + + if(nin >= 5) { + if(!uIsRealScalar(in[IN_R])) { + mexErrMsgTxt("R must be a real scalar.") ; + } + r = *mxGetPr(in[IN_R]) ; + } + + dimensions = mxGetDimensions(in[IN_D]) ; + M = dimensions[0] ; + N = dimensions[1] ; + S = dimensions[2] ; + smin = (int)(*mxGetPr(in[IN_SMIN])) ; + + if(S < 3 || M < 3 || N < 3) { + mexErrMsgTxt("All dimensions of DOG must be not less than 3.") ; + } + + K = mxGetN(in[IN_P]) ; + P_pt = mxGetPr(in[IN_P]) ; + D_pt = mxGetPr(in[IN_D]) ; + + /* If the input array is empty, then output an empty array as well. */ + if( K == 0) { + out[OUT_Q] = mxDuplicateArray(in[IN_P]) ; + return ; + } + + /* ----------------------------------------------------------------- + * Do the job + * -------------------------------------------------------------- */ + { + double* buffer = (double*) mxMalloc(K*3*sizeof(double)) ; + double* buffer_iterator = buffer ; + int p ; + const int yo = 1 ; + const int xo = M ; + const int so = M*N ; + +/* printf("Actual values = %d\n\n", K); +*/ + for(p = 0 ; p < K ; ++p) { + int x = ((int)*P_pt++) ; + int y = ((int)*P_pt++) ; +/* printf("%d\t%d\n", ((int)*P_pt), smin); +*/ + int s = ((int)*P_pt++) - smin ; + int iter ; + double b[3] ; + + /* Local maxima extracted from the DOG + * have coorrinates 1<=x<=N-2, 1<=y<=M-2 + * and 1<=s-mins<=S-2. This is also the range of the points + * that we can refine. + */ +/* printf("%d\t%d\t%d\t%d\t%d\t%d\n", x, N-2,y,M-2, s, S-2); +*/ + if(x < 1 || x > N-2 || + y < 1 || y > M-2 || + s < 1 || s > S-2) { + continue ; + } + +#define at(dx,dy,ds) (*(pt + (dx)*xo + (dy)*yo + (ds)*so)) + + { + const double* pt = D_pt + y*yo + x*xo + s*so ; + + double Dx=0,Dy=0,Ds=0,Dxx=0,Dyy=0,Dss=0,Dxy=0,Dxs=0,Dys=0 ; + int dx = 0 ; + int dy = 0 ; +/* printf("%d\t%d\t%d\t%d\t%d\t%d\t%d\t%f\t%f\n",S, y, yo, x, xo, s, so, *D_pt, *pt); +*/ + for(iter = 0 ; iter < max_iter ; ++iter) { + + double A[3*3] ; + int ipiv[3] ; + int n = 3 ; + int one = 1 ; + int info = 0 ; + +#define Aat(i,j) (A[(i)+(j)*3]) + + x += dx ; + y += dy ; + pt = D_pt + y*yo + x*xo + s*so ; + + /* Compute the gradient. */ + Dx = 0.5 * (at(+1,0,0) - at(-1,0,0)) ; + Dy = 0.5 * (at(0,+1,0) - at(0,-1,0)); + Ds = 0.5 * (at(0,0,+1) - at(0,0,-1)) ; + + /* Compute the Hessian. */ + Dxx = (at(+1,0,0) + at(-1,0,0) - 2.0 * at(0,0,0)) ; + Dyy = (at(0,+1,0) + at(0,-1,0) - 2.0 * at(0,0,0)) ; + Dss = (at(0,0,+1) + at(0,0,-1) - 2.0 * at(0,0,0)) ; + + Dxy = 0.25 * ( at(+1,+1,0) + at(-1,-1,0) - at(-1,+1,0) - at(+1,-1,0) ) ; + Dxs = 0.25 * ( at(+1,0,+1) + at(-1,0,-1) - at(-1,0,+1) - at(+1,0,-1) ) ; + Dys = 0.25 * ( at(0,+1,+1) + at(0,-1,-1) - at(0,-1,+1) - at(0,+1,-1) ) ; + + /* Solve linear system. */ + Aat(0,0) = Dxx ; + Aat(1,1) = Dyy ; + Aat(2,2) = Dss ; + Aat(0,1) = Aat(1,0) = Dxy ; + Aat(0,2) = Aat(2,0) = Dxs ; + Aat(1,2) = Aat(2,1) = Dys ; + + b[0] = - Dx ; + b[1] = - Dy ; + b[2] = - Ds ; + + /* DGESV (&n, &one, A, &n, ipiv, b, &n, &info) ; + */ + /* If the translation of the keypoint is big, move the keypoint + * and re-iterate the computation. Otherwise we are all set. + */ + dx= ((b[0] > 0.6 && x < N-2) ? 1 : 0 ) + + ((b[0] < -0.6 && x > 1 ) ? -1 : 0 ) ; + + dy= ((b[1] > 0.6 && y < M-2) ? 1 : 0 ) + + ((b[1] < -0.6 && y > 1 ) ? -1 : 0 ) ; + + if( dx == 0 && dy == 0 ) break ; + + } + + { + double val = at(0,0,0) + 0.5 * (Dx * b[0] + Dy * b[1] + Ds * b[2]) ; + double score = (Dxx+Dyy)*(Dxx+Dyy) / (Dxx*Dyy - Dxy*Dxy) ; + double xn = x + b[0] ; + double yn = y + b[1] ; + double sn = s + b[2] ; + +/* printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", fabs(val),threshold,score,(r+1)*(r+1)/r,fabs(b[0]), fabs(b[1]), fabs(b[2]),xn,yn,sn,r); +*/ + if(fabs(val) > threshold && + score < (r+1)*(r+1)/r && + score >= 0 && + fabs(b[0]) < 1.5 && + fabs(b[1]) < 1.5 && + fabs(b[2]) < 1.5 && + xn >= 0 && + xn <= N-1 && + yn >= 0 && + yn <= M-1 && + sn >= 0 && + sn <= S-1) { + *buffer_iterator++ = xn ; + *buffer_iterator++ = yn ; + *buffer_iterator++ = sn+smin ; + } + } + } + } + + /* Copy the result into an array. */ + { + int i; + int NL = (buffer_iterator - buffer)/3 ; +/* printf("%NL VALUE = %d\t%d\t%d\n", NL, buffer_iterator, buffer); +*/ + out[OUT_Q] = mxCreateDoubleMatrix(3, NL, mxREAL) ; + result = mxGetPr(out[OUT_Q]); + for(i=0; i<(3*NL); i++) + { + result[i] = buffer[i]; +/* printf("%f\t", buffer[i]); +*/ + } +/* printf("\n"); + memcpy(result, buffer, sizeof(double) * 3 * NL) ; +*/ + } + mxFree(buffer) ; + } + +} + + diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.m b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.m new file mode 100644 index 0000000..222d610 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.m @@ -0,0 +1,61 @@ +% SIFTREFINEMX Subpixel localization, thresholding and on-edge test +% Q = SIFTREFINEMX(P, OCTAVE, SMIN) refines, thresholds and performs +% the on-edge test for the SIFT frames P extracted from the DOG +% octave OCTAVE with parameter SMIN (see GAUSSIANSS()). +% +% Q = SIFTREFINEMX(P, OCTAVE, SMIN, THRESH, R) specifies custom +% values for the local maximum strength threshold THRESH and the +% local maximum peakedeness threshold R. +% +% OCTAVE is an octave of the Difference Of Gaussian scale space. P +% is a 3xK matrix specifying the indexes (X,Y,S) of the points of +% extremum of the octave OCTAVE. The spatial indexes X,Y are integer +% with base zero. The scale index S is integer with base SMIN and +% represents a scale sublevel in the specified octave. +% +% The function returns a matrix Q containing the refined keypoints. +% The matrix has the same format as P, except that the indexes are +% now fractional. The function drops the points that do not satisfy +% the strength and peakedness tests. +% +% See also SIFT(). + +% AUTORIGHTS +% Copyright (c) 2006 The Regents of the University of California. +% All Rights Reserved. +% +% Created by Andrea Vedaldi +% UCLA Vision Lab - Department of Computer Science +% +% Permission to use, copy, modify, and distribute this software and its +% documentation for educational, research and non-profit purposes, +% without fee, and without a written agreement is hereby granted, +% provided that the above copyright notice, this paragraph and the +% following three paragraphs appear in all copies. +% +% This software program and documentation are copyrighted by The Regents +% of the University of California. The software program and +% documentation are supplied "as is", without any accompanying services +% from The Regents. The Regents does not warrant that the operation of +% the program will be uninterrupted or error-free. The end-user +% understands that the program was developed for research purposes and +% is advised not to rely exclusively on the program for any reason. +% +% This software embodies a method for which the following patent has +% been issued: "Method and apparatus for identifying scale invariant +% features in an image and use of same for locating an object in an +% image," David G. Lowe, US Patent 6,711,293 (March 23, +% 2004). Provisional application filed March 8, 1999. Asignee: The +% University of British Columbia. +% +% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. + diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.mexa64 b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.mexa64 new file mode 100755 index 0000000..c91556b Binary files /dev/null and b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.mexa64 differ diff --git a/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.mexglx b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.mexglx new file mode 100755 index 0000000..74b312c Binary files /dev/null and b/SD-VBS/benchmarks/sift/src/matlab/siftrefinemx.mexglx differ diff --git a/SD-VBS/benchmarks/sift/src/matlab/tightsubplot.m b/SD-VBS/benchmarks/sift/src/matlab/tightsubplot.m new file mode 100644 index 0000000..b790ef0 --- /dev/null +++ b/SD-VBS/benchmarks/sift/src/matlab/tightsubplot.m @@ -0,0 +1,151 @@ +function H = tightsubplot(varargin) +% TIGHTSUBPLOT Tiles axes without wasting space +% H = TIGHTSUBPLOT(K,P) returns an handle to the P-th axis in a +% regular grid of K axes. The K axes are numbered from left to right +% and from top to bottom. The function operates similarly to +% SUBPLOT(), but by default it does not put any margin between axes. +% +% H = TIGHTSUBPLOT(M,N,P) retursn an handle to the P-th axes in a +% regular subdivision with M rows and N columns. +% +% The function accepts the following option-value pairs: +% +% 'Spacing' [0] +% Set extra spacing between axes. The space is added between the +% inner or outer boxes, depending on the setting below. +% +% 'Box' ['inner'] (** ONLY >R14 **) +% If set to 'outer', the function displaces the axes by their +% outer box, thus protecting title and labels. Unfortunately +% MATLAB typically picks unnecessarily large insets, so that a bit +% of space is wasted in this case. If set to 'inner', the +% function uses the inner box. This causes the instets of nearby +% axes to overlap, but it is very space conservative. +% +% REMARK. While SUBPLOT kills any pre-existing axes that overalps a +% new one, this function does not. +% +% See also SUBPLOT(). + +% AUTORIGHTS +% Copyright (c) 2006 The Regents of the University of California. +% All Rights Reserved. +% +% Created by Andrea Vedaldi +% UCLA Vision Lab - Department of Computer Science +% +% Permission to use, copy, modify, and distribute this software and its +% documentation for educational, research and non-profit purposes, +% without fee, and without a written agreement is hereby granted, +% provided that the above copyright notice, this paragraph and the +% following three paragraphs appear in all copies. +% +% This software program and documentation are copyrighted by The Regents +% of the University of California. The software program and +% documentation are supplied "as is", without any accompanying services +% from The Regents. The Regents does not warrant that the operation of +% the program will be uninterrupted or error-free. The end-user +% understands that the program was developed for research purposes and +% is advised not to rely exclusively on the program for any reason. +% +% This software embodies a method for which the following patent has +% been issued: "Method and apparatus for identifying scale invariant +% features in an image and use of same for locating an object in an +% image," David G. Lowe, US Patent 6,711,293 (March 23, +% 2004). Provisional application filed March 8, 1999. Asignee: The +% University of British Columbia. +% +% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY +% FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, +% INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND +% ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN +% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF +% CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT +% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +% A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" +% BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE +% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. + +sp=0.0 ; +use_outer=0 ; + +% -------------------------------------------------------------------- +% Parse arguments +% -------------------------------------------------------------------- +K=varargin{1} ; +p=varargin{2} ; +N = ceil(sqrt(K)) ; +M = ceil(K/N) ; + +a=3 ; +NA = length(varargin) ; +if NA > 2 + if isa(varargin{3},'char') + % Called with K and p + else + % Called with M,N and p + a = 4 ; + M = K ; + N = p ; + p = varargin{3} ; + end +end + +for a=a:2:NA + switch varargin{a} + case 'Spacing' + sp=varargin{a+1} ; + case 'Box' + switch varargin{a+1} + case 'inner' + use_outer = 0 ; + case 'outer' + if ~strcmp(version('-release'), '14') + %warning(['Box option supported only on MATALB 14']) ; + continue; + end + use_outer = 1 ; + otherwise + error(['Box is either ''inner'' or ''outer''']) ; + end + otherwise + error(['Uknown parameter ''', varargin{a}, '''.']) ; + end +end + +% -------------------------------------------------------------------- +% Check the arguments +% -------------------------------------------------------------------- + +[j,i]=ind2sub([N M],p) ; +i=i-1 ; +j=j-1 ; + +dt = sp/2 ; +db = sp/2 ; +dl = sp/2 ; +dr = sp/2 ; + +pos = [ j*1/N+dl,... + 1-i*1/M-1/M-db,... + 1/N-dl-dr,... + 1/M-dt-db] ; + +switch use_outer + case 0 + H = findobj(gcf, 'Type', 'axes', 'Position', pos) ; + if(isempty(H)) + H = axes('Position', pos) ; + else + axes(H) ; + end + + case 1 + H = findobj(gcf, 'Type', 'axes', 'OuterPosition', pos) ; + if(isempty(H)) + H = axes('ActivePositionProperty', 'outerposition',... + 'OuterPosition', pos) ; + else + axes(H) ; + end +end -- cgit v1.2.2