// Copyright (c) 2008-2011, Guoshen Yu // Copyright (c) 2008-2011, Jean-Michel Morel // // WARNING: // This file implements an algorithm possibly linked to the patent // // Jean-Michel Morel and Guoshen Yu, Method and device for the invariant // affine recognition recognition of shapes (WO/2009/150361), patent pending. // // This file is made available for the exclusive aim of serving as // scientific tool to verify of the soundness and // completeness of the algorithm description. Compilation, // execution and redistribution of this file may violate exclusive // patents rights in certain countries. // The situation being different for every country and changing // over time, it is your responsibility to determine which patent // rights restrictions apply to you before you compile, use, // modify, or redistribute this file. A patent lawyer is qualified // to make this determination. // If and only if they don't conflict with any patent terms, you // can benefit from the following license terms attached to this // file. // // This program is provided for scientific and educational only: // you can use and/or modify it for these purposes, but you are // not allowed to redistribute this work or derivative works in // source or executable form. A license must be obtained from the // patent right holders for any other use. // // //*------------------------ compute_asift_keypoints -------------------------*/ // Compute the ASIFT keypoints on the input image. // // Please report bugs and/or send comments to Guoshen Yu yu@cmap.polytechnique.fr // // Reference: J.M. Morel and G.Yu, ASIFT: A New Framework for Fully Affine Invariant Image // Comparison, SIAM Journal on Imaging Sciences, vol. 2, issue 2, pp. 438-469, 2009. // Reference: ASIFT online demo (You can try ASIFT with your own images online.) // http://www.ipol.im/pub/algo/my_affine_sift/ /*---------------------------------------------------------------------------*/ #include #include #include #include #include #include "compute_asift_keypoints.h" #ifdef _OPENMP #include #endif #define ABS(x) (((x) > 0) ? (x) : (-(x))) /* InitSigma gives the amount of smoothing applied to the image at the first level of each octave. In effect, this determines the sampling needed in the image domain relative to amount of smoothing. Good values determined experimentally are in the range 1.2 to 1.8. */ /* float InitSigma_aa = 1.0;*/ static float InitSigma_aa = 1.6; #define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5)) /* Gaussian convolution kernels are truncated at this many sigmas from the center. While it is more efficient to keep this value small, experiments show that for consistent scale-space analysis it needs a value of about 3.0, at which point the Gaussian has fallen to only 1% of its central value. A value of 2.0 greatly reduces keypoint consistency, and a value of 4.0 is better than 3.0. */ const float GaussTruncate1 = 4.0; /* --------------------------- Blur image --------------------------- */ /* Same as ConvBuffer, but implemented with loop unrolling for increased speed. This is the most time intensive routine in keypoint detection, so deserves careful attention to efficiency. Loop unrolling simply sums 5 multiplications at a time to allow the compiler to schedule operations better and avoid loop overhead. This almost triples speed of previous version on a Pentium with gcc. */ void ConvBufferFast(float *buffer, float *kernel, int rsize, int ksize) { int i; float *bp, *kp, *endkp; float sum; for (i = 0; i < rsize; i++) { sum = 0.0; bp = &buffer[i]; kp = &kernel[0]; endkp = &kernel[ksize]; /* Loop unrolling: do 5 multiplications at a time. */ // while (kp + 4 < endkp) { // sum += (double) bp[0] * (double) kp[0] + (double) bp[1] * (double) kp[1] + (double) bp[2] * (double) kp[2] + // (double) bp[3] * (double) kp[3] + (double) bp[4] * (double) kp[4]; // bp += 5; // kp += 5; // } // /* Do 2 multiplications at a time on remaining items. */ // while (kp + 1 < endkp) { // sum += (double) bp[0] * (double) kp[0] + (double) bp[1] * (double) kp[1]; // bp += 2; // kp += 2; // } // /* Finish last one if needed. */ // if (kp < endkp) { // sum += (double) *bp * (double) *kp; // } while (kp < endkp) { sum += *bp++ * *kp++; } buffer[i] = sum; } } /* Convolve image with the 1-D kernel vector along image rows. This is designed to be as efficient as possible. Pixels outside the image are set to the value of the closest image pixel. */ void ConvHorizontal(vector& image, int width, int height, float *kernel, int ksize) { int rows, cols, r, c, i, halfsize; float buffer[4000]; vector pixels(width*height); rows = height; cols = width; halfsize = ksize / 2; pixels = image; assert(cols + ksize < 4000); for (r = 0; r < rows; r++) { /* Copy the row into buffer with pixels at ends replicated for half the mask size. This avoids need to check for ends within inner loop. */ for (i = 0; i < halfsize; i++) buffer[i] = pixels[r*cols]; for (i = 0; i < cols; i++) buffer[halfsize + i] = pixels[r*cols+i]; for (i = 0; i < halfsize; i++) buffer[halfsize + cols + i] = pixels[r*cols+cols-1]; ConvBufferFast(buffer, kernel, cols, ksize); for (c = 0; c < cols; c++) pixels[r*cols+c] = buffer[c]; } image = pixels; } /* Same as ConvHorizontal, but apply to vertical columns of image. */ void ConvVertical(vector& image, int width, int height, float *kernel, int ksize) { int rows, cols, r, c, i, halfsize; float buffer[4000]; vector pixels(width*height); rows = height; cols = width; halfsize = ksize / 2; pixels = image; assert(rows + ksize < 4000); for (c = 0; c < cols; c++) { for (i = 0; i < halfsize; i++) buffer[i] = pixels[c]; for (i = 0; i < rows; i++) buffer[halfsize + i] = pixels[i*cols+c]; for (i = 0; i < halfsize; i++) buffer[halfsize + rows + i] = pixels[(rows - 1)*cols+c]; ConvBufferFast(buffer, kernel, rows, ksize); for (r = 0; r < rows; r++) pixels[r*cols+c] = buffer[r]; } image = pixels; } /* 1D Convolve image with a Gaussian of width sigma and store result back in image. This routine creates the Gaussian kernel, and then applies it in horizontal (flag_dir=0) OR vertical directions (flag_dir!=0). */ void GaussianBlur1D(vector& image, int width, int height, float sigma, int flag_dir) { float x, kernel[100], sum = 0.0; int ksize, i; /* The Gaussian kernel is truncated at GaussTruncate sigmas from center. The kernel size should be odd. */ ksize = (int)(2.0 * GaussTruncate1 * sigma + 1.0); ksize = MAX(3, ksize); /* Kernel must be at least 3. */ if (ksize % 2 == 0) /* Make kernel size odd. */ ksize++; assert(ksize < 100); /* Fill in kernel values. */ for (i = 0; i <= ksize; i++) { x = i - ksize / 2; kernel[i] = exp(- x * x / (2.0 * sigma * sigma)); sum += kernel[i]; } /* Normalize kernel values to sum to 1.0. */ for (i = 0; i < ksize; i++) kernel[i] /= sum; if (flag_dir == 0) { ConvHorizontal(image, width, height, kernel, ksize); } else { ConvVertical(image, width, height, kernel, ksize); } } void compensate_affine_coor1(float *x0, float *y0, int w1, int h1, float t1, float t2, float Rtheta) { float x_ori, y_ori; float x_tmp, y_tmp; float x1 = *x0; float y1 = *y0; Rtheta = Rtheta*PI/180; if ( Rtheta <= PI/2 ) { x_ori = 0; y_ori = w1 * sin(Rtheta) / t1; } else { x_ori = -w1 * cos(Rtheta) / t2; y_ori = ( w1 * sin(Rtheta) + h1 * sin(Rtheta-PI/2) ) / t1; } float sin_Rtheta = sin(Rtheta); float cos_Rtheta = cos(Rtheta); /* project the coordinates of im1 to original image before tilt-rotation transform */ /* Get the coordinates with respect to the 'origin' of the original image before transform */ x1 = x1 - x_ori; y1 = y1 - y_ori; /* Invert tilt */ x1 = x1 * t2; y1 = y1 * t1; /* Invert rotation (Note that the y direction (vertical) is inverse to the usual concention. Hence Rtheta instead of -Rtheta to inverse the rotation.) */ x_tmp = cos_Rtheta*x1 - sin_Rtheta*y1; y_tmp = sin_Rtheta*x1 + cos_Rtheta*y1; x1 = x_tmp; y1 = y_tmp; *x0 = x1; *y0 = y1; } /* -------------- MAIN FUNCTION ---------------------- */ int compute_asift_keypoints(vector& image, int width, int height, int num_of_tilts, int verb, vector< vector< keypointslist > >& keys_all, siftPar &siftparameters) // Compute ASIFT keypoints in the input image. // Input: // image: input image // width, height: width and height of the input image. // num_of_tilts: number of tilts to simulate. // verb: 1/0 --> show/don not show verbose messages. (1 for debugging) // keys_all (output): ASIFT keypoints. It is a 2D matrix with varying rows and columns. Each entry keys_all[tt][rr] // stores the SIFT keypoints calculated on the image with the simulated tilt index tt and simulated rotation index rr (see the code below). In the coordinates of the keypoints, // the affine distortions have been compensated. // siftparameters: SIFT parameters. // // Output: the number of keypoints { vector image_t, image_tmp1, image_tmp; float t_min, t_k; int num_tilt, tt, num_rot_t2, rr; int fproj_o; float fproj_p, fproj_bg; char fproj_i; float *fproj_x4, *fproj_y4; // float frot_b=0; float frot_b=128; char *frot_k; int counter_sim=0, num_sim; int flag_dir = 1; float BorderFact=6*sqrt(2.); int num_keys_total=0; fproj_o = 3; fproj_p = 0; fproj_i = 0; fproj_bg = 0; fproj_x4 = 0; fproj_y4 = 0; frot_k = 0; num_rot_t2 = 10; t_min = 1; t_k = sqrt(2.); num_tilt = num_of_tilts; if ( num_tilt < 1) { printf("Number of tilts num_tilt should be equal or larger than 1. \n"); exit(-1); } image_tmp1 = image; /* Calculate the number of simulations, and initialize keys_all */ keys_all = std::vector< vector< keypointslist > >(num_tilt); for (tt = 1; tt <= num_tilt; tt++) { float t = t_min * pow(t_k, tt-1); if ( t == 1 ) { counter_sim ++; keys_all[tt-1] = std::vector< keypointslist >(1); } else { int num_rot1 = round(num_rot_t2*t/2); if ( num_rot1%2 == 1 ) { num_rot1 = num_rot1 + 1; } num_rot1 = num_rot1 / 2; counter_sim += num_rot1; keys_all[tt-1] = std::vector< keypointslist >(num_rot1); } } num_sim = counter_sim; if ( verb ) { printf("%d affine simulations will be performed. \n", num_sim); } counter_sim = 0; /* Affine simulation (rotation+tilt simulation) */ // Loop on tilts. #ifdef _OPENMP omp_set_nested(1); #endif #pragma omp parallel for private(tt) for (tt = 1; tt <= num_tilt; tt++) { float t = t_min * pow(t_k, tt-1); float t1 = 1; float t2 = 1/t; // If tilt t = 1, do not simulate rotation. if ( t == 1 ) { // copy the image from vector to array as compute_sift_keypoints uses only array. float *image_tmp1_float = new float[width*height]; for (int cc = 0; cc < width*height; cc++) image_tmp1_float[cc] = image_tmp1[cc]; compute_sift_keypoints(image_tmp1_float,keys_all[tt-1][0],width,height,siftparameters); delete[] image_tmp1_float; } else { // The number of rotations to simulate under the current tilt. int num_rot1 = round(num_rot_t2*t/2); if ( num_rot1%2 == 1 ) { num_rot1 = num_rot1 + 1; } num_rot1 = num_rot1 / 2; float delta_theta = PI/num_rot1; // Loop on rotations. #pragma omp parallel for private(rr) for ( int rr = 1; rr <= num_rot1; rr++ ) { float theta = delta_theta * (rr-1); theta = theta * 180 / PI; vector image_t; int width_r, height_r; // simulate a rotation: rotate the image with an angle theta. (the outside of the rotated image are padded with the value frot_b) frot(image, image_t, width, height, &width_r, &height_r, &theta, &frot_b , frot_k); /* Tilt */ int width_t = (int) (width_r * t1); int height_t = (int) (height_r * t2); int fproj_sx = width_t; int fproj_sy = height_t; float fproj_x1 = 0; float fproj_y1 = 0; float fproj_x2 = width_t; float fproj_y2 = 0; float fproj_x3 = 0; float fproj_y3 = height_t; /* Anti-aliasing filtering along vertical direction */ /* sigma_aa = InitSigma_aa * log2(t);*/ float sigma_aa = InitSigma_aa * t / 2; GaussianBlur1D(image_t,width_r,height_r,sigma_aa,flag_dir); // simulate a tilt: subsample the image along the vertical axis by a factor of t. vector image_tmp(width_t*height_t); fproj (image_t, image_tmp, width_r, height_r, &fproj_sx, &fproj_sy, &fproj_bg, &fproj_o, &fproj_p, &fproj_i , fproj_x1 , fproj_y1 , fproj_x2 , fproj_y2 , fproj_x3 , fproj_y3, fproj_x4, fproj_y4); vector image_tmp1 = image_tmp; if ( verb ) { printf("Rotation theta = %.2f, Tilt t = %.2f. w=%d, h=%d, sigma_aa=%.2f, \n", theta, t, width_t, height_t, sigma_aa); } float *image_tmp1_float = new float[width_t*height_t]; for (int cc = 0; cc < width_t*height_t; cc++) image_tmp1_float[cc] = image_tmp1[cc]; // compute SIFT keypoints on simulated image. keypointslist keypoints; keypointslist keypoints_filtered; compute_sift_keypoints(image_tmp1_float,keypoints,width_t,height_t,siftparameters); delete[] image_tmp1_float; /* check if the keypoint is located on the boundary of the parallelogram (i.e., the boundary of the distorted input image). If so, remove it to avoid boundary artifacts. */ if ( keypoints.size() != 0 ) { for ( int cc = 0; cc < (int) keypoints.size(); cc++ ) { float x0, y0, x1, y1, x2, y2, x3, y3 ,x4, y4, d1, d2, d3, d4, scale1, theta1, sin_theta1, cos_theta1, BorderTh; x0 = keypoints[cc].x; y0 = keypoints[cc].y; scale1= keypoints[cc].scale; theta1 = theta * PI / 180; sin_theta1 = sin(theta1); cos_theta1 = cos(theta1); /* the coordinates of the 4 submits of the parallelogram */ if ( theta <= 90 ) { x1 = height * sin_theta1; y1 = 0; y2 = width * sin_theta1; x3 = width * cos_theta1; x4 = 0; y4 = height * cos_theta1; x2 = x1 + x3; y3 = y2 + y4; /* note that the vertical direction goes from top to bottom!!! The calculation above assumes that the vertical direction goes from the bottom to top. Thus the vertical coordinates need to be reversed!!! */ y1 = y3 - y1; y2 = y3 - y2; y4 = y3 - y4; y3 = 0; y1 = y1 * t2; y2 = y2 * t2; y3 = y3 * t2; y4 = y4 * t2; } else { y1 = -height * cos_theta1; x2 = height * sin_theta1; x3 = 0; y3 = width * sin_theta1; x4 = -width * cos_theta1; y4 = 0; x1 = x2 + x4; y2 = y1 + y3; /* note that the vertical direction goes from top to bottom!!! The calculation above assumes that the vertical direction goes from the bottom to top. Thus the vertical coordinates need to be reversed!!! */ y1 = y2 - y1; y3 = y2 - y3; y4 = y2 - y4; y2 = 0; y1 = y1 * t2; y2 = y2 * t2; y3 = y3 * t2; y4 = y4 * t2; } /* the distances from the keypoint to the 4 sides of the parallelogram */ d1 = ABS((x2-x1)*(y1-y0)-(x1-x0)*(y2-y1)) / sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); d2 = ABS((x3-x2)*(y2-y0)-(x2-x0)*(y3-y2)) / sqrt((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)); d3 = ABS((x4-x3)*(y3-y0)-(x3-x0)*(y4-y3)) / sqrt((x4-x3)*(x4-x3)+(y4-y3)*(y4-y3)); d4 = ABS((x1-x4)*(y4-y0)-(x4-x0)*(y1-y4)) / sqrt((x1-x4)*(x1-x4)+(y1-y4)*(y1-y4)); BorderTh = BorderFact*scale1; if (!((d1