C++ IpPairVec类代码示例

int mainStaticMatch()
{
  IplImage *img1, *img2;
  img1 = cvLoadImage("imgs/img1.jpg");
  img2 = cvLoadImage("imgs/img2.jpg");

  IpVec ipts1, ipts2;
  surfDetDes(img1,ipts1,false,4,4,2,0.0001f);
  surfDetDes(img2,ipts2,false,4,4,2,0.0001f);

  IpPairVec matches;
  getMatches(ipts1,ipts2,matches);

  for (unsigned int i = 0; i < matches.size(); ++i)
  {
    drawPoint(img1,matches[i].first);
    drawPoint(img2,matches[i].second);
  
    const int & w = img1->width;
    cvLine(img1,cvPoint(matches[i].first.x,matches[i].first.y),cvPoint(matches[i].second.x+w,matches[i].second.y), cvScalar(255,255,255),1);
    cvLine(img2,cvPoint(matches[i].first.x-w,matches[i].first.y),cvPoint(matches[i].second.x,matches[i].second.y), cvScalar(255,255,255),1);
  }

  std::cout<< "Matches: " << matches.size();

  cvNamedWindow("1", CV_WINDOW_AUTOSIZE );
  cvNamedWindow("2", CV_WINDOW_AUTOSIZE );
  cvShowImage("1", img1);
  cvShowImage("2",img2);
  cvWaitKey(0);

  return 0;
}

//! Populate IpPairVec with matched ipts 
void getMatches(IpVec &ipts1, IpVec &ipts2, IpPairVec &matches) {
	float dist, d1, d2;
	Ipoint *match;

	matches.clear();

	for (unsigned int i = 0; i < ipts1.size(); i++) {
		d1 = d2 = 1000;

		for (unsigned int j = 0; j < ipts2.size(); j++) {
			dist = ipts1[i] - ipts2[j];

			if (dist < d1) // if this feature matches better than current best
			{
				d2 = d1;
				d1 = dist;
				match = &ipts2[j];
			} else if (dist < d2) // this feature matches better than second best
			{
				d2 = dist;
			}
		}

		// If match has a d1:d2 ratio < 0.65 ipoints are a match
		if (d1 / d2 < 0.65) {
			// Store the change in position
			ipts1[i].dx = match->x - ipts1[i].x;
			ipts1[i].dy = match->y - ipts1[i].y;
			match->dx =  ipts1[i].x -match->x;
			match->dy =  ipts1[i].y -match->y;

			matches.push_back(std::make_pair(ipts1[i], *match));
		}
	}
}

int ofxSurfObjCorners(IpPairVec & matches,const ofPoint src_crn[4],ofPoint dst_crn[4]) {
    double h[9];
    CvMat _h = cvMat(3,3,CV_64F,h);
    vector<CvPoint2D32f> pt1,pt2;
    CvMat _pt1,_pt2;

    int n = (int)(matches.size());
    if(n<4)return 0;

    pt1.resize(n);
    pt2.resize(n);

    for(int i=0; i<n; i++) {
        pt1[i] = cvPoint2D32f(matches[i].second.x,matches[i].second.y);
        pt2[i] = cvPoint2D32f(matches[i].first.x,matches[i].first.y);
    }
    _pt1 = cvMat(1,n,CV_32F,&pt1[0]);
    _pt2 = cvMat(1,n,CV_32F,&pt2[0]);



    //if(!cvFindHomography(&_pt1,&_pt2,&_h,CV_RANSAC,5))return 0;
    /*for(int i=0;i<4;i++){
        double x = (double)src_crn[i].x;
        double y = (double)src_crn[i].y;
        double Z = 1./(h[6]*x + h[7]*y + h[8]);
        double X = (h[0]*x + h[1]*y + h[2])*Z;
        double Y = (h[3]*x + h[4]*y + h[5])*Z;
        dst_crn[i].set(cvRound(X),cvRound(Y));
    }*/

    return 1;
}

