mainSift.c

00001 //********************************************************//
00002 // CUDA SIFT extractor by Marten Bjorkman aka Celebrandil //
00003 //              celle @ nada.kth.se                       //
00004 //********************************************************//
00005 
00006 #include <cutil.h>
00007 #include <iostream>
00008 #include <cmath>
00009 
00010 #include "tpimage.h"
00011 #include "tpimageutil.h"
00012 
00013 #include "cudaImage.h"
00014 #include "cudaSift.h"
00015 
00016 
00017 ///////////////////////////////////////////////////////////////////////////////
00018 // Reference CPU convolution
00019 ///////////////////////////////////////////////////////////////////////////////
00020 extern "C" void ConvRowCPU(float *h_Result, float *h_Data,
00021                            float *h_Kernel, int w, int h, int kernelR);
00022 extern "C" void ConvColCPU(float *h_Result, float *h_Data,
00023                            float *h_Kernel, int w, int h, int kernelR);
00024 extern "C" double Find3DMinMaxCPU(CudaImage *res, CudaImage *data1,
00025                                   CudaImage *data2, CudaImage *data3);
00026 
00027 ///////////////////////////////////////////////////////////////////////////////
00028 // Main program
00029 ///////////////////////////////////////////////////////////////////////////////
00030 int main(int argc, char **argv)
00031 {
00032   CudaImage img1, img2;
00033   Image<float> limg(1280, 960);
00034   Image<float> rimg(1280, 960);
00035 
00036   limg.Load("data/left2.pgm");
00037   rimg.Load("data/righ2.pgm");
00038   ReScale(limg, 1.0/256.0f);
00039   ReScale(rimg, 1.0/256.0f);
00040   unsigned int w = limg.GetWidth();
00041   unsigned int h = limg.GetHeight();
00042 
00043   std::cout << "Image size = (" << w << "," << h << ")" << std::endl;
00044 
00045   InitCuda(argc,argv);
00046   //CUT_DEVICE_INIT(argc,argv);
00047   std::cout << "Initializing data..." << std::endl;
00048   AllocCudaImage(&img1, w, h, w, false, true);
00049   img1.h_data = limg.GetData();
00050   AllocCudaImage(&img2, w, h, w, false, true);
00051   img2.h_data = rimg.GetData();
00052   Download(&img1);
00053   Download(&img2);
00054 
00055   SiftData siftData1, siftData2;
00056   InitSiftData(&siftData1, 128, true, true);
00057   ExtractSift(&siftData1, &img1, 3, 3, 0.3f, 0.04);
00058   InitSiftData(&siftData2, 128, true, true);
00059   ExtractSift(&siftData2, &img2, 3, 3, 0.3f, 0.04);
00060   std::cout << "Number of original features: " <<  siftData1.numPts << " "
00061             << siftData2.numPts << std::endl;
00062   MatchSiftData(&siftData1, &siftData2);
00063   float homography[9];
00064   int numMatches;
00065   FindHomography(&siftData1, homography, &numMatches, 1000, 0.85f, 0.95f, 5.0);
00066 
00067   int numPts = siftData1.numPts;
00068   SiftPoint *sift1 = siftData1.h_data;
00069   SiftPoint *sift2 = siftData2.h_data;
00070   float *h_img = img1.h_data;
00071   for (int j=0;j<numPts;j++) {
00072     int k = sift1[j].match;
00073     float x = sift1[j].xpos;
00074     float y = sift1[j].ypos;
00075     float den = homography[6]*x + homography[7]*y + homography[8];
00076     float x2 = (homography[0]*x + homography[1]*y + homography[2]) / den;
00077     float y2 = (homography[3]*x + homography[4]*y + homography[5]) / den;
00078     float erx = x2 - sift2[k].xpos;
00079     float ery = y2 - sift2[k].ypos;
00080     float er2 = erx*erx + ery*ery;
00081     if (er2<25.0f && sift1[j].score>0.85f && sift1[j].ambiguity<0.95f) {
00082       float dx = sift2[k].xpos - sift1[j].xpos;
00083       float dy = sift2[k].ypos - sift1[j].ypos;
00084       int len = (int)(fabs(dx)>fabs(dy) ? fabs(dx) : fabs(dy));
00085       for (int l=0;l<len;l++) {
00086         int x = (int)(sift1[j].xpos + dx*l/len);
00087         int y = (int)(sift1[j].ypos + dy*l/len);
00088         h_img[y*w+x] = 1.0f;
00089       }
00090     }
00091     int p = (int)(sift1[j].ypos+0.5)*w + (int)(sift1[j].xpos+0.5);
00092     p += (w+1);
00093     for (int k=0;k<(int)(1.41*sift1[j].scale);k++)
00094       h_img[p-k] = h_img[p+k] = h_img[p-k*w] =h_img[p+k*w] = 0.0f;
00095     p -= (w+1);
00096     for (int k=0;k<(int)(1.41*sift1[j].scale);k++)
00097       h_img[p-k] = h_img[p+k] = h_img[p-k*w] =h_img[p+k*w] = 1.0f;
00098   }
00099 
00100   FreeSiftData(&siftData1);
00101   FreeSiftData(&siftData2);
00102   limg.Store("data/limg_pts.pgm", true, false);
00103 
00104   img1.h_data = NULL;
00105   FreeCudaImage(&img1);
00106   img2.h_data = NULL;
00107   FreeCudaImage(&img2);
00108   CUT_EXIT(argc, argv);
00109 }
Generated on Sun May 8 08:40:37 2011 for iLab Neuromorphic Vision Toolkit by  doxygen 1.6.3