contrib/oxl/mvl/ComputeGRIC.cxx
Go to the documentation of this file.
00001 #include "ComputeGRIC.h"
00002 #include <mvl/FMatrixComputeMLESAC.h>
00003 #include <mvl/HMatrix2DComputeMLESAC.h>
00004 #include <mvl/HomgInterestPointSet.h>
00005 #include <vcl_cmath.h>
00006 #include <vcl_iostream.h>
00007 #include <vcl_memory.h>
00008 
00009 ComputeGRIC::ComputeGRIC(double std) : std_(std) {}
00010 
00011 ComputeGRIC::~ComputeGRIC(){}
00012 
00013 bool ComputeGRIC::compute(PairMatchSetCorner* matches) {
00014 
00015   // Make some copies of the corner matches so they are
00016   // not violated
00017   PairMatchSetCorner matches_copy1(*matches);
00018   PairMatchSetCorner matches_copy2(*matches);
00019 
00020   // Compute the two differing relations using MLE
00021   F_.reset(new FMatrix);
00022   vcl_auto_ptr<FMatrixComputeRobust> computor1(new FMatrixComputeMLESAC(true, std_));
00023   computor1->compute(matches_copy1, F_.get());
00024   residualsF_ = computor1->get_residuals();
00025   inliersF_ = computor1->get_inliers();
00026   basisF_ = computor1->get_basis();
00027   computor1.reset();
00028 
00029   H_.reset(new HMatrix2D);
00030   vcl_auto_ptr<HMatrix2DComputeRobust> computor2(new HMatrix2DComputeMLESAC(std_));
00031   computor2->compute(matches_copy2, H_.get());
00032   residualsH_ = computor2->get_residuals();
00033   inliersH_ = computor2->get_inliers();
00034   basisH_ = computor2->get_basis();
00035   computor2.reset();
00036 
00037   vcl_cerr << "Finished calculations" << vcl_endl;
00038 
00039   // Compare the GRIC scores of the two different models
00040   int inf = 0, inh = 0;
00041   double stdf = 0.0, stdh = 0.0;
00042   int n = 0;
00043   for (unsigned int i = 0; i < inliersF_.size(); i++) {
00044     n++;
00045     if (inliersF_[i]) {
00046       inf++;
00047       stdf += residualsF_[i];
00048     }
00049   }
00050   for (unsigned int i = 0; i < inliersH_.size(); i++) {
00051     if (inliersH_[i]) {
00052       inh++;
00053       stdh += residualsH_[i];
00054     }
00055   }
00056   stdf /= inf;
00057   stdh /= inh;
00058   vcl_cerr << "inf : " << inf << vcl_endl;
00059   vcl_cerr << "inh : " << inh << vcl_endl;
00060   vcl_cerr << "stdf : " << stdf << vcl_endl;
00061   vcl_cerr << "stdh : " << stdh << vcl_endl;
00062   int df = 3, dh = 2, r = 4, kf = 7, kh = 8;
00063   double l1 = vcl_log(4.0), l2 = vcl_log(4.0*n), l3 = 2.0;
00064   double GRICF = 0.0;
00065   double thresh = l3*(r - df);
00066   for (unsigned int i = 0; i < residualsF_.size(); i++) {
00067     double t = residualsF_[i] / std_;
00068     if (t < thresh)
00069       GRICF += t;
00070     else
00071       GRICF += thresh;
00072   }
00073   GRICF += l1*(df*n) + l2*kf;
00074   vcl_cerr << "GRICF : " << GRICF << vcl_endl;
00075   double GRICH = 0.0;
00076   thresh = l3*(r - dh);
00077   for (unsigned int i = 0; i < residualsH_.size(); i++) {
00078     double t = residualsH_[i] / std_;
00079     if (t < thresh)
00080       GRICH += t;
00081     else
00082       GRICH += thresh;
00083   }
00084   GRICH += l1*(dh*n) + l2*kh;
00085   vcl_cerr << "GRICH : " << GRICH << vcl_endl;
00086 
00087   // Determine the winner
00088   if (GRICH < GRICF) {
00089     degenerate_ = true;
00090     return false;
00091   } else {
00092     degenerate_ = false;
00093     return true;
00094   }
00095 }