jift.cpp

SIFT的c++实现

                             middle(xi+1,yi-1) -
                             middle(xi-1,yi+1))/4.0 ;
                        
                            // Compute the trace and the determinant of the Hessian
                        trH  = dxx + dyy;
                        detH = dxx*dyy - dxy*dxy;
                        
                            // Compute the ratio of the principal curvatures
                        curvature_ratio = (trH*trH)/detH;
                        if (detH<0 || curvature_ratio>curvature_threshold)
                        {
                            m_extrema[i][j-1](xi,yi) = false;
                            num --;
                        }
                        
                    }
                }
        }
    }

    m_kpnum = num;
    cout << "Number of keypoints detected: " << m_kpnum << endl;
}



void JIFT::assignOrientation(void)
{
        // This function is to assign orientation to the keypoints.
        // For this, we histogram the gradient orientation over a region about
        // each keypoints
    
    
    cout << endl << "Assigning orientation to keypoints..." << endl;
    
    unsigned int i,j,k,xi,yi;
    int kk, tt;
    
    vil_image_view<double> ** magnitude = new vil_image_view<double> * [m_octaves];
    for (i=0;i<m_octaves;i++)
        magnitude[i] = new vil_image_view<double>[m_intervals];
    vil_image_view<double> ** orientation = new vil_image_view<double> * [m_octaves];
    for (i=0;i<m_octaves;i++)
        orientation[i] = new vil_image_view<double>[m_intervals];
    

    
        // Compute the gradient direction and magnitude of
        // the gaussian pyramid images
    for (i=0;i<m_octaves;i++)
        for (j=1;j<m_intervals+1;j++)
        {

            magnitude[i][j-1].set_size(m_glist[i][j].ni(),m_glist[i][j].nj());
            orientation[i][j-1].set_size(m_glist[i][j].ni(),m_glist[i][j].nj());
                // fill the view with value 0
            magnitude[i][j-1].fill(0.0);
            orientation[i][j-1].fill(0.0);
            
            
            
                // fetch the gaussian pyramid image
            for (xi=1; xi<m_glist[i][j].ni()-1; xi++)
                for (yi=1; yi<m_glist[i][j].nj()-1;yi++)
                {
                        // compute x and y derivatives using pixel differences
                    double dx = m_glist[i][j](xi+1,yi) - m_glist[i][j](xi-1,yi);
                    double dy = m_glist[i][j](xi,yi+1) - m_glist[i][j](xi,yi-1);

                        // compute the magnitude and orientation of the gradient
                    magnitude[i][j-1](xi,yi)   = sqrt(dx * dx + dy * dy);
                    orientation[i][j-1](xi,yi) = (atan2(dy,dx)==M_PI)? -M_PI:atan2(dy,dx); 
                }
                // zero padding, assign the boundary to 0.0

                /*
                  for (xi=0;xi<magnitude[i][j-1].ni();xi++)
                  {
                  magnitude[i][j-1](xi,0)                          = 0.0;
                  magnitude[i][j-1](xi,magnitude[i][j-1].nj()-1)   = 0.0;
                  orientation[i][j-1](xi,0)                        = 0.0;
                  orientation[i][j-1](xi,magnitude[i][j-1].nj()-1) = 0.0;
                  }
                  
                  for (yi=0;yi<magnitude[i][j-1].nj();yi++)
                  {
                  magnitude[i][j-1](0,yi)                          = 0.0;
                  magnitude[i][j-1](magnitude[i][j-1].ni()-1,yi)   = 0.0;
                  orientation[i][j-1](0,yi)                        = 0.0;
                  orientation[i][j-1](magnitude[i][j-1].ni()-1,yi) = 0.0;
                  }
                */
            
            
        }

    
        // The next step of the algorithm is to assign orientation to the keypoints
        // that have been located.
        // This is done by looking for peaks in the histograms of gradient orientations in
        // the regions surrounding each keypoint. A keypoint may be assigned more than one
        // orientations. If this is the case, then two identical descriptors
        // are added to the database with different orientations


    vector<double>hist_orient(NUM_BINS);
    
    for (i=0;i<m_octaves;i++)
    {
        unsigned int scale = static_cast<unsigned int> (pow(2.0f,static_cast<float>(i)));
        unsigned int ni = m_glist[i][0].ni();
        unsigned int nj = m_glist[i][0].nj();
        
        for (j=1; j<m_intervals+1; j++)
        {
                // first, fetch the image from the absolute sigma
            double abs_sigma = m_abssigmas[i][j];

                // Each sample added to the histogram is weighted by its gradient magnitude
                // and by a Gaussian-weighted circular window with a sigma that is 1.5 times
                // that of the scale of the keypoints
            vil_image_view<double> weight;
            Gaussian::gaussian_filter(magnitude[i][j-1], weight, 1.5*abs_sigma);

                // get the gaussian kernel size
            unsigned int hfsz = Gaussian::getGaussianKernelSize(1.5*abs_sigma) / 2 ;
            
            
                // iterate over all the keypoints at this octave and orientation
            for (xi=0;xi<ni;xi++)
                for (yi=0;yi<nj;yi++)
                {
                        // if it is a keypoint
                    if (m_extrema[i][j-1](xi,yi)==true)
                    {
                            // histogram the gradient orientations for this keypoint weighted by
                            // the gradient magnitude and the gaussian weighting mask

                            // clear, set the histogram to all 0.0 
                        for (k=0;k<NUM_BINS;k++)
                            hist_orient[k] = 0.0;
                        
                        for (kk=-hfsz; kk<=static_cast<int>(hfsz); kk++)
                            for (tt=-hfsz; tt<=static_cast<int>(hfsz); tt++)
                            {
                                    // if out of boundary do nothing
                                if (static_cast<int>(xi)+kk<0     ||
                                    static_cast<int>(xi)+kk>=static_cast<int>(ni)   ||
                                    static_cast<int>(yi)+tt<0     ||
                                    static_cast<int>(yi)+tt>=static_cast<int>(nj)  ) continue;
                                
                                    //fetch the orientation of sample point
                                double sampleorient = orientation[i][j-1](xi+kk,yi+tt);
                                    // note that orient[xx][yy] is within the range [-M_PI,M_PI)
                                    // should be normalised to [0, 2*M_PI)

                                assert (sampleorient >=-M_PI && sampleorient < M_PI);
                                sampleorient += M_PI;
                                
                                
                                    // convert to degree
                                unsigned sampleorient_d = static_cast<unsigned int>(sampleorient * 180 / M_PI);
                                    // accumulate the histogram bins
                                hist_orient[sampleorient_d / (360/NUM_BINS)] += weight(xi+kk,yi+tt);
                            }
                        
                            // find the maximum peak in the orientation histogram
                        double max_peak_val       = hist_orient[0];
                        unsigned int max_peak_idx = 0;
                        for (k=1; k<NUM_BINS; k++)
                            if (hist_orient[k]>max_peak_val)
                            {
                                max_peak_val = hist_orient[k];
                                max_peak_idx = k;
                            }

                            // iterate over all peaks within 80% of the largest peak and add
                            // keypoints with the orientation corresponding to those peaks
                            // to the keypoint list
                        vector<double> orien;
                        vector<double> mag;
                        
                        for (k=0; k<NUM_BINS; k++)
                            if (hist_orient[k] > 0.8 * max_peak_val)
                            {
                                    // if this is the maxmium peak
                                    // interpolate the peak by fitting a parabola to the three
                                    // histogram values closest to each peak, for better accuracy
                                    // Parabola??? oh my god ...
                                double x1 = static_cast<double>(k-1);
                                double y1;
                                double x2 = static_cast<double>(k);
                                double y2 = hist_orient[k];
                                double x3 = static_cast<double>(k+1);
                                double y3;

                                if (k==0)
                                {
                                    y1 = hist_orient[NUM_BINS-1];
                                    y3 = hist_orient[1];
                                }
                                else if (k==NUM_BINS-1)
                                {
                                    y1 = hist_orient[NUM_BINS-2];
                                    y3 = hist_orient[0];
                                }
                                else
                                {
                                    y1 = hist_orient[k-1];
                                    y3 = hist_orient[k+1];
                                }

                                    // A downward parabola has the general form

                                    // y = a * x^2 + bx + c
                                
                                    // Now the three equations stem from the three points
                                    // (x1,y1) (x2,y2) (x3.y3) are
                                
                                    // y1 = a * x1^2 + b * x1 + c
                                    // y2 = a * x2^2 + b * x2 + c
                                    // y3 = a * x3^2 + b * x3 + c
                                
                                    // in Matrix notation, this is y = Xb, where
                                    // y = (y1 y2 y3)' b = (a b c)' and
                                    // 
                                    //     x1^2 x1 1
                                    // X = x2^2 x2 1
                                    //     x3^2 x3 1
                                    //
                                    // OK, we need to solve this equation for b
                                    // this is done by inverse the matrix X
                                    //
                                    // b = inv(X) Y
                                vnl_vector<double> b(3);
                                vnl_matrix<double> X(3,3);

                                X(0,0) = x1*x1; X(0,1) = x1; X(0,2) = 1;
                                X(1,0) = x2*x2; X(1,1) = x2; X(1,2) = 1;
                                X(2,0) = x3*x3; X(2,1) = x3; X(2,2) = 1;
                                
                                X = vnl_matrix_inverse<double>(X);

                                b[0] = X(0,0)*y1 + X(0,1)*y2 + X(0,2)*y3;
                                b[1] = X(1,0)*y1 + X(1,1)*y2 + X(1,2)*y3;
                                b[2] = X(2,0)*y1 + X(2,1)*y2 + X(2,2)*y3;

                                    // the vertex of parabola is (-b/(2a),c-b^2/4a)
                                double x0 = -b[1]/(2*b[0]);

                                
                                    // abnormal situation
                                if (fabs(x0)>2*NUM_BINS)
                                    x0 = x2;
                                
                                
                                    // make sure it is within the range
                                while (x0 < 0)
                                    x0 = x0 + NUM_BINS;
                                while (x0 >= NUM_BINS)
                                    x0 = x0 - NUM_BINS;
                                
                                    // normalise
                                double x0_n = x0*(2*M_PI/NUM_BINS);
                                    // since x0_n is [0,2*M_PI) we should set it within the range [-M_PI,M_PI)
                                
                                assert (x0_n>=0 && x0_n<2*M_PI);
                                x0_n -= M_PI;
                                assert (x0_n >= -M_PI && x0_n < M_PI);
                                
                                    // store the keypoint orientations
                                    // allow multiple orientations
                                orien.push_back(x0_n);
                                mag.push_back(hist_orient[k]);
                            }
                        
                        m_kps.push_back(
                            Keypoint(static_cast<float>(xi*scale)/2.0,  // x position
                                     static_cast<float>(yi*scale)/2.0,  // y position
                                     mag,                               // magnitude
                                     orien,                             // orientation
                                     i*m_intervals+j-1)                 // scale
                            );
                    }
                    
                }
            
        }
        
    }
    
    
    assert (m_kps.size() == m_kpnum);

        // free the memory space
    for (i=0;i<m_octaves;i++)
    {
        delete [] magnitude[m_octaves-i-1];
        delete [] orientation[m_octaves-i-1];
    }
    delete [] magnitude;
    delete [] orientation;
    
    
    
}

