jift.cpp

SIFT的c++实现

            InterpolatedOrientation[i][j-1].set_size(ni+1,nj+1);

                // starting from 0.5
            for (float ii=0.5;ii<ni-0.5;ii++)
                for(float jj=0.5;jj<nj-0.5;jj++)
                {
                        //compute x and y derivatives using pixel differences
                    double dx =  vil_bilin_interp_safe_extend(m_glist[i][j],ii+1.0,jj)
                        - vil_bilin_interp_safe_extend(m_glist[i][j],ii-1.0,jj);
                    double dy =  vil_bilin_interp_safe_extend(m_glist[i][j],ii,jj+1.0)
                        - vil_bilin_interp_safe_extend(m_glist[i][j],ii,jj-1.0);
                    unsigned int iii = static_cast<unsigned int>(ii+1.0);
                    unsigned int jjj = static_cast<unsigned int>(jj+1.0);
                    
                        // vcl_cout << "iii = " << iii << ", jjj = " << jjj << vcl_endl;
                    assert (iii <= ni && jjj <= nj);
                    
                    InterpolatedMagnitude[i][j-1](iii,jjj)   = sqrt(dx*dx + dy*dy);
                        // vcl_cout << InterpolatedMagnitude[i][j-1](iii,jjj) << vcl_endl;
                    InterpolatedOrientation[i][j-1](iii,jjj) = (atan2(dy,dx)==M_PI)? -M_PI:atan2(dy,dx);
                        // vcl_cout << InterpolatedOrientation[i][j-1](iii,jjj) << vcl_endl;
                }
            
                // zero padding, assign the boundary to 0.0
            for(unsigned int iii=0;iii<ni+1;iii++)
            {
                InterpolatedMagnitude[i][j-1](iii,0)      = 0.0;
                InterpolatedMagnitude[i][j-1](iii,nj)     = 0.0;
                InterpolatedOrientation[i][j-1](iii,0)    = 0.0;
                InterpolatedOrientation[i][j-1](iii,nj)   = 0.0;
            }
            for(unsigned int jjj=0; jjj<nj+1; jjj++)
            {
                InterpolatedMagnitude[i][j-1](0,jjj)      = 0.0;
                InterpolatedMagnitude[i][j-1](ni,jjj)     = 0.0;
                InterpolatedOrientation[i][j-1](0,jjj)    = 0.0;
                InterpolatedOrientation[i][j-1](ni,jjj)   = 0.0;
            }
            
        }
    
        // Then, build a Interpolated gaussian table of size FEATURE_WINDOW_SIZE
        // Lowe suggest the sigma to be one half of the window width
    vnl_matrix<double> G = Gaussian::buildInterpolatedGaussianTable(FEATURE_WINDOW_SIZE,0.5*FEATURE_WINDOW_SIZE);

        // the histogram
    vector<double> hist(DESC_NUM_BINS);

    
        // Loop over all the keypoints
    for (unsigned int ikp = 0; ikp<m_kpnum; ikp++)
    {
            // vcl_cout << "ikp = " << ikp << vcl_endl;
        
        unsigned int scale = this->m_kps[ikp].scale;
        float kpxi         = this->m_kps[ikp].xi;
        float kpyi         = this->m_kps[ikp].yi;

            // the descriptor's position information
        float descxi       = kpxi;
        float descyi       = kpyi;
        
        unsigned int ii = static_cast<unsigned int>(kpxi*2) /
            static_cast<unsigned int>(pow(2.0f,static_cast<float>(scale/m_intervals)));
        unsigned int jj = static_cast<unsigned int>(kpyi*2) /
            static_cast<unsigned int>(pow(2.0f,static_cast<float>(scale/m_intervals)));
        
        unsigned int ni = m_glist[scale/m_intervals][0].ni();
        unsigned int nj = m_glist[scale/m_intervals][0].nj();

            // In order to achieve orientation invariance, the coordinates of the descriptor and the
            // gradient orientations are rotated relative to the keypoint orientation.
            // find the main orientation,
            // which is the orientation corresponding to the largest magnitude
        vector<double> orien = this->m_kps[ikp].orien;
        vector<double> mag   = this->m_kps[ikp].mag;
        
        double main_mag = mag[0];
        double main_orien = orien[0];

        for (unsigned int orient_count=1; orient_count<mag.size(); orient_count++)
            if (mag[orient_count] > main_mag)
            {
                main_orien = orien[orient_count];
                main_mag = mag[orient_count];
            }

            // vcl_cout << "main_orient = " << main_orien << vcl_endl;
            // vcl_cout << "main_mag    = " << main_mag    << vcl_endl;
        
        
            // sample points in FEATURE_WINDOW_SIZE x FEATURE_WINDOW_SIZE region 
            // assign a weight to the magnitude of each sample point
        unsigned int hfsz = FEATURE_WINDOW_SIZE/2;
        vnl_matrix<double> weight(FEATURE_WINDOW_SIZE,FEATURE_WINDOW_SIZE);
        vcl_vector<double> fv(FVSIZE);

        for(i=0;i<FEATURE_WINDOW_SIZE;i++)
            for(j=0;j<FEATURE_WINDOW_SIZE;j++)
            {
                    // vcl_cout << "ii = "   << ii   << ", "
                    //     << "i = "    << i    << ", "
                    //     << "ni = "   << ni   << ", "
                    //     << "jj = "   << jj   << ", "
                    //     << "j ="     << j    << ", "
                    //     << "nj = "   << nj   << ", "
                    //    << "hfsz = " << hfsz << endl;

                if(ii+i+1<hfsz || ii+i+1>ni+hfsz || jj+j+1<hfsz || jj+j+1>nj+hfsz)
                    weight[i][j] = 0;
                else
                {
                        // vcl_cout << scale/m_intervals << " " << scale%m_intervals << vcl_endl;
                        //  vcl_cout << ii +i+1 - hfsz << " " << jj+j+1-hfsz << vcl_endl;

                        // vcl_cout << InterpolatedMagnitude[scale/m_intervals][scale%m_intervals].ni() << " "
                        //     << InterpolatedMagnitude[scale/m_intervals][scale%m_intervals].nj() << " "
                        //     << ii+i+1-hfsz << " "
                        //     << jj+j+1-hfsz << vcl_endl;
                    
                    
                    weight[i][j] = G[i][j] *
                        InterpolatedMagnitude[scale/m_intervals][scale%m_intervals](ii+i+1-hfsz,jj+j+1-hfsz);
                    
                    
                }
                
                
                
//                vcl_cout << "G["      << i << "][" << j << "] = " << G[i][j]      << vcl_endl;
//                vcl_cout << "weight[" << i << "][" << j << "] = " << weight[i][j] << vcl_endl;
            }
        
        
        
        
            // OK, now generate the 4x4x8 = 128 feature vector for each keypoint
            // the 8-feature vectors are in the region
        
            // (kprow-hfsz+1, kpcol-hfsz+1) -- (kprow-hfsz/2, kpcol-hfsz/2)
            // (kprow-hfsz/2+1,kpcol-hfsz+1) --  (kprow, kpcol-hfsz/2)
            // (kprow+1,kpcol-hfsz+1) -- (kprow+hfsz/2,kpcol-hfsz/2)
            // (kprow+hfsz/2+1,kpcol-hfsz+1) -- (kprow+hfsz,kpcol-hfsz/2)
        
            // (kprow-hfsz+1, kpcol-hfsz/2+1) -- (kprow-hfsz/2, kpcol)
            // (kprow-hfsz/2+1,kpcol-hfsz/2+1) --  (kprow, kpcol)
            // (kprow+1,kpcol-hfsz/2+1) -- (kprow+hfsz/2,kpcol)
            // (kprow+hfsz/2+1,kpcol-hfsz/2+1) -- (kprow+hfsz,kpcol)

            // (kprow-hfsz+1, kpcol+1) -- (kprow-hfsz/2, kpcol+hfsz/2)
            // (kprow-hfsz/2+1,kpcol+1) --  (kprow, kpcol+hfsz/2)
            // (kprow+1,kpcol+1) -- (kprow+hfsz/2,kpcol+hfsz/2)
            // (kprow+hfsz/2+1,kpcol+1) -- (kprow+hfsz,kpcol+hfsz/2)

            // (kprow-hfsz+1, kpcol+hfsz/2+1) -- (kprow-hfsz/2, kpcol+hfsz)
            // (kprow-hfsz/2+1,kpcol+hfsz/2+1) --  (kprow, kpcol+hfsz)
            // (kprow+1,kpcol+hfsz/2+1) -- (kprow+hfsz/2,kpcol+hfsz)
            // (kprow+hfsz/2+1,kpcol+hfsz/2+1) -- (kprow+hfsz,kpcol+hfsz)

            // so in general for i=0..3 j=0..3
        
            // (kprow-hfsz+1+(hfsz/2)*i, kpcol-hfsz+1+(hfsz/2)*j) -- (kprow+(hfsz/2)*(i-1), kpcol+(hfsz/2)*(j-1))

//        vcl_cout << "start" << vcl_endl;
        
        for (i=0;i<FEATURE_WINDOW_SIZE/4;i++)
            for(j=0;j<FEATURE_WINDOW_SIZE/4;j++)
            {
                    // clear the histogram
                for(unsigned int t=0;t<DESC_NUM_BINS;t++)
                    hist[t] = 0.0;
                
                    // find the start and the limit
                int starti = static_cast<int>(ii)-static_cast<int>(hfsz)+1+static_cast<int>((hfsz/2)*i);
                int startj = static_cast<int>(jj)-static_cast<int>(hfsz)+1+static_cast<int>((hfsz/2)*j);
                int limiti = static_cast<int>(ii)+static_cast<int>(hfsz/2)*(static_cast<int>(i)-1);
                int limitj = static_cast<int>(jj)+static_cast<int>(hfsz/2)*(static_cast<int>(j)-1);

                    // vcl_cout << "starti = " << starti << ", "
                    //     << "startj = " << startj << ", "
                    //     << "limiti = " << limiti << ", "
                    //     << "limitj = " << limitj << vcl_endl;
                    // loop over the region
                for(int k=starti;k<=limiti;k++)
                    for(int t=startj;t<=limitj;t++)
                    {
                            // if out of boundary
                        if(k<0 || k>static_cast<int>(ni) || t<0 || t>static_cast<int>(nj))
                            continue;
                            // fetch the sample point orientation
                        double sample_orien = InterpolatedOrientation[scale/m_intervals][scale%m_intervals](k,t);
                            // rotate relative to keypoint orientation
                        sample_orien -= main_orien;
                        
                            // as sample_orient is within [-M_PI,M_PI) and main_orient is also within [-M_PI,M_PI)
                            // the subtraction will result in (-2*M_PI,2*M_PI)
                            // so we need to normalise it to [0,2*M_PI)
                        if(sample_orien<0)
                            sample_orien += 2*M_PI;

                            // vcl_cout << "sample_orien = " << sample_orien << vcl_endl;
                        
                        assert (sample_orien >= 0 && sample_orien < 2*M_PI);
                        
                        
                            // now convert into degree [0,360)
                        unsigned int sample_orien_d = static_cast<unsigned int>(sample_orien*180/M_PI);

                        assert (sample_orien_d<360);
                        
                        
                            // bin represent the bin index
                            // while bin_f is the actual entry
                        unsigned int bin = sample_orien_d/(360/DESC_NUM_BINS);
                        double bin_f = static_cast<double>(sample_orien_d)/static_cast<double>(360/DESC_NUM_BINS);
                        
                        
                            // To avoid boundary effects, each sample is accumulated into neighbouring bins
                            // weighted by 1-d in all dimensions, where d is the distance from the center of
                            // bin measured in units of bin spacing
                        
                            // the central value of bin entry is (bin+0.5)
                            // so the distance of sample from the central value of the bin is
                            // d = fabs(bin_f-(bin+0.5))
                            // so the weight 1-d is
                            // w = 1 - fabs(bin_f-(bin+0.5))
                        
                            // the weight of (k,t) correspond to weight[?][?]
                            // index of (k,t) corresponding to keypoint position is
                            // (k-kprow, t-kpcol)
                            // which should correspond to weight[k-kprow+hfsz-1][t-kpcol+hfsz-1]
                        
                            // then accumulate into the histogram
                            //  vcl_cout << k-ii+hfsz-1 << " " << t-jj+hfsz-1 << vcl_endl;

                        assert (bin < DESC_NUM_BINS);
                        assert (k+hfsz-1-ii<FEATURE_WINDOW_SIZE && t+hfsz-1-jj < FEATURE_WINDOW_SIZE);
                        
                        
                        hist[bin] += (1-fabs(bin_f-(bin+0.5))) * weight[k+hfsz-1-ii][t+hfsz-1-jj];
                    }
                    // assign the histogram to feature vector
                    // the index of feature vector is (i*FEATURE_WINDOW_SIZE+j)*DESC_NUM_BINS + bin
                
                for(unsigned int t=0;t<DESC_NUM_BINS;t++)
                {
                    fv[(i*FEATURE_WINDOW_SIZE/4+j)*DESC_NUM_BINS + t] = hist[t];
                        // vcl_cout << "hist[" << t << "] = " << hist[t] << vcl_endl;
                }
                
            }


            //        vcl_cout << "end" << vcl_endl;
        
        
            // then, normalise ...
        double norm = 0.0;
        for(unsigned int t=0;t<FVSIZE;t++)
            norm += pow(fv[t],2.0);
        norm = sqrt(norm);
        for(unsigned int t=0;t<FVSIZE;t++)
        {
            fv[t] /= norm;
                // vcl_cout << fv[t] << vcl_endl;
        }
        
        
            // for(t=0;t<FVSIZE;t++)
            // printf("ikp=%d,t=%d,fv=%f\n",ikp,t,this -> m_desc[ikp].fv[t]);

            // thresholding the vector to no more than FV_THRESHOLD 
        for(unsigned int t=0;t<FVSIZE;t++)
            if(fv[t] > FV_THRESHOLD)
                fv[t] =  FV_THRESHOLD;
        
            //renomalise
        
        norm = 0.0;
        for(unsigned int t=0;t<FVSIZE;t++)
            norm += pow(fv[t],2.0);
        norm = sqrt(norm);
        for(unsigned int t=0;t<FVSIZE;t++)
        {
            fv[t] /= norm;
                // vcl_cout << fv[t] << vcl_endl;
            
        }
        
        
            // vcl_cout << "start" << vcl_endl;
        
        
            // finally push back the descriptor
        m_desc.push_back(Descriptor(descxi, descyi, fv));

            // vcl_cout << "end" << vcl_endl;
    }
    
    
        // make sure it is OK
    assert (m_desc.size()==m_kpnum);
    

        // free the storage
    for (i=0;i<m_octaves;i++)
    {
        delete [] InterpolatedMagnitude[m_octaves-i-1];
        delete [] InterpolatedOrientation[m_octaves-i-1];
    }
    delete [] InterpolatedMagnitude;
    delete [] InterpolatedOrientation;

}



void JIFT::doJIFT(void)
{
    this -> buildPyramid();
    this -> detectLocalExtrema();
    this -> assignOrientation();
    this -> extractKeypointDescriptor();
}

const vector<Descriptor> & JIFT::descriptor() const
{
    return m_desc;
}

const vector<Keypoint>   & JIFT::keypoint() const
{
    return m_kps;
}


const vil_image_view<vxl_byte> & JIFT::srcimg()const
{
    return m_srcimg;
}