int mainStaticMatch( IplImage *img1, IplImage *img2)
{
//IplImage *img1, *img2;
//img1 = cvLoadImage("../imgs/img1.jpg");
//img2 = cvLoadImage("../imgs/img2.jpg");

  IpVec ipts1, ipts2;

  LARGE_INTEGER llPerfCount = {0};
  QueryPerformanceCounter(&llPerfCount);
  __int64 beginPerfCount = llPerfCount.QuadPart;

  //surfDetDes(img1,ipts1,false,4,4,2,0.0001f);
  //surfDetDes(img2,ipts2,false,4,4,2,0.0001f);
  surfDetDes(img1,ipts1,true,4,4,2,0.0001f);
  surfDetDes(img2,ipts2,true,4,4,2,0.0001f);

  IpPairVec matches;
  getMatches(ipts1,ipts2,matches);

  QueryPerformanceCounter(&llPerfCount);
  __int64 endPerfCount = llPerfCount.QuadPart;
  LARGE_INTEGER liPerfFreq={0};
  QueryPerformanceFrequency(&liPerfFreq);
  std::cout << __FUNCTION__ << " excute time: " 
	  <<  float(endPerfCount - beginPerfCount) * 1000 / liPerfFreq.QuadPart  << " millisecond(ºÁÃë)" << std::endl;

  for (unsigned int i = 0; i < matches.size(); ++i)
  {
    drawPoint(img1,matches[i].first);
    drawPoint(img2,matches[i].second);
  
    const int & w = img1->width;
    cvLine(img1,cvPoint(matches[i].first.x,matches[i].first.y),cvPoint(matches[i].second.x+w,matches[i].second.y), cvScalar(255,255,255),1);
    cvLine(img2,cvPoint(matches[i].first.x-w,matches[i].first.y),cvPoint(matches[i].second.x,matches[i].second.y), cvScalar(255,255,255),1);
  }

  std::cout<< "Matches: " << matches.size();

  cvNamedWindow("1", CV_WINDOW_AUTOSIZE );
  cvNamedWindow("2", CV_WINDOW_AUTOSIZE );
  cvShowImage("1", img1);
  cvShowImage("2",img2);
  cvWaitKey(0);

  return 0;
}

int mainMotionPoints(void)
{
  // Initialise capture device
  CvCapture* capture = cvCaptureFromCAM( CV_CAP_ANY );
  if(!capture) error("No Capture");

  // Create a window 
  cvNamedWindow("OpenSURF", CV_WINDOW_AUTOSIZE );

  // Declare Ipoints and other stuff
  IpVec ipts, old_ipts, motion;
  IpPairVec matches;
  IplImage *img;

  // Main capture loop
  while( 1 ) 
  {
    // Grab frame from the capture source
    img = cvQueryFrame(capture);

    // Detect and describe interest points in the image
    old_ipts = ipts;
    surfDetDes(img, ipts, true, 3, 4, 2, 0.0004f);

    // Fill match vector
    getMatches(ipts,old_ipts,matches);
    for (unsigned int i = 0; i < matches.size(); ++i) 
    {
      const float & dx = matches[i].first.dx;
      const float & dy = matches[i].first.dy;
      float speed = sqrt(dx*dx+dy*dy);
      if (speed > 5 && speed < 30) 
        drawIpoint(img, matches[i].first, 3);
    }
        
    // Display the result
    cvShowImage("OpenSURF", img);

    // If ESC key pressed exit loop
    if( (cvWaitKey(10) & 255) == 27 ) break;
  }

  // Release the capture device
  cvReleaseCapture( &capture );
  cvDestroyWindow( "OpenSURF" );
  return 0;
}