void JIFT::extractKeypointDescriptor(void)
{
    cout << endl << "Extracting feature descriptors..." << endl;
    
        // The final of the SIFT algorithm is to extract feature descriptors for the keypoints
        // The descriptors are a grid of gradient orientation histograms, where the sampling
        // grid for the histograms is rotated to the main orientation of each keypoint.
        // The grid is a 4x4 array of 4x4 sample cells of a bin orientation histograms.
        // This produces 128 dimensional feature vectors

        // First, compute the gradient magnitude and orientation at the sample point (X.5)
        // As the sample point is not an integer number, it should be interpolated using
        // bilinear interpolation


    unsigned int i,j;
    vil_image_view<double> ** InterpolatedMagnitude   = new vil_image_view<double>*[m_octaves];
    vil_image_view<double> ** InterpolatedOrientation = new vil_image_view<double>*[m_octaves];
    for (i=0;i<m_octaves;i++)
    {
        InterpolatedMagnitude[i]   = new vil_image_view<double>[m_intervals];
        InterpolatedOrientation[i] = new vil_image_view<double>[m_intervals];
    }
    
    
    for (i=0; i<m_octaves; i++)
        for (j=1; j<m_intervals+1; j++)
        {
                // fetch the width and height of the pyramid image
            unsigned int ni = m_glist[i][j].ni();
            unsigned int nj = m_glist[i][j].nj();

                // vcl_cout << "ni = " << ni << ", nj = " << nj << vcl_endl;
            
                // note that n grids will have n+1 conjunctions
            InterpolatedMagnitude[i][j-1].set_size(ni+1,nj+1);