[Opencv] Sift principle and source code analysis: Key POint Description

Index of Articles in the sift principle and source code analysis series: Middleware. The next step is to describe the key point, that is, using a set of vectors to describe the key point. This describe includes not only the key point, but also the pixels that contribute to the key point. Used Target matching basis(So the descriptor should have a high uniqueness to ensure the matching rate.) It can also make the key points more unchanged, such as illumination change and 3D Viewpoint Change.

Sift description sub-h (x, y, θ) is the result of Gaussian image gradient statistics near the key point. It is a three-dimensional matrix, but is usually represented by a vector. Vector is obtained by regular arrangement of three-dimensional matrices.

The description subsampling area feature description is related to the scale of the key points. Therefore, the gradient calculation should be performed on the Gaussian image corresponding to the feature points. Divide the key points into sub-areas d x D, and the size of each sub-area is m σ pixel (D = 4, M = 3, σ is the scale value of the feature points ). Considering that bilinear interpolation is required for actual calculation, the calculated image region is m σ (D + 1). If rotation is considered, the actually calculated image region is, as shown in: source code

Point Pt (cvround (PTF. x), cvround (PTF. y); // calculate the cosine, sine, cv_pi/180: Convert the angle value to the amplitude value float cos_t = cosf (ORI * (float) (cv_pi/180 )); float sin_t = sinf (ORI * (float) (cv_pi/180); float bins_per_rad = N/360.f; float exp_scale =-1.f/ (D * 0.5f); // D: sift_descr_width 4 float hist_width = sift_descr_scl_fctr * SCL; // sift_descr_scl_fctr: 3 // SCL: size * 0.5f // calculate the image region radius M σ (D + 1) /2 * SQRT (2) // 1.4142135623730951f is the root number 2 int radius = cvround (hist_width * 1.4142135623730951f * (D + 1) * 0.5f); cos_t/= hist_width; sin_t/= hist_width;

In order to ensure the feature vectors have rotation immutability, the coordinate axis rotation must take the feature point as the center and rotate the θ angle in the nearby area, that is, the rotation is the direction of the feature point. The new coordinates of the sampling points in the area after rotation are: source code

// Calculate the coordinate rotation for (I =-radius, K = 0; I <= radius; I ++) for (j =-radius; j <= radius; j ++) {/* calculate sample's histogram array coords rotated relative to Ori. subtract 0.5 so samples that fall e.g. in the center of Row 1 (I. e. r_rot = 1.5) have full weight placed in Row 1 after interpolation. */float c_rot = J * cos_t-I * sin_t; float r_rot = J * sin_t + I * cos_t; float rbin = r_rot + D/2-0.5f; float cbin = c_rot + D/2-0.5f; int r = pt. Y + I, C = pt. X + J; if (rbin>-1 & rbin <D & cbin>-1 & cbin <D & R> 0 & R <rows-1 & C> 0 & & C <Cols-1) {float dx = (float) (IMG. at <short> (R, C + 1)-IMG. at <short> (R, C-1); float DY = (float) (IMG. at <short> (R-1, c)-IMG. at <short> (R + 1, C); X [k] = DX; y [k] = Dy; rbin [k] = rbin; cbin [k] = cbin; W [k] = (c_rot * c_rot + r_rot * r_rot) * exp_scale; k ++ ;}}

Calculate the gradient histogram of the sampling area and divide the area after rotation into D × d sub-areas (the interval of each area is m σ pixel). Calculate the gradient histogram of eight directions in the sub-areas, draw the cumulative value of each direction gradient to form a seed point. The difference is that, in this case, the histogram of the gradient direction of each subarea is 0 ° ~ 360 ° is divided into eight direction intervals, each of which is 45 °. That is, each seed point has eight gradient intensity intervals. Because d x D exists, that is, 4x4 subareas, a total of 4x4x8 = 128 data records are generated to form a 128-dimensional sift feature vector. The feature vectors need to be weighted, and mσ D/2 standard Gaussian function is used for weighting. In addition to the effect of illumination changes, we also need to normalize the image. Source code

// Calculate the gradient histogram for (k = 0; k <Len; k ++) {float rbin = rbin [K], cbin = cbin [k]; float obin = (ORI [k]-Ori) * bins_per_rad; float mag = mag [k] * W [k]; int R0 = cvfloor (rbin ); int C0 = cvfloor (cbin); int O0 = cvfloor (obin); rbin-= R0; cbin-= C0; obin-= O0; // n is sift_descr_hist_bins: 8, about 360 ° is divided into eight intervals: If (O0 <0) O0 + = N; If (O0> = N) O0-= N; // histogram update using tri-linear interpolation // bilinear interpolation float v_r1 = mag * rbin, v_r0 = mag-v_r1; float v_rc11 = v_r1 * cbin, v_rc10 = v_r1-v_rc11; float v_rc01 = v_r0 * cbin, v_rc00 = v_r0-v_rc01; float ranges = v_rc11 * obin, v_rco110 = v_rc11-rows; float ranges = v_rc10 * obin, v_rco100 = v_rc10-v_rco101; float v_rco011 = v_rc01 * obin, v_rco010 = v_rc01-v_rco011; float v_rco001 = v_rc00 * obin, v_rco000 = v_rc00-v_rco001; int idx = (R0 + 1) * (D + 2) + C0 + 1) * (N + 2) + O0; hist [idx] + = v_rco000; hist [idx + 1] + = v_rco001; hist [idx + (n + 2)] + = v_rco010; hist [idx + (n + 3)] + = v_rco011; hist [idx + (D + 2) * (N + 2)] + = v_rco100; hist [idx + (D + 2) * (N + 2) + 1] + = v_rco101; hist [idx + (D + 3) * (N + 2)] + = v_rco110; hist [idx + (D + 3) * (N + 2) + 1] + = v_rco111 ;}

Key description source code

// Sift key point feature description // sift descriptor is a form of static void calcsiftdescriptor (const mat & IMG, point2f Ptf, float Ori, float SCL, int D, int N, float * DST) {point Pt (cvround (PTF. x), cvround (PTF. y); // calculate the cosine, sine, cv_pi/180: Convert the angle value to the amplitude value float cos_t = cosf (ORI * (float) (cv_pi/180 )); float sin_t = sinf (ORI * (float) (cv_pi/180); float bins_per_rad = N/360.f; float exp_scale =-1.f/ (D * 0.5f); // D: sift_descr_width 4 float hist_width = sift_descr_scl_fctr * SCL; // sift_descr_scl_fctr: 3 // SCL: size * 0.5f // calculate the image region radius M σ (D + 1) /2 * SQRT (2) // 1.4142135623730951f is the root number 2 int radius = cvround (hist_width * 1.4142135623730951f * (D + 1) * 0.5f); cos_t/= hist_width; sin_t/= hist_width; int I, J, K, Len = (radius * 2 + 1) * (radius * 2 + 1), histlen = (D + 2) * (D + 2) * (N + 2); int rows = IMG. rows, cols = IMG. cols; autobuffer <float> Buf (LEN * 6 + histlen); float * x = Buf, * Y = x + Len, * mag = Y, * Ori = Mag + Len, * w = Ori + Len; float * rbin = W + Len, * cbin = rbin + Len, * hist = cbin + Len; // initialize the histogram for (I = 0; I <D + 2; I ++) {for (j = 0; j <D + 2; j ++) for (k = 0; k <n + 2; k ++) hist [(I * (D + 2) + J) * (N + 2) + K] = 0 .;} // calculate the coordinate rotation for (I =-radius, K = 0; I <= radius; I ++) for (j =-radius; j <= radius; j ++) {/* calculate sample's histogram array coords rotated relative to Ori. subtract 0.5 so samples that fall e.g. in the center of Row 1 (I. e. r_rot = 1.5) have full weight placed in Row 1 after interpolation. */float c_rot = J * cos_t-I * sin_t; float r_rot = J * sin_t + I * cos_t; float rbin = r_rot + D/2-0.5f; float cbin = c_rot + D/2-0.5f; int r = pt. Y + I, C = pt. X + J; if (rbin>-1 & rbin <D & cbin>-1 & cbin <D & R> 0 & R <rows-1 & C> 0 & & C <Cols-1) {float dx = (float) (IMG. at <short> (R, C + 1)-IMG. at <short> (R, C-1); float DY = (float) (IMG. at <short> (R-1, c)-IMG. at <short> (R + 1, C); X [k] = DX; y [k] = Dy; rbin [k] = rbin; cbin [k] = cbin; W [k] = (c_rot * c_rot + r_rot * r_rot) * exp_scale; k ++;} Len = K; fastatan2 (Y, X, ori, Len, true); magnitude (X, Y, Mag, Len); exp (W, W, Len); // calculate the gradient histogram for (k = 0; k <Len; k ++) {float rbin = rbin [K], cbin = cbin [k]; float obin = (ORI [k]-Ori) * bins_per_rad; float mag = mag [k] * W [k]; int R0 = cvfloor (rbin); int C0 = cvfloor (cbin); int O0 = cvfloor (obin ); rbin-= R0; cbin-= C0; obin-= O0; // n is sift_descr_hist_bins: 8, that is, 360 ° is divided into eight intervals if (O0 <0) o0 + = N; If (O0> = N) O0-= N; // histogram update using tri-linear interpolation // bilinear interpolation float v_r1 = mag * rbin, v_r0 = mag-v_r1; float v_rc11 = v_r1 * cbin, v_rc10 = v_r1-v_rc11; float v_rc01 = v_r0 * cbin, v_rc00 = v_r0-v_rc01; float v_rco111 = v_rc11 * obin, v_rco110 = v_rc11-rows; float v_rco101 = v_rc10 * obin, rows = v_rc10-v_rco101; float rows = v_rc01 * obin, rows = v_rc01-rows; float v_rco001 = v_rc00 * obin, v_rco000 = v_rc00-v_rco001; int idx = (R0 + 1) * (D + 2) + C0 + 1) * (N + 2) + O0; hist [idx] + = v_rco000; hist [idx + 1] + = v_rco001; hist [idx + (n + 2)] + = v_rco010; hist [idx + (n + 3)] + = v_rco011; hist [idx + (D + 2) * (N + 2)] + = v_rco100; hist [idx + (D + 2) * (N + 2) + 1] + = v_rco101; hist [idx + (D + 3) * (N + 2)] + = v_rco110; hist [idx + (D + 3) * (N + 2) + 1] + = v_rco111;} // finalize histogram, since the orientation histograms are circular // finally determine the histogram. The target direction histogram is circular for (I = 0; I <D; I ++) for (j = 0; j <D; j ++) {int idx = (I + 1) * (D + 2) + (J + 1) * (N + 2 ); hist [idx] + = Hist [idx + N]; hist [idx + 1] + = Hist [idx + n + 1]; for (k = 0; k <N; k ++) DST [(I * D + J) * n + k] = Hist [idx + k];} // copy histogram to the descriptor, // apply hysteresis thresholding // and scale the result, so that it can be easily converted // to byte array float nrm2 = 0; Len = D * N; for (k = 0; k <Len; k ++) nrm2 + = DST [k] * DST [k]; float thr = STD: SQRT (nrm2) * sift_descr_mag_thr; for (I = 0, nrm2 = 0; I <K; I ++) {float val = STD: min (DST [I], Thr ); DST [I] = val; nrm2 + = Val * val;} nrm2 = sift_int_descr_fctr/STD: max (STD: SQRT (nrm2), flt_epsilon ); for (k = 0; k <Len; k ++) {DST [k] = saturate_cast <uchar> (DST [k] * nrm2 );}}

Now, the sift description sub-generation is complete ~ See "sift principle and source code analysis" (reprinted please indicate the author and Source: http://blog.csdn.net/xiaowei_cqu is not allowed for commercial purposes)