void
getMatches(IpVec &ipts1, IpVec &ipts2, IpPairVec &matches)
{
	// Allocate host memory
	float *lin_descs1 = new float[ipts1.size() * 64];
	float *lin_descs2 = new float[ipts2.size() * 64];
	int *indices = new int[ipts1.size()];
	float *dists = new float[ipts1.size()];

	linearizeDescriptors(lin_descs1, ipts1);
	linearizeDescriptors(lin_descs2, ipts2);

	// Allocate GPU memory
	float *d_lin_descs1;
	float *d_lin_descs2;
	int *d_indices;
	float *d_dists;
	CUDA_SAFE_CALL( cudaMalloc((void **)&d_lin_descs1, ipts1.size() * 64 * sizeof(float)) );
	CUDA_SAFE_CALL( cudaMalloc((void **)&d_lin_descs2, ipts2.size() * 64 * sizeof(float)) );
	CUDA_SAFE_CALL( cudaMalloc((void **)&d_indices, ipts1.size() * sizeof(int)) );
	CUDA_SAFE_CALL( cudaMalloc((void **)&d_dists, ipts1.size() * sizeof(float)) );

	CUDA_SAFE_CALL( cudaMemcpy(d_lin_descs1, lin_descs1, ipts1.size() * 64 * sizeof(float), cudaMemcpyHostToDevice) );
	CUDA_SAFE_CALL( cudaMemcpy(d_lin_descs2, lin_descs2, ipts2.size() * 64 * sizeof(float), cudaMemcpyHostToDevice) );

	prepare_matchSURFKeypointsGPU(0.77f);
	matchSURFKeypointsGPU(d_indices, d_dists,
		d_lin_descs1, ipts1.size(), 64 * sizeof(float),
		d_lin_descs2, ipts2.size(), 64 * sizeof(float),
		64);

	CUDA_SAFE_CALL( cudaMemcpy(indices, d_indices, ipts1.size() * sizeof(int), cudaMemcpyDeviceToHost) );
	CUDA_SAFE_CALL( cudaMemcpy(dists, d_dists, ipts1.size() * sizeof(float), cudaMemcpyDeviceToHost) );

	for (size_t i = 0; i < ipts1.size(); i++)
	{
		if (indices[i] != -1)
		{
			Ipoint &ipt1 = ipts1[i];
			const Ipoint &match = ipts2[indices[i]];

			ipt1.dx = match.x - ipt1.x;
			ipt1.dy = match.y - ipt1.y;
			matches.push_back(std::make_pair(ipt1, match));
		}
	}

	delete[] lin_descs1;
	delete[] lin_descs2;
	delete[] indices;
	delete[] dists;

	CUDA_SAFE_CALL( cudaFree(d_lin_descs1) );
	CUDA_SAFE_CALL( cudaFree(d_lin_descs2) );
	CUDA_SAFE_CALL( cudaFree(d_indices) );
	CUDA_SAFE_CALL( cudaFree(d_dists) );
}

int mainStaticMatch()
{
  // Make images as Mats; convert to IplImage for OpenSURF library actions
  cv::Mat mimg1, mimg2;
  mimg1=cv::imread("OpenSURF/imgs/img1.jpg", CV_LOAD_IMAGE_COLOR);
  mimg2=cv::imread("OpenSURF/imgs/img2.jpg", CV_LOAD_IMAGE_COLOR);

  IplImage iimg1, iimg2;
  iimg1=mimg1;
  iimg2=mimg2;

  IplImage *img1, *img2;
  img1 = &iimg1;
  img2 = &iimg2;

  IpVec ipts1, ipts2;
  surfDetDes(img1,ipts1,false,4,4,2,0.0001f);
  surfDetDes(img2,ipts2,false,4,4,2,0.0001f);

  IpPairVec matches;
  getMatches(ipts1,ipts2,matches);

  for (unsigned int i = 0; i < matches.size(); ++i)
  {
    drawPoint(img1,matches[i].first);
    drawPoint(img2,matches[i].second);
  
    const int & w = img1->width;
    cvLine(img1,cvPoint(matches[i].first.x,matches[i].first.y),cvPoint(matches[i].second.x+w,matches[i].second.y), cvScalar(255,255,255),1);
    cvLine(img2,cvPoint(matches[i].first.x-w,matches[i].first.y),cvPoint(matches[i].second.x,matches[i].second.y), cvScalar(255,255,255),1);
  }

  std::cout<< "Matches: " << matches.size();

  cvNamedWindow("1", CV_WINDOW_AUTOSIZE );
  cvNamedWindow("2", CV_WINDOW_AUTOSIZE );
  cvShowImage("1", img1);
  cvShowImage("2",img2);
  cvWaitKey(0);

  return 0;
}

