From 14c854543b1a3cf344a371a5b45595657f95786b Mon Sep 17 00:00:00 2001 From: leochanj105 Date: Tue, 26 Oct 2021 13:22:32 -0400 Subject: add deadlines and costs --- 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 | 80 ------ 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, 1357 deletions(-) delete mode 100644 SD-VBS/benchmarks/sift/src/c/diffss.c delete mode 100644 SD-VBS/benchmarks/sift/src/c/doubleSize.c delete mode 100644 SD-VBS/benchmarks/sift/src/c/filterBoundaryPoints.c delete mode 100644 SD-VBS/benchmarks/sift/src/c/gaussianss.c delete mode 100644 SD-VBS/benchmarks/sift/src/c/halveSize.c delete mode 100644 SD-VBS/benchmarks/sift/src/c/imsmooth.c delete mode 100644 SD-VBS/benchmarks/sift/src/c/script_sift.c delete mode 100644 SD-VBS/benchmarks/sift/src/c/sift.c delete mode 100644 SD-VBS/benchmarks/sift/src/c/sift.h delete mode 100644 SD-VBS/benchmarks/sift/src/c/siftlocalmax.c delete mode 100644 SD-VBS/benchmarks/sift/src/c/siftrefinemx.c (limited to 'SD-VBS/benchmarks/sift/src/c') diff --git a/SD-VBS/benchmarks/sift/src/c/diffss.c b/SD-VBS/benchmarks/sift/src/c/diffss.c deleted file mode 100644 index 683cbe7..0000000 --- a/SD-VBS/benchmarks/sift/src/c/diffss.c +++ /dev/null @@ -1,53 +0,0 @@ -/******************************** -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 deleted file mode 100644 index 52dd95b..0000000 --- a/SD-VBS/benchmarks/sift/src/c/gaussianss.c +++ /dev/null @@ -1,185 +0,0 @@ -/******************************** -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 deleted file mode 100644 index fe1e536..0000000 --- a/SD-VBS/benchmarks/sift/src/c/halveSize.c +++ /dev/null @@ -1,35 +0,0 @@ -/******************************** -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<<24 -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]; - scanf("%s", imSrc); - im = readImage(imSrc); - image = fiDeepCopy(im); - rows = image->height; - cols = image->width; - - - for_each_job{ - image = fiDeepCopy(im); - normalizeImage(image); - /** Extract sift features for the normalized image **/ - frames = sift(image); - } - -#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 deleted file mode 100644 index ea0d308..0000000 --- a/SD-VBS/benchmarks/sift/src/c/sift.c +++ /dev/null @@ -1,314 +0,0 @@ -/******************************** -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 deleted file mode 100644 index d75bcbd..0000000 --- a/SD-VBS/benchmarks/sift/src/c/siftlocalmax.c +++ /dev/null @@ -1,247 +0,0 @@ -#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