xform.c
上传用户:lwxipeng
上传日期:2022-05-16
资源大小:15982k
文件大小:16k
源码类别:

视频捕捉/采集

开发平台:

Visual C++

  1. /*
  2. This file contains definitions for functions to compute transforms from
  3. image feature correspondences
  4. Copyright (C) 2006  Rob Hess <hess@eecs.oregonstate.edu>
  5. @version 1.1.1-20070330
  6. */
  7. #include "xform.h"
  8. #include "imgfeatures.h"
  9. #include "utils.h"
  10. #include <cxcore.h>
  11. #include <gslgsl_sf.h>
  12. #include <gslgsl_rng.h>
  13. #include <gslgsl_randist.h>
  14. #include <time.h>
  15. /************************* Local Function Prototypes *************************/
  16. static __inline struct feature* get_match( struct feature*, int );
  17. int get_matched_features( struct feature*, int, int, struct feature*** );
  18. int calc_min_inliers( int, int, double, double );
  19. struct feature** draw_ransac_sample( struct feature**, int, int, gsl_rng* );
  20. void extract_corresp_pts( struct feature**, int, int, CvPoint2D64f**,
  21.  CvPoint2D64f** );
  22. int find_consensus( struct feature**, int, int, CvMat*, ransac_err_fn,
  23.    double, struct feature*** );
  24. static __inline void release_mem( CvPoint2D64f*, CvPoint2D64f*,
  25. struct feature** );
  26. /********************** Functions prototyped in xform.h **********************/
  27. /*
  28. Calculates a best-fit image transform from image feature correspondences
  29. using RANSAC.
  30. For more information refer to:
  31. Fischler, M. A. and Bolles, R. C.  Random sample consensus: a paradigm for
  32. model fitting with applications to image analysis and automated cartography.
  33. <EM>Communications of the ACM, 24</EM>, 6 (1981), pp. 381--395.
  34. @param features an array of features; only features with a non-NULL match
  35. of type mtype are used in homography computation
  36. @param n number of features in feat
  37. @param mtype determines which of each feature's match fields to use
  38. for model computation; should be one of FEATURE_FWD_MATCH,
  39. FEATURE_BCK_MATCH, or FEATURE_MDL_MATCH; if this is FEATURE_MDL_MATCH,
  40. correspondences are assumed to be between a feature's img_pt field
  41. and its match's mdl_pt field, otherwise correspondences are assumed to
  42. be between the the feature's img_pt field and its match's img_pt field
  43. @param xform_fn pointer to the function used to compute the desired
  44. transformation from feature correspondences
  45. @param m minimum number of correspondences necessary to instantiate the
  46. model computed by xform_fn
  47. @param p_badxform desired probability that the final transformation
  48. returned by RANSAC is corrupted by outliers (i.e. the probability that
  49. no samples of all inliers were drawn)
  50. @param err_fn pointer to the function used to compute a measure of error
  51. between putative correspondences and a computed model
  52. @param err_tol correspondences within this distance of a computed model are
  53. considered as inliers
  54. @param inliers if not NULL, output as an array of pointers to the final
  55. set of inliers
  56. @param n_in if not NULL and a inliers is not NULL, output as the final
  57. number of inliers
  58. @return Returns a transformation matrix computed using RANSAC or NULL
  59. on error or if an acceptable transform could not be computed.
  60. */
  61. CvMat* ransac_xform( struct feature* features, int n, int mtype,
  62. ransac_xform_fn xform_fn, int m, double p_badxform,
  63. ransac_err_fn err_fn, double err_tol,
  64. struct feature*** inliers, int* n_in )
  65. {
  66. struct feature** matched, ** sample, ** consensus, ** consensus_max = NULL;
  67. struct ransac_data* rdata;
  68. CvPoint2D64f* pts, * mpts;
  69. CvMat* M = NULL;
  70. gsl_rng* rng;
  71. double p, in_frac = RANSAC_INLIER_FRAC_EST;
  72. int i, nm, in, in_min, in_max = 0, k = 0;
  73. nm = get_matched_features( features, n, mtype, &matched );
  74. if( nm < m )
  75. {
  76. fprintf( stderr, "Warning: not enough matches to compute xform, %s" 
  77. " line %dn", __FILE__, __LINE__ );
  78. goto end;
  79. }
  80. /* initialize random number generator */
  81. rng = gsl_rng_alloc( gsl_rng_mt19937 );
  82. gsl_rng_set( rng, time(NULL) );
  83. in_min = calc_min_inliers( nm, m, RANSAC_PROB_BAD_SUPP, p_badxform );
  84. p = pow( 1.0 - pow( in_frac, m ), k );
  85. i = 0;
  86. while( p > p_badxform )
  87. {
  88. sample = draw_ransac_sample( matched, nm, m, rng );
  89. extract_corresp_pts( sample, m, mtype, &pts, &mpts );
  90. M = xform_fn( pts, mpts, m );
  91. if( ! M )
  92. goto iteration_end;
  93. in = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus);
  94. if( in > in_max )
  95. {
  96. if( consensus_max )
  97. free( consensus_max );
  98. consensus_max = consensus;
  99. in_max = in;
  100. in_frac = (double)in_max / nm;
  101. }
  102. else
  103. free( consensus );
  104. cvReleaseMat( &M );
  105. iteration_end:
  106. release_mem( pts, mpts, sample );
  107. p = pow( 1.0 - pow( in_frac, m ), ++k );
  108. }
  109. /* calculate final transform based on best consensus set */
  110. if( in_max >= in_min )
  111. {
  112. extract_corresp_pts( consensus_max, in_max, mtype, &pts, &mpts );
  113. M = xform_fn( pts, mpts, in_max );
  114. in = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus);
  115. cvReleaseMat( &M );
  116. release_mem( pts, mpts, consensus_max );
  117. extract_corresp_pts( consensus, in, mtype, &pts, &mpts );
  118. M = xform_fn( pts, mpts, in );
  119. if( inliers )
  120. {
  121. *inliers = consensus;
  122. consensus = NULL;
  123. }
  124. if( n_in )
  125. *n_in = in;
  126. release_mem( pts, mpts, consensus );
  127. }
  128. else if( consensus_max )
  129. {
  130. if( inliers )
  131. *inliers = NULL;
  132. if( n_in )
  133. *n_in = 0;
  134. free( consensus_max );
  135. }
  136. gsl_rng_free( rng );
  137. end:
  138. for( i = 0; i < nm; i++ )
  139. {
  140. rdata = feat_ransac_data( matched[i] );
  141. matched[i]->feature_data = rdata->orig_feat_data;
  142. free( rdata );
  143. }
  144. free( matched );
  145. return M;
  146. }
  147. /*
  148. Calculates a least-squares planar homography from point correspondeces.
  149. @param pts array of points
  150. @param mpts array of corresponding points; each pts[i], i=0..n-1, corresponds
  151. to mpts[i]
  152. @param n number of points in both pts and mpts; must be at least 4
  153. @return Returns the 3 x 3 least-squares planar homography matrix that
  154. transforms points in pts to their corresponding points in mpts or NULL if
  155. fewer than 4 correspondences were provided
  156. */
  157. CvMat* lsq_homog( CvPoint2D64f* pts, CvPoint2D64f* mpts, int n )
  158. {
  159. CvMat* H, * A, * B, X;
  160. double x[9];
  161. int i;
  162. if( n < 4 )
  163. {
  164. fprintf( stderr, "Warning: too few points in lsq_homog(), %s line %dn",
  165. __FILE__, __LINE__ );
  166. return NULL;
  167. }
  168. /* set up matrices so we can unstack homography into X; AX = B */
  169. A = cvCreateMat( 2*n, 8, CV_64FC1 );
  170. B = cvCreateMat( 2*n, 1, CV_64FC1 );
  171. X = cvMat( 8, 1, CV_64FC1, x );
  172. H = cvCreateMat(3, 3, CV_64FC1);
  173. cvZero( A );
  174. for( i = 0; i < n; i++ )
  175. {
  176. cvmSet( A, i, 0, pts[i].x );
  177. cvmSet( A, i+n, 3, pts[i].x );
  178. cvmSet( A, i, 1, pts[i].y );
  179. cvmSet( A, i+n, 4, pts[i].y );
  180. cvmSet( A, i, 2, 1.0 );
  181. cvmSet( A, i+n, 5, 1.0 );
  182. cvmSet( A, i, 6, -pts[i].x * mpts[i].x );
  183. cvmSet( A, i, 7, -pts[i].y * mpts[i].x );
  184. cvmSet( A, i+n, 6, -pts[i].x * mpts[i].y );
  185. cvmSet( A, i+n, 7, -pts[i].y * mpts[i].y );
  186. cvmSet( B, i, 0, mpts[i].x );
  187. cvmSet( B, i+n, 0, mpts[i].y );
  188. }
  189. cvSolve( A, B, &X, CV_SVD );
  190. x[8] = 1.0;
  191. X = cvMat( 3, 3, CV_64FC1, x );
  192. cvConvert( &X, H );
  193. cvReleaseMat( &A );
  194. cvReleaseMat( &B );
  195. return H;
  196. }
  197. /*
  198. Calculates the transfer error between a point and its correspondence for
  199. a given homography, i.e. for a point x, it's correspondence x', and
  200. homography H, computes d(x', Hx)^2.
  201. @param pt a point
  202. @param mpt pt's correspondence
  203. @param H a homography matrix
  204. @return Returns the transfer error between pt and mpt given H
  205. */
  206. double homog_xfer_err( CvPoint2D64f pt, CvPoint2D64f mpt, CvMat* H )
  207. {
  208. CvPoint2D64f xpt = persp_xform_pt( pt, H );
  209. return sqrt( dist_sq_2D( xpt, mpt ) );
  210. }
  211. /*
  212. Performs a perspective transformation on a single point.  That is, for a
  213. point (x, y) and a 3 x 3 matrix T this function returns the point
  214. (u, v), where
  215. [x' y' w']^T = T * [x y 1]^T,
  216. and
  217. (u, v) = (x'/w', y'/w').
  218. Note that affine transforms are a subset of perspective transforms.
  219. @param pt a 2D point
  220. @param T a perspective transformation matrix
  221. @return Returns the point (u, v) as above.
  222. */
  223. CvPoint2D64f persp_xform_pt( CvPoint2D64f pt, CvMat* T )
  224. {
  225. CvMat XY, UV;
  226. double xy[3] = { pt.x, pt.y, 1.0 }, uv[3] = { 0 };
  227. CvPoint2D64f rslt;
  228. cvInitMatHeader( &XY, 3, 1, CV_64FC1, xy, CV_AUTOSTEP );
  229. cvInitMatHeader( &UV, 3, 1, CV_64FC1, uv, CV_AUTOSTEP );
  230. cvMatMul( T, &XY, &UV );
  231. rslt = cvPoint2D64f( uv[0] / uv[2], uv[1] / uv[2] );
  232. return rslt;
  233. }
  234. /************************ Local funciton definitions *************************/
  235. /*
  236. Returns a feature's match according to a specified match type
  237. @param feat feature
  238. @param mtype match type, one of FEATURE_FWD_MATCH, FEATURE_BCK_MATCH, or
  239. FEATURE_MDL_MATCH
  240. @return Returns feat's match corresponding to mtype or NULL for bad mtype
  241. */
  242. static __inline struct feature* get_match( struct feature* feat, int mtype )
  243. {
  244. if( mtype == FEATURE_MDL_MATCH )
  245. return feat->mdl_match;
  246. if( mtype == FEATURE_BCK_MATCH )
  247. return feat->bck_match;
  248. if( mtype == FEATURE_FWD_MATCH )
  249. return feat->fwd_match;
  250. return NULL;
  251. }
  252. /*
  253. Finds all features with a match of a specified type and stores pointers
  254. to them in an array.  Additionally initializes each matched feature's
  255. feature_data field with a ransac_data structure.
  256. @param features array of features
  257. @param n number of features in features
  258. @param mtype match type, one of FEATURE_{FWD,BCK,MDL}_MATCH
  259. @param matched output as an array of pointers to features with a match of
  260. the specified type
  261. @return Returns the number of features output in matched.
  262. */
  263. int get_matched_features( struct feature* features, int n, int mtype,
  264. struct feature*** matched )
  265. {
  266. struct feature** _matched;
  267. struct ransac_data* rdata;
  268. int i, m = 0;
  269. _matched = calloc( n, sizeof( struct feature* ) );
  270. for( i = 0; i < n; i++ )
  271. if( get_match( features + i, mtype ) )
  272. {
  273. rdata = malloc( sizeof( struct ransac_data ) );
  274. memset( rdata, 0, sizeof( struct ransac_data ) );
  275. rdata->orig_feat_data = features[i].feature_data;
  276. _matched[m] = features + i;
  277. _matched[m]->feature_data = rdata;
  278. m++;
  279. }
  280. *matched = _matched;
  281. return m;
  282. }
  283. /*
  284. Calculates the minimum number of inliers as a function of the number of
  285. putative correspondences.  Based on equation (7) in
  286. Chum, O. and Matas, J.  Matching with PROSAC -- Progressive Sample Consensus.
  287. In <EM>Conference on Computer Vision and Pattern Recognition (CVPR)</EM>,
  288. (2005), pp. 220--226.
  289. @param n number of putative correspondences
  290. @param m min number of correspondences to compute the model in question
  291. @param p_badsupp prob. that a bad model is supported by a correspondence
  292. @param p_badxform desired prob. that the final transformation returned is bad
  293. @return Returns the minimum number of inliers required to guarantee, based
  294. on p_badsupp, that the probability that the final transformation returned
  295. by RANSAC is less than p_badxform
  296. */
  297. int calc_min_inliers( int n, int m, double p_badsupp, double p_badxform )
  298. {
  299. double pi, sum;
  300. int i, j;
  301. for( j = m+1; j <= n; j++ )
  302. {
  303. sum = 0;
  304. for( i = j; i <= n; i++ )
  305. {
  306. pi = pow( p_badsupp, i - m ) * pow( 1.0 - p_badsupp, n - i + m ) *
  307. gsl_sf_choose( n - m, i - m );
  308. sum += pi;
  309. }
  310. if( sum < p_badxform )
  311. break;
  312. }
  313. return j;
  314. }
  315. /*
  316. Draws a RANSAC sample from a set of features.
  317. @param features array of pointers to features from which to sample
  318. @param n number of features in features
  319. @param m size of the sample
  320. @param rng random number generator used to sample
  321. @return Returns an array of pointers to the sampled features; the sampled
  322. field of each sampled feature's ransac_data is set to 1
  323. */
  324. struct feature** draw_ransac_sample( struct feature** features, int n,
  325. int m, gsl_rng* rng )
  326. {
  327. struct feature** sample, * feat;
  328. struct ransac_data* rdata;
  329. int i, x;
  330. for( i = 0; i < n; i++ )
  331. {
  332. rdata = feat_ransac_data( features[i] );
  333. rdata->sampled = 0;
  334. }
  335. sample = calloc( m, sizeof( struct feature* ) );
  336. for( i = 0; i < m; i++ )
  337. {
  338. do
  339. {
  340. x = gsl_rng_uniform_int( rng, n );
  341. feat = features[x];
  342. rdata = feat_ransac_data( feat );
  343. }
  344. while( rdata->sampled );
  345. sample[i] = feat;
  346. rdata->sampled = 1;
  347. }
  348. return sample;
  349. }
  350. /*
  351. Extrancs raw point correspondence locations from a set of features
  352. @param features array of features from which to extract points and match
  353. points; each of these is assumed to have a match of type mtype
  354. @param n number of features
  355. @param mtype match type; if FEATURE_MDL_MATCH correspondences are assumed
  356. to be between each feature's img_pt field and it's match's mdl_pt field,
  357. otherwise, correspondences are assumed to be between img_pt and img_pt
  358. @param pts output as an array of raw point locations from features
  359. @param mpts output as an array of raw point locations from features' matches
  360. */
  361. void extract_corresp_pts( struct feature** features, int n, int mtype,
  362.  CvPoint2D64f** pts, CvPoint2D64f** mpts )
  363. {
  364. struct feature* match;
  365. CvPoint2D64f* _pts, * _mpts;
  366. int i;
  367. _pts = calloc( n, sizeof( CvPoint2D64f ) );
  368. _mpts = calloc( n, sizeof( CvPoint2D64f ) );
  369. if( mtype == FEATURE_MDL_MATCH )
  370. for( i = 0; i < n; i++ )
  371. {
  372. match = get_match( features[i], mtype );
  373. if( ! match )
  374. fatal_error( "feature does not have match of type %d, %s line %d",
  375. mtype, __FILE__, __LINE__ );
  376. _pts[i] = features[i]->img_pt;
  377. _mpts[i] = match->mdl_pt;
  378. }
  379. else
  380. for( i = 0; i < n; i++ )
  381. {
  382. match = get_match( features[i], mtype );
  383. if( ! match )
  384. fatal_error( "feature does not have match of type %d, %s line %d",
  385. mtype, __FILE__, __LINE__ );
  386. _pts[i] = features[i]->img_pt;
  387. _mpts[i] = match->img_pt;
  388. }
  389. *pts = _pts;
  390. *mpts = _mpts;
  391. }
  392. /*
  393. For a given model and error function, finds a consensus from a set of
  394. feature correspondences.
  395. @param features set of pointers to features; every feature is assumed to
  396. have a match of type mtype
  397. @param n number of features in features
  398. @param mtype determines the match field of each feature against which to
  399. measure error; if this is FEATURE_MDL_MATCH, correspondences are assumed
  400. to be between the feature's img_pt field and the match's mdl_pt field;
  401. otherwise matches are assumed to be between img_pt and img_pt
  402. @param M model for which a consensus set is being found
  403. @param err_fn error function used to measure distance from M
  404. @param err_tol correspondences within this distance of M are added to the
  405. consensus set
  406. @param consensus output as an array of pointers to features in the
  407. consensus set
  408. @return Returns the number of points in the consensus set
  409. */
  410. int find_consensus( struct feature** features, int n, int mtype,
  411.    CvMat* M, ransac_err_fn err_fn, double err_tol,
  412.    struct feature*** consensus )
  413. {
  414. struct feature** _consensus;
  415. struct feature* match;
  416. CvPoint2D64f pt, mpt;
  417. double err;
  418. int i, in = 0;
  419. _consensus = calloc( n, sizeof( struct feature* ) );
  420. if( mtype == FEATURE_MDL_MATCH )
  421. for( i = 0; i < n; i++ )
  422. {
  423. match = get_match( features[i], mtype );
  424. if( ! match )
  425. fatal_error( "feature does not have match of type %d, %s line %d",
  426. mtype, __FILE__, __LINE__ );
  427. pt = features[i]->img_pt;
  428. mpt = match->mdl_pt;
  429. err = err_fn( pt, mpt, M );
  430. if( err <= err_tol )
  431. _consensus[in++] = features[i];
  432. }
  433. else
  434. for( i = 0; i < n; i++ )
  435. {
  436. match = get_match( features[i], mtype );
  437. if( ! match )
  438. fatal_error( "feature does not have match of type %d, %s line %d",
  439. mtype, __FILE__, __LINE__ );
  440. pt = features[i]->img_pt;
  441. mpt = match->img_pt;
  442. err = err_fn( pt, mpt, M );
  443. if( err <= err_tol )
  444. _consensus[in++] = features[i];
  445. }
  446. *consensus = _consensus;
  447. return in;
  448. }
  449. /*
  450. Releases memory and reduces code size above
  451. @param pts1 an array of points
  452. @param pts2 an array of points
  453. @param features an array of pointers to features; can be NULL
  454. */
  455. static __inline void release_mem( CvPoint2D64f* pts1, CvPoint2D64f* pts2,
  456. struct feature** features )
  457. {
  458. free( pts1 );
  459. free( pts2 );
  460. if( features )
  461. free( features );
  462. }