From d17b33131c14864bd1eae275f49a3f148e21cf29 Mon Sep 17 00:00:00 2001 From: Leo Chan Date: Thu, 22 Oct 2020 01:53:21 -0400 Subject: Squashed commit of the sb-vbs branch. Includes the SD-VBS benchmarks modified to: - Use libextra to loop as realtime jobs - Preallocate memory before starting their main computation - Accept input via stdin instead of via argc Does not include the SD-VBS matlab code. Fixes libextra execution in LITMUS^RT. --- 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 +++++++++++++ 11 files changed, 1360 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 (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