//! Find homography between matched points and translate src_corners to dst_corners
int translateCorners(IpPairVec &matches, const CvPoint src_corners[4], CvPoint dst_corners[4])
{
#ifndef LINUX
  double h[9];
  cv::Mat _h = cv::Mat(3, 3, CV_64F, h);
  std::vector<CvPoint2D32f> pt1, pt2;
  cv::Mat _pt1, _pt2;
  
  int n = (int)matches.size();
  if( n < 4 ) return 0;

  // Set vectors to correct size
  pt1.resize(n);
  pt2.resize(n);

  // Copy Ipoints from match vector into cvPoint vectors
  for(int i = 0; i < n; i++ )
  {
    pt1[i] = cvPoint2D32f(matches[i].second.x, matches[i].second.y);
    pt2[i] = cvPoint2D32f(matches[i].first.x, matches[i].first.y);
  }
  _pt1 = cv::Mat(1, n, CV_32FC2, &pt1[0] );
  _pt2 = cv::Mat(1, n, CV_32FC2, &pt2[0] );

  cv::InputArray _pt1_ia(_pt1);
  cv::InputArray _pt2_ia(_pt2);
  cv::OutputArray _h_ia(_h);

  cv::Mat hMat = cv::findHomography(_pt1_ia, _pt2_ia, _h_ia, CV_RANSAC, 5.0);
   // Find the homography (transformation) between the two sets of points
  if(hMat.empty())  // this line requires opencv 1.1
    return 0;

  // Translate src_corners to dst_corners using homography
  for(int i = 0; i < 4; i++ )
  {
    double x = src_corners[i].x, y = src_corners[i].y;
    double Z = 1./(h[6]*x + h[7]*y + h[8]);
    double X = (h[0]*x + h[1]*y + h[2])*Z;
    double Y = (h[3]*x + h[4]*y + h[5])*Z;
    dst_corners[i] = cvPoint(cvRound(X), cvRound(Y));
  }
#endif
  return 1;
}

int mainMatch(void)
{
  // Initialise capture device
  CvCapture* capture = cvCaptureFromCAM( CV_CAP_ANY );
  if(!capture) error("No Capture");

  // Declare Ipoints and other stuff
  IpPairVec matches;
  IpVec ipts, ref_ipts;
  
  // This is the reference object we wish to find in video frame
  // Replace the line below with IplImage *img = cvLoadImage("imgs/object.jpg"); 
  // where object.jpg is the planar object to be located in the video
  IplImage *img = cvLoadImage("../imgs/object.jpg"); 
  if (img == NULL) error("Need to load reference image in order to run matching procedure");
  CvPoint src_corners[4] = {{0,0}, {img->width,0}, {img->width, img->height}, {0, img->height}};
  CvPoint dst_corners[4];

  // Extract reference object Ipoints
  surfDetDes(img, ref_ipts, false, 3, 4, 3, 0.004f);
  drawIpoints(img, ref_ipts);
  showImage(img);

  // Create a window 
  cvNamedWindow("OpenSURF", CV_WINDOW_AUTOSIZE );

  // Main capture loop
  while( true ) 
  {
    // Grab frame from the capture source
    img = cvQueryFrame(capture);
     
    // Detect and describe interest points in the frame
    surfDetDes(img, ipts, false, 3, 4, 3, 0.004f);

    // Fill match vector
    getMatches(ipts,ref_ipts,matches);
    
    // This call finds where the object corners should be in the frame
    if (translateCorners(matches, src_corners, dst_corners))
    {
      // Draw box around object
      for(int i = 0; i < 4; i++ )
      {
        CvPoint r1 = dst_corners[i%4];
        CvPoint r2 = dst_corners[(i+1)%4];
        cvLine( img, cvPoint(r1.x, r1.y),
          cvPoint(r2.x, r2.y), cvScalar(255,255,255), 3 );
      }

      for (unsigned int i = 0; i < matches.size(); ++i)
        drawIpoint(img, matches[i].first);
    }

    // Draw the FPS figure
    drawFPS(img);

    // Display the result
    cvShowImage("OpenSURF", img);

    // If ESC key pressed exit loop
    if( (cvWaitKey(10) & 255) == 27 ) break;
  }

  // Release the capture device
  cvReleaseCapture( &capture );
  cvDestroyWindow( "OpenSURF" );
  return 0;
}

float visualStaticMatch(TestImage image1, TestImage image2, bool visual, IplImage* &half)
{
	// Get matches
	IpPairVec matches;
	//	double time = (double)cvGetTickCount();
	Score result = getMatchesRANSAC(image1.ipts,image2.ipts,matches);
	float	match_score = result.score;
	//	double ttmatch = (double)cvGetTickCount() - time;

	if(visual){

		IplImage* img1 = image1.img;
		IplImage* img2 = image2.img;
		// Display the horizon lines on the images
		int integral_height = WIDTH;

		//Draws the horizon lines
		cvLine(img1, cvPoint(0, image1.left+WIDTH/2),cvPoint(img1->width, image1.right+WIDTH/2), cvScalar(0, 0, 255),1);
		cvLine(img1, cvPoint(0, image1.left-WIDTH/2),cvPoint(img1->width, image1.right-WIDTH/2), cvScalar(0, 0, 255),1);
		cvLine(img2, cvPoint(0, image2.left+WIDTH/2),cvPoint(img2->width, image2.right+WIDTH/2), cvScalar(0, 0, 255),1);
		cvLine(img2, cvPoint(0, image2.left-WIDTH/2),cvPoint(img2->width, image2.right-WIDTH/2), cvScalar(0, 0, 255),1);

		// Copy the original images to the output
		IplImage* visual = cvCreateImage( cvSize(img1->width + img2->height + integral_height,
				img2->width + img1->height + integral_height), img1->depth, img1->nChannels );
		cvZero(visual);
		cvSetImageROI( visual, cvRect( img2->height + integral_height, img2->width+ integral_height, img1->width, img1->height ) );
		cvCopy( img1, visual );
		cvSetImageROI( visual, cvRect( 0, 0, img2->height, img2->width ) );

		cvFlip(img2, NULL, 0);
		//cvFlip(img2, NULL, 1);
		cvTranspose(img2, visual);
		//cvFlip(img2, NULL, 1);
		cvFlip(img2, NULL, 0);

		// Display the integral image above each image
		std::vector<int> y_coords;
		IplImage *int_img1 = getGrayHorizon(img1, image1.left, image1.right, y_coords, WIDTH, 1, true);
		IplImage *im8 = cvCreateImage(cvSize(int_img1->width, int_img1->height), visual->depth, int_img1->nChannels);
		cvConvertScale(int_img1, im8, 255.);
		IplImage *colint_img1 = cvCreateImage( cvSize(int_img1->width, int_img1->height), visual->depth, visual->nChannels );
		cvCvtColor( im8, colint_img1, CV_GRAY2BGR );
		cvSetImageROI( visual, cvRect( img2->height + integral_height, img2->width, int_img1->width, integral_height ) );
		cvResize(colint_img1, visual);
		cvReleaseImage(&int_img1);
		cvReleaseImage(&im8);
		cvReleaseImage(&colint_img1);

		IplImage *int_img2 = getGrayHorizon(img2, image2.left, image2.right, y_coords, WIDTH, 1, true);
		im8 = cvCreateImage(cvSize(int_img2->width, int_img2->height), visual->depth, int_img2->nChannels);
		cvConvertScale(int_img2, im8, 255.);
		IplImage *colint_img2 = cvCreateImage( cvSize(int_img2->width, int_img2->height), visual->depth, visual->nChannels );
		cvCvtColor( im8, colint_img2, CV_GRAY2BGR );
		cvSetImageROI( visual, cvRect( img2->height , 0, integral_height, int_img2->width  ) );

		IplImage *transpose = cvCreateImage( cvSize(colint_img2->height, colint_img2->width), colint_img2->depth, colint_img2->nChannels );
		cvFlip(colint_img2, NULL, 0);
		//cvFlip(colint_img2, NULL, 1);
		cvTranspose(colint_img2, transpose);
		cvResize(transpose, visual);

		cvReleaseImage(&transpose);
		cvReleaseImage(&int_img2);
		cvReleaseImage(&im8);
		cvReleaseImage(&colint_img2);

		// Draw match points
		cvSetImageROI( visual, cvRect( img2->height+integral_height, 0, img1->width, img2->width ) );
		cvSet(visual,cvScalar(255, 255, 255));
		for (unsigned int i = 0; i < matches.size(); ++i){
			drawMatch(visual, matches[i], img2->width, visual);
		}

		// Display text
		cvSetImageROI( visual, cvRect( 0 , img2->width, img2->height + integral_height, img2->height + integral_height  ) );
		char text[200];
		CvFont font;
		cvInitFont(&font,CV_FONT_HERSHEY_SIMPLEX, 0.9,0.9,0,2);
		sprintf(text,"Features: %d/%d",(int)image2.ipts.size(),(int)image1.ipts.size());
		cvPutText (visual, text, cvPoint(20,50), &font, cvScalar(255,255,255));
		sprintf(text,"Valid matches: %d",(int)matches.size());
		cvPutText (visual, text, cvPoint(20,100), &font, cvScalar(255,255,255));
		sprintf(text,"Recognition score: %.2f",match_score);
		cvPutText (visual, text, cvPoint(20,150), &font, cvScalar(255,255,255));
		sprintf(text,"Extraction time: %.2fms",Ipoint::extractionTime);//image2.tt_extraction, image1.tt_extraction);
		cvPutText (visual, text, cvPoint(20,200), &font, cvScalar(255,255,255));
		sprintf(text,"Matching time: %.2fms",Ipoint::detectionTime);
		cvPutText (visual, text, cvPoint(20,250), &font, cvScalar(255,255,255));

		// Resize and Show output
		cvResetImageROI( visual );
		half = cvCreateImage( cvSize(visual->width/2, visual->height/2), visual->depth, visual->nChannels );
		cvResize(visual, half);
		cvReleaseImage(&visual);

	}
	return match_score;
}

void msrmtupdate(unsigned char *img,double updt[][9],IpVec *refpts,double
*refhist,int w,int h)
{
    long double tot=0,dist2=0;
    double hist[257 * 10]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
    static double ptr[4],p;
    int key[10][2];
    IpVec pts,tpts;
    IplImage *cl,*tmp2;
    cl=cvCreateImage(cvSize(w,h),IPL_DEPTH_8U,3);
    memcpy(cl->imageData,(img),cl->imageSize);
    tmp2=cvCloneImage(cl);
    surfDetDes(cl,tpts ,false, 5, 4, 2, 0.00004f);
    IpPairVec matches;
    IpVec ipts, ref_ipts;
    CvPoint src_corners[4] = {cvPoint(0,0), cvPoint(80,0), cvPoint(80, 60),
    cvPoint(0, 60)};
    CvPoint dst_corners[4];
    getMatches(tpts,*refpts,matches);
    int tt=0;
    tt=translateCorners(matches, src_corners, dst_corners);
    if (translateCorners(matches, src_corners, dst_corners))
     {
         // Draw box around object
         for(int i = 0; i < 4; i++ )

         {
             CvPoint r1 = dst_corners[i%4];
             CvPoint r2 = dst_corners[(i+1)%4];
             cvLine( cl, cvPoint(r1.x, r1.y), cvPoint(r2.x, r2.y), cvScalar(255,255,255),3 );
         }
         for (unsigned int i = 0; i < matches.size(); ++i)
         drawIpoint(cl, matches[i].first);
     }
    CvPoint cpt;
    cpt.x=((dst_corners[0].x)+(dst_corners[2].x))/2;
    cpt.y=((dst_corners[0].y)+(dst_corners[2].y))/2;
    p++;
    if(tt)
    {
        if((abs(ptr[2]-abs(dst_corners[0].x-dst_corners[1].x))>=30 ||
        abs(ptr[3]-abs(dst_corners[0].y-dst_corners[3].y))>=30 ||
        !isrect(dst_corners)) && p>3 )
        {
             tt=0;
        }
        else
        {
            cvCvtColor(tmp2,cl ,CV_RGB2HSV);
            ptr[0]=cpt.x;ptr[1]=cpt.y;ptr[2]=abs(dst_corners[0].xst_corners[1].x);ptr[3]=abs(dst_corners[0].y-dst_corners[3].y);
            crhist((unsigned char *)cl->imageData,hist,w,h,ptr);
            dist2=.1*(double)exp(-2*pow(comphist(hist,refhist),2));
        }
    }
    for(int i=0;i<N;i++)
    {
    if(tt && dist2>.05 )
    {
        updt[i][0]=cpt.x;
        updt[i][1]=cpt.y;
        updt[i][2]=ptr[2];
        updt[i][3]=ptr[3];
        updt[i][4]=1;
        updt[i][5]=1;
        updt[i][8]=1;
        tot++;
    }
    else
     {
        double pt[4];
        for(int k=0;k<4;k++)
        {
            pt[k]=updt[i][k];
        }
        cvCvtColor(tmp2,cl, CV_RGB2HSV);
        crhist((unsigned char *)cl->imageData,hist,w,h,pt);
        dist2=.1*(double)exp(-100*pow(comphist(hist,refhist),2));
        updt[i][8]=dist2;
        tot+=updt[i][8];
    }
    }
    for(int i=0;i<N;i++)
      updt[i][8]/=(double)tot;

}

//! Populate IpPairVec with matched ipts using nearest neighbour and RANSAC
Score getMatchesRANSAC(IpVec &ipts1, IpVec &ipts2, IpPairVec &matches)
{
#if RUNSWIFT
#else
	timespec matchings, matchinge, verifys, verifye;

	clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &matchings);
#endif
	float dist, d1, d2;
	Ipoint *match;
	float matching_score = 0;

	matches.clear();

	for(unsigned int i = 0; i < ipts1.size(); i++)
	{
		ipts1[i].x = ipts1[i].x;
		d1 = d2 = FLT_MAX;
		match = &ipts2[0]; // to avoid unitialized warning

		for(unsigned int j = 0; j < ipts2.size(); j++)
		{
			ipts2[j].x = ipts2[j].x;
			dist = ipts1[i] - ipts2[j];

			if(dist<d1) // if this feature matches better than current best
			{
				d2 = d1;
				d1 = dist;
				match = &ipts2[j];
			}
			else if(dist<d2) // this feature matches better than second best
			{
				d2 = dist;
			}
		}

		// If match has a d1:d2 ratio < 0.75 ipoints are a match
		if(d1/d2 < 0.75)
		{
			// Store the match
			matches.push_back(std::make_pair(ipts1[i], *match));
			//Increment the matching score
			matching_score += 1/d1;
		}
	}

	float best_score = matching_score;
	float best_b = -1;
	float best_m = -1;
#if RUNSWIFT
#else
	Ipoint::totalNumMatches = matches.size();
	//At this point we have the total matches before the final number of matches
	clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &verifys);
#endif

	if(matches.size()>1){

		best_score = 0;

		for(int i=0; i<ITERATIONS; i++){
			//Choose random matches
			int pos1 = rand() % (int)matches.size();
			int pos2 = rand() % (int)matches.size();
			while(pos1 == pos2) {
				//Make sure that both matches are different
				pos2 = rand() % (int)matches.size();
			}
			//Should generate a positive value
			float m = (matches.at(pos2).second.x - matches.at(pos1).second.x)/(matches.at(pos2).first.x - matches.at(pos1).first.x);
			//If a gradient is discarded
			if (m <= 0){
				continue;
			}
			//Calculate the translation component
			float b = matches.at(pos2).second.x - m*matches.at(pos2).first.x;
			float score = 0;
			for(int j=0; j<(int)matches.size(); j++){
				//Calculate the function x_stored,i = b_s * x_test,i + b_d
				if( fabs(matches.at(j).second.x - (m*matches.at(j).first.x + b)) < PIXEL_ERROR_MARGIN)
					score += 1/fabs(matches.at(j).first - matches.at(j).second);
			}
			if (score > best_score){
				best_score = score;
				best_b = b;
				best_m = m;
			}	
		}
	}

	// Now remove all matches who are not within this pixel error margin
	//if(best_m > 0){
	for(int j=0; j<(int)matches.size(); j++){
		if( fabs(matches.at(j).second.x - (best_m*matches.at(j).first.x + best_b)) >= PIXEL_ERROR_MARGIN) {
			matches.erase(matches.begin() + j);
			j--;
		}
	}
//.........这里部分代码省略.........

int mainStaticMatch()
{

    time_t start,end1,end2,end3,end4,end5;
    start = clock();

    IplImage *img1, *img2;
    img1 = cvLoadImage("../data/1.JPG");
    img2 = cvLoadImage("../data/2.JPG");


    end1 = clock();

    IpVec ipts1, ipts2;
    surfDetDes(img1,ipts1,false,4,4,2,0.0008f);
    surfDetDes(img2,ipts2,false,4,4,2,0.0008f);

    std::cout << "im1" << std::endl;
    std::cout << "Size:" << ipts1.size() << std::endl;

    std::cout << "im2" << std::endl;
    std::cout << "Size:" << ipts2.size() << std::endl;
    end2 = clock();

    IpPairVec matches;
    getMatches(ipts1,ipts2,matches);

    end3 = clock();

    for (unsigned int i = 0; i < matches.size(); ++i)
    {
        drawPoint(img1,matches[i].first);
        drawPoint(img2,matches[i].second);

        const int & w = img1->width;
        cvLine(img1,cvPoint(matches[i].first.x,matches[i].first.y),cvPoint(matches[i].second.x+w,matches[i].second.y), cvScalar(255,255,255),1);
        cvLine(img2,cvPoint(matches[i].first.x-w,matches[i].first.y),cvPoint(matches[i].second.x,matches[i].second.y), cvScalar(255,255,255),1);
    }

    std::cout << "Matches: " << matches.size() << std::endl;
    /*
      cvNamedWindow("1", CV_WINDOW_AUTOSIZE );
      cvNamedWindow("2", CV_WINDOW_AUTOSIZE );
      cvShowImage("1", img1);
      cvShowImage("2", img2);
      cvWaitKey(0);
    */
    end4 = clock();

//  cvSaveImage("result_gpu1.jpg",img1);
//	cvSaveImage("result_gpu2.jpg",img2);

    // Stitch two images
    IplImage *img = cvCreateImage(cvSize(img1->width + img2->width,
                                         img1->height),img1->depth,img1->nChannels);
    cvSetImageROI( img, cvRect( 0, 0, img1->width, img1->height ) );
    cvCopy(img1, img);
    cvSetImageROI( img, cvRect(img1->width,0, img2->width, img2->height) );
    cvCopy(img2, img);
    cvResetImageROI(img);
    cvSaveImage("result_gpu.jpg",img);

    end5 = clock();
    double dif1 = (double)(end1 - start) / CLOCKS_PER_SEC;
    double dif2 = (double)(end2 - end1) / CLOCKS_PER_SEC;
    double dif3 = (double)(end3 - end2) / CLOCKS_PER_SEC;
    double dif4 = (double)(end4 - end3) / CLOCKS_PER_SEC;
    double dif5 = (double)(end5 - end4) / CLOCKS_PER_SEC;
    double total = (double)(end5 - start) / CLOCKS_PER_SEC;
    std::cout.setf(std::ios::fixed,std::ios::floatfield);
    std::cout.precision(5);
    std::cout << "Time(load):" << dif1 << std::endl;
    std::cout << "Time(descriptor):" << dif2 << std::endl;
    std::cout << "Time(match):" << dif3 << std::endl;
    std::cout << "Time(plot):" << dif4 << std::endl;
    std::cout << "Time(save):" << dif5 << std::endl;
    std::cout << "Time(Total):" << total << std::endl;
    return 0;
}