-
Notifications
You must be signed in to change notification settings - Fork 5
/
compute_asift_matches.cpp
699 lines (600 loc) · 22.3 KB
/
compute_asift_matches.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
// Copyright (c) 2008-2011, Guoshen Yu <yu@cmap.polytechnique.fr>
// Copyright (c) 2008-2011, Jean-Michel Morel <morel@cmla.ens-cachan.fr>
//
// WARNING:
// This file implements an algorithm possibly linked to the patent
//
// Jean-Michel Morel and Guoshen Yu, Method and device for the invariant
// affine recognition recognition of shapes (WO/2009/150361), patent pending.
//
// This file is made available for the exclusive aim of serving as
// scientific tool to verify of the soundness and
// completeness of the algorithm description. Compilation,
// execution and redistribution of this file may violate exclusive
// patents rights in certain countries.
// The situation being different for every country and changing
// over time, it is your responsibility to determine which patent
// rights restrictions apply to you before you compile, use,
// modify, or redistribute this file. A patent lawyer is qualified
// to make this determination.
// If and only if they don't conflict with any patent terms, you
// can benefit from the following license terms attached to this
// file.
//
// This program is provided for scientific and educational only:
// you can use and/or modify it for these purposes, but you are
// not allowed to redistribute this work or derivative works in
// source or executable form. A license must be obtained from the
// patent right holders for any other use.
//
//
//*------------------------ compute_asift_matches-- -------------------------*/
// Match the ASIFT keypoints.
//
// Please report bugs and/or send comments to Guoshen Yu yu@cmap.polytechnique.fr
//
// Reference: J.M. Morel and G.Yu, ASIFT: A New Framework for Fully Affine Invariant Image
// Comparison, SIAM Journal on Imaging Sciences, vol. 2, issue 2, pp. 438-469, 2009.
// Reference: ASIFT online demo (You can try ASIFT with your own images online.)
// http://www.ipol.im/pub/algo/my_affine_sift/
/*---------------------------------------------------------------------------*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#ifdef _OPENMP
#include <omp.h>
#endif
#include "compute_asift_matches.h"
#include "lib/match/match.h"
#include "orsa.h"
#define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5))
/* Remove the repetitive matches that appear in different simulations and retain only one */
void unique_match1(matchingslist & seg_in, matchingslist & seg_out, vector < vector < float > >&Minfoall_in,
vector < vector < float > >&Minfoall_out)
{
int i_in, i_out;
float x1_in, x2_in, y1_in, y2_in, x1_out, x2_out, y1_out, y2_out;
int flag_unique;
float d1, d2;
int Th2 = 2;
seg_out.push_back(seg_in[0]);
Minfoall_out.push_back(Minfoall_in[0]);
/* For other matches */
if (seg_in.size() > 1) {
/* check if a match is unique. if yes, copy */
matchingslist::iterator ptr_in = seg_in.begin();
for (i_in = 1; i_in < (int)seg_in.size(); i_in++, ptr_in++) {
x1_in = ptr_in->first.x;
y1_in = ptr_in->first.y;
x2_in = ptr_in->second.x;
y2_in = ptr_in->second.y;
flag_unique = 1;
matchingslist::iterator ptr_out = seg_out.begin();
for (i_out = 0; i_out < (int)seg_out.size(); i_out++, ptr_out++) {
x1_out = ptr_out->first.x;
y1_out = ptr_out->first.y;
x2_out = ptr_out->second.x;
y2_out = ptr_out->second.y;
d1 = (x1_in - x1_out) * (x1_in - x1_out) + (y1_in - y1_out) * (y1_in - y1_out);
d2 = (x2_in - x2_out) * (x2_in - x2_out) + (y2_in - y2_out) * (y2_in - y2_out);
if ((d1 <= Th2) && (d2 <= Th2)) {
flag_unique = 0;
continue;
}
}
if (flag_unique == 1) {
seg_out.push_back(seg_in[i_in]);
Minfoall_out.push_back(Minfoall_in[i_in]);
}
}
}
}
/* Remove the ALL one-to-multiple matches. */
void clean_match1(matchingslist & seg_in, matchingslist & seg_out, vector < vector < float > >&Minfoall_in,
vector < vector < float > >&Minfoall_out)
{
int i1, i2;
float x1_in, x2_in, y1_in, y2_in, x1_out, x2_out, y1_out, y2_out;
// Guoshen Yu, 2010.09.22, Windows version
// int flag_unique[seg_in.size()];
int tmp_size = seg_in.size();
int *flag_unique = new int[tmp_size];
int sum_flag = 0;
float d1, d2;
int Th1 = 1;
int Th2 = 4;
for (i1 = 0; i1 < (int)seg_in.size(); i1++) {
flag_unique[i1] = 1;
}
/* Set the flag of redundant matches to 0. */
matchingslist::iterator ptr_in = seg_in.begin();
for (i1 = 0; i1 < (int)seg_in.size() - 1; i1++, ptr_in++) {
x1_in = ptr_in->first.x;
y1_in = ptr_in->first.y;
x2_in = ptr_in->second.x;
y2_in = ptr_in->second.y;
matchingslist::iterator ptr_out = ptr_in + 1;
for (i2 = i1 + 1; i2 < (int)seg_in.size(); i2++, ptr_out++) {
x1_out = ptr_out->first.x;
y1_out = ptr_out->first.y;
x2_out = ptr_out->second.x;
y2_out = ptr_out->second.y;
d1 = (x1_in - x1_out) * (x1_in - x1_out) + (y1_in - y1_out) * (y1_in - y1_out);
d2 = (x2_in - x2_out) * (x2_in - x2_out) + (y2_in - y2_out) * (y2_in - y2_out);
/* If redundant, set flags of both elements to 0. */
if (((d1 <= Th1) && (d2 > Th2)) || ((d1 > Th2) && (d2 <= Th1))) {
flag_unique[i1] = 0;
flag_unique[i2] = 0;
}
}
}
for (i1 = 0; i1 < (int)seg_in.size(); i1++) {
sum_flag += flag_unique[i1];
}
/* Copy the matches that are not redundant */
if (sum_flag > 0) {
for (i1 = 0; i1 < (int)seg_in.size(); i1++) {
if (flag_unique[i1] == 1) {
seg_out.push_back(seg_in[i1]);
Minfoall_out.push_back(Minfoall_in[i1]);
}
}
} else {
printf
("Warning: all matches are redundant and are thus removed! This step of match cleaning is short circuited. (Normally this should not happen...)\n");
}
// Guoshen Yu, 2010.09.22, Windows version
delete[]flag_unique;
}
/* Remove the ALL multiple-to-one matches */
void clean_match2(matchingslist & seg_in, matchingslist & seg_out, vector < vector < float > >&Minfoall_in,
vector < vector < float > >&Minfoall_out)
{
int i1, i2;
float x1_in, x2_in, y1_in, y2_in, x1_out, x2_out, y1_out, y2_out;
// Guoshen Yu, 2010.09.22, Windows version
// int flag_unique[seg_in.size()];
int tmp_size = seg_in.size();
int *flag_unique = new int[tmp_size];
int sum_flag = 0;
float d1, d2;
int Th1 = 1;
int Th2 = 4;
for (i1 = 0; i1 < (int)seg_in.size(); i1++) {
flag_unique[i1] = 1;
}
/* Set the flag of redundant matches to 0. */
matchingslist::iterator ptr_in = seg_in.begin();
for (i1 = 0; i1 < (int)seg_in.size() - 1; i1++, ptr_in++) {
x1_in = ptr_in->first.x;
y1_in = ptr_in->first.y;
x2_in = ptr_in->second.x;
y2_in = ptr_in->second.y;
matchingslist::iterator ptr_out = ptr_in + 1;
for (i2 = i1 + 1; i2 < (int)seg_in.size(); i2++, ptr_out++) {
x1_out = ptr_out->first.x;
y1_out = ptr_out->first.y;
x2_out = ptr_out->second.x;
y2_out = ptr_out->second.y;
d1 = (x1_in - x1_out) * (x1_in - x1_out) + (y1_in - y1_out) * (y1_in - y1_out);
d2 = (x2_in - x2_out) * (x2_in - x2_out) + (y2_in - y2_out) * (y2_in - y2_out);
/* If redundant, set flags of both elements to 0. */
if ((d1 > Th2) && (d2 <= Th1)) {
flag_unique[i1] = 0;
flag_unique[i2] = 0;
}
}
}
for (i1 = 0; i1 < (int)seg_in.size(); i1++) {
sum_flag += flag_unique[i1];
}
/* Copy the matches that are not redundant */
if (sum_flag > 0) {
for (i1 = 0; i1 < (int)seg_in.size(); i1++) {
if (flag_unique[i1] == 1) {
seg_out.push_back(seg_in[i1]);
Minfoall_out.push_back(Minfoall_in[i1]);
}
}
} else {
printf
("Warning: all matches are redundant and are thus removed! This step of match cleaning is short circuited. (Normally this should not happen...)\n");
}
// Guoshen Yu, 2010.09.22, Windows version
delete[]flag_unique;
}
// Normalize the coordinates of the matched points by compensating the simulate affine transformations
void compensate_affine_coor(matching & matching1, int w1, int h1, int w2, int h2, float t1, float t2, float Rtheta,
float t_im2_1, float t_im2_2, float Rtheta2)
{
float x_ori, y_ori;
float x_ori2, y_ori2, x_tmp, y_tmp;
float x1, y1, x2, y2;
Rtheta = Rtheta * PI / 180;
if (Rtheta <= PI / 2) {
x_ori = 0;
y_ori = w1 * sin(Rtheta) / t1;
} else {
x_ori = -w1 * cos(Rtheta) / t2;
y_ori = (w1 * sin(Rtheta) + h1 * sin(Rtheta - PI / 2)) / t1;
}
Rtheta2 = Rtheta2 * PI / 180;
if (Rtheta2 <= PI / 2) {
x_ori2 = 0;
y_ori2 = w2 * sin(Rtheta2) / t_im2_1;
} else {
x_ori2 = -w2 * cos(Rtheta2) / t_im2_2;
y_ori2 = (w2 * sin(Rtheta2) + h2 * sin(Rtheta2 - PI / 2)) / t_im2_1;
}
float sin_Rtheta = sin(Rtheta);
float cos_Rtheta = cos(Rtheta);
float sin_Rtheta2 = sin(Rtheta2);
float cos_Rtheta2 = cos(Rtheta2);
x1 = matching1.first.x;
y1 = matching1.first.y;
x2 = matching1.second.x;
y2 = matching1.second.y;
/* project the coordinates of im1 to original image before tilt-rotation transform */
/* Get the coordinates with respect to the 'origin' of the original image before transform */
x1 = x1 - x_ori;
y1 = y1 - y_ori;
/* Invert tilt */
x1 = x1 * t2;
y1 = y1 * t1;
/* Invert rotation (Note that the y direction (vertical) is inverse to the usual concention. Hence Rtheta instead of -Rtheta to inverse the rotation.) */
x_tmp = cos_Rtheta * x1 - sin_Rtheta * y1;
y_tmp = sin_Rtheta * x1 + cos_Rtheta * y1;
x1 = x_tmp;
y1 = y_tmp;
/* Coordinate projection on image2 */
/* Get the coordinates with respect to the 'origin' of the original image before transform */
x2 = x2 - x_ori2;
y2 = y2 - y_ori2;
/* Invert tilt */
x2 = x2 * t_im2_2;
y2 = y2 * t_im2_1;
/* Invert rotation (Note that the y direction (vertical) is inverse to the usual concention. Hence Rtheta instead of -Rtheta to inverse the rotation.) */
x_tmp = cos_Rtheta2 * x2 - sin_Rtheta2 * y2;
y_tmp = sin_Rtheta2 * x2 + cos_Rtheta2 * y2;
x2 = x_tmp;
y2 = y_tmp;
matching1.first.x = x1;
matching1.first.y = y1;
matching1.second.x = x2;
matching1.second.y = y2;
}
int compute_asift_matches(int num_of_tilts1, int num_of_tilts2, int w1, int h1, int w2, int h2, int verb,
vector < vector < keypointslist > >&keys1, vector < vector < keypointslist > >&keys2,
matchingslist & matchings, siftPar & siftparameters)
// Match the ASIFT keypoints.
// Input:
// num_of_tilts1, num_of_tilts2: number of tilts that have been simulated on the two images. (They can be different.)
// w1, h1, w2, h2: widht/height of image1/image2.
// verb: 1/0 --> show/don not show verbose messages. (1 for debugging)
// keys1, keys2: ASIFT keypoints of image1/image2. (They should be calculated with compute_asift_keypoints.)
// matchings (output): the coordinates (col1, row1, col2, row2) of all the matching points.
//
// Output: the number of matching points.
{
float t_min, t_k, t;
int num_tilt1, num_tilt2, tt, num_rot_t2, num_rot1, rr;
int cc;
int tt2, rr2, num_rot1_2;
float t_im2;
/* It stores the coordinates of ALL matches points of ALL affine simulations */
vector < vector < float > >Minfoall;
int Tmin = 8;
float nfa_max = -2;
num_rot_t2 = 10;
t_min = 1;
t_k = sqrt(2.);
num_tilt1 = num_of_tilts1;
num_tilt2 = num_of_tilts2;
if ((num_tilt1 < 1) || (num_tilt2 < 1)) {
printf("Number of tilts num_tilt should be equal or larger than 1. \n");
exit(-1);
}
/* Initialize the vector structure for the matching points */
std::vector < vector < vector < vector < matchingslist > > > >matchings_vec(num_tilt1);
std::vector < vector < vector < vector < vector < vector < float > > > > > >Minfoall_vec(num_tilt1);
for (tt = 1; tt <= num_tilt1; tt++) {
t = t_min * pow(t_k, tt - 1);
if (t == 1) {
num_rot1 = 1;
} else {
num_rot1 = round(num_rot_t2 * t / 2);
if (num_rot1 % 2 == 1) {
num_rot1 = num_rot1 + 1;
}
num_rot1 = num_rot1 / 2;
}
matchings_vec[tt - 1].resize(num_rot1);
Minfoall_vec[tt - 1].resize(num_rot1);
for (rr = 1; rr <= num_rot1; rr++) {
matchings_vec[tt - 1][rr - 1].resize(num_tilt2);
Minfoall_vec[tt - 1][rr - 1].resize(num_tilt2);
for (tt2 = 1; tt2 <= num_tilt2; tt2++) {
t_im2 = t_min * pow(t_k, tt2 - 1);
if (t_im2 == 1) {
num_rot1_2 = 1;
} else {
num_rot1_2 = round(num_rot_t2 * t_im2 / 2);
if (num_rot1_2 % 2 == 1) {
num_rot1_2 = num_rot1_2 + 1;
}
num_rot1_2 = num_rot1_2 / 2;
}
matchings_vec[tt - 1][rr - 1][tt2 - 1].resize(num_rot1_2);
Minfoall_vec[tt - 1][rr - 1][tt2 - 1].resize(num_rot1_2);
}
}
}
///*
// * setup the tilt and rotation parameters
// * for all the loops, this vector will hold
// * the following parameters:
// * tt, num_rot1, rr, tt2, num_rot1_2, rr2
// */
//vector<int> tilt_rot;
///* loop on tilts for image 1 */
//for (int tt = 1; tt <= num_tilt1; tt++)
//{
// float t = t_min * pow(t_k, tt-1);
// int num_rot1;
// /* if tilt t = 1, do not simulate rotation. */
// if ( 1 == tt )
// num_rot1 = 1;
// else
// {
// /* number of rotations to simulate */
// num_rot1 = round(num_rot_t2 * t / 2);
// if ( num_rot1%2 == 1 )
// num_rot1 = num_rot1 + 1;
// num_rot1 = num_rot1 / 2;
// }
// /* loop on rotations for image 1 */
// for (int rr = 1; rr <= num_rot1; rr++ )
// {
// /* loop on tilts for image 2 */
// for (int tt2 = 1; tt2 <= num_tilt2; tt2++)
// {
// float t_im2 = t_min * pow(t_k, tt2-1);
// int num_rot1_2;
// if ( tt2 == 1 )
// num_rot1_2 = 1;
// else
// {
// num_rot1_2 = round(num_rot_t2 * t_im2 / 2);
// if ( num_rot1_2%2 == 1 )
// num_rot1_2 = num_rot1_2 + 1;
// num_rot1_2 = num_rot1_2 / 2;
// }
// /* loop on rotations for image 2 */
// for (int rr2 = 1; rr2 <= num_rot1_2; rr2++ )
// {
// tilt_rot.push_back(tt);
// tilt_rot.push_back(num_rot1);
// tilt_rot.push_back(rr);
// tilt_rot.push_back(tt2);
// tilt_rot.push_back(num_rot1_2);
// tilt_rot.push_back(rr2);
// }
// }
// }
//}
/* Calculate the number of simulations */
#ifdef _OPENMP
omp_set_nested(1);
#endif
// loop on tilts for image 1.
#pragma omp parallel for private(tt)
for (int tt = 1; tt <= num_tilt1; tt++) {
float t = t_min * pow(t_k, tt - 1);
/* Attention: the t1, t2 do not follow the same convention as in compute_asift_keypoints */
float t1 = t;
float t2 = 1;
int num_rot1;
// If tilt t = 1, do not simulate rotation.
if (tt == 1) {
num_rot1 = 1;
} else {
// The number of rotations to simulate under the current tilt.
num_rot1 = round(num_rot_t2 * t / 2);
if (num_rot1 % 2 == 1) {
num_rot1 = num_rot1 + 1;
}
num_rot1 = num_rot1 / 2;
}
float delta_theta = PI / num_rot1;
// Loop on rotations for image 1.
#pragma omp parallel for private(rr)
for (int rr = 1; rr <= num_rot1; rr++) {
float theta = delta_theta * (rr - 1);
theta = theta * 180 / PI;
/* Read the keypoints of image 1 */
keypointslist keypoints1 = keys1[tt - 1][rr - 1];
// loop on tilts for image 2.
#pragma omp parallel for private(tt2)
for (int tt2 = 1; tt2 <= num_tilt2; tt2++) {
float t_im2 = t_min * pow(t_k, tt2 - 1);
/* Attention: the t1, t2 do not follow the same convention as in asift_v1.c */
float t_im2_1 = t_im2;
float t_im2_2 = 1;
int num_rot1_2;
if (tt2 == 1) {
num_rot1_2 = 1;
} else {
num_rot1_2 = round(num_rot_t2 * t_im2 / 2);
if (num_rot1_2 % 2 == 1) {
num_rot1_2 = num_rot1_2 + 1;
}
num_rot1_2 = num_rot1_2 / 2;
}
float delta_theta2 = PI / num_rot1_2;
#pragma omp parallel for private(rr2)
// Loop on rotations for image 2.
for (int rr2 = 1; rr2 <= num_rot1_2; rr2++) {
float theta2 = delta_theta2 * (rr2 - 1);
theta2 = theta2 * 180 / PI;
/* Read the keypoints of image2. */
keypointslist keypoints2 = keys2[tt2 - 1][rr2 - 1];
// Match the keypoints of image1 and image2.
matchingslist matchings1;
compute_sift_matches(keypoints1, keypoints2, matchings1, siftparameters);
if (verb) {
printf
("t1=%.2f, theta1=%.2f, num keys1 = %d, t2=%.2f, theta2=%.2f, num keys2 = %d, num matches=%d\n",
t, theta, (int)keypoints1.size(), t_im2, theta2, (int)keypoints2.size(),
(int)matchings1.size());
}
/* Store the matches */
if (matchings1.size() > 0) {
matchings_vec[tt - 1][rr - 1][tt2 - 1][rr2 - 1] =
matchingslist(matchings1.size());
Minfoall_vec[tt - 1][rr - 1][tt2 - 1][rr2 - 1].resize(matchings1.size());
for (int cc = 0; cc < (int)matchings1.size(); cc++) {
///// In the coordinates the affine transformations have been normalized already in compute_asift_keypoints. So no need to normalize here.
// Normalize the coordinates of the matched points by compensating the simulate affine transformations
// compensate_affine_coor(matchings1[cc], w1, h1, w2, h2, t1, t2, theta, t_im2_1, t_im2_2, theta2);
matchings_vec[tt - 1][rr - 1][tt2 - 1][rr2 - 1][cc] = matchings1[cc];
vector < float >Minfo_1match(6);
Minfo_1match[0] = t1;
Minfo_1match[1] = t2;
Minfo_1match[2] = theta;
Minfo_1match[3] = t_im2_1;
Minfo_1match[4] = t_im2_2;
Minfo_1match[5] = theta2;
Minfoall_vec[tt - 1][rr - 1][tt2 - 1][rr2 - 1][cc] = Minfo_1match;
}
}
}
}
}
}
// Move the matches to a 1D vector
for (tt = 1; tt <= num_tilt1; tt++) {
t = t_min * pow(t_k, tt - 1);
if (t == 1) {
num_rot1 = 1;
} else {
num_rot1 = round(num_rot_t2 * t / 2);
if (num_rot1 % 2 == 1) {
num_rot1 = num_rot1 + 1;
}
num_rot1 = num_rot1 / 2;
}
for (rr = 1; rr <= num_rot1; rr++) {
for (tt2 = 1; tt2 <= num_tilt2; tt2++) {
t_im2 = t_min * pow(t_k, tt2 - 1);
if (t_im2 == 1) {
num_rot1_2 = 1;
} else {
num_rot1_2 = round(num_rot_t2 * t_im2 / 2);
if (num_rot1_2 % 2 == 1) {
num_rot1_2 = num_rot1_2 + 1;
}
num_rot1_2 = num_rot1_2 / 2;
}
for (rr2 = 1; rr2 <= num_rot1_2; rr2++) {
for (cc = 0; cc < (int)matchings_vec[tt - 1][rr - 1][tt2 - 1][rr2 - 1].size(); cc++) {
matchings.push_back(matchings_vec[tt - 1][rr - 1][tt2 - 1][rr2 - 1][cc]);
Minfoall.push_back(Minfoall_vec[tt - 1][rr - 1][tt2 - 1][rr2 - 1][cc]);
}
}
}
}
}
if (verb) {
printf("The number of matches is %d \n", (int)matchings.size());
}
if (matchings.size() > 0) {
/* Remove the repetitive matches that appear in different simulations and retain only one. */
// Since tilts are simuated on both image 1 and image 2, it is normal to have repetitive matches.
matchingslist matchings_unique;
vector < vector < float > >Minfoall_unique;
unique_match1(matchings, matchings_unique, Minfoall, Minfoall_unique);
matchings = matchings_unique;
Minfoall = Minfoall_unique;
if (verb) {
printf("The number of unique matches is %d \n", (int)matchings.size());
}
// There often appear to be some one-to-multiple/multiple-to-one matches (one point in image 1 matches with many points in image 2/vice versa).
// This is an artifact of SIFT on interpolated images, as the interpolation tends to create some auto-similar structures (steps for example).
// These matches need to be removed.
/* Separating the removal of multiple-to-one and one-to-multiple in two steps:
- first remove multiple-to-one
- then remove one-to-multiple
This allows to avoid removing some good matches: multiple-to-one matches is much more frequent than one-to-multiple. Sometimes some of the feature points in image 1 that take part in "multiple-to-one" bad matches have also correct matches in image 2. The modified scheme avoid removing these good matches. */
// Remove to multiple-to-one matches
matchings_unique.clear();
Minfoall_unique.clear();
clean_match2(matchings, matchings_unique, Minfoall, Minfoall_unique);
matchings = matchings_unique;
Minfoall = Minfoall_unique;
// Remove to one-to-multiple matches
matchings_unique.clear();
Minfoall_unique.clear();
clean_match1(matchings, matchings_unique, Minfoall, Minfoall_unique);
matchings = matchings_unique;
Minfoall = Minfoall_unique;
if (verb) {
printf("The number of final matches is %d \n", (int)matchings.size());
}
// If enough matches to do epipolar filtering
if ((int)matchings.size() >= Tmin) {
//////// Use ORSA to filter out the incorrect matches.
// store the coordinates of the matching points
vector < Match > match_coor;
for (cc = 0; cc < (int)matchings.size(); cc++) {
Match match1_coor;
match1_coor.x1 = matchings[cc].first.x;
match1_coor.y1 = matchings[cc].first.y;
match1_coor.x2 = matchings[cc].second.x;
match1_coor.y2 = matchings[cc].second.y;
match_coor.push_back(match1_coor);
}
std::vector < float >index;
// Guoshen Yu, 2010.09.23
// index.clear();
int t_value_orsa = 10000;
int verb_value_orsa = 0;
int n_flag_value_orsa = 0;
int mode_value_orsa = 2;
int stop_value_orsa = 0;
// epipolar filtering with the Moisan-Stival ORSA algorithm.
// float nfa = orsa(w1, h1, match_coor, index, t_value_orsa, verb_value_orsa, n_flag_value_orsa, mode_value_orsa, stop_value_orsa);
float nfa =
orsa((w1 + w2) / 2, (h1 + h2) / 2, match_coor, index, t_value_orsa, verb_value_orsa,
n_flag_value_orsa, mode_value_orsa, stop_value_orsa);
// if the matching is significant, register the good matches
if (nfa < nfa_max) {
// extract meaningful matches
matchings_unique.clear();
Minfoall_unique.clear();
for (cc = 0; cc < (int)index.size(); cc++) {
matchings_unique.push_back(matchings[(int)index[cc]]);
Minfoall_unique.push_back(Minfoall[(int)index[cc]]);
}
matchings = matchings_unique;
Minfoall = Minfoall_unique;
cout << "The two images match! " << matchings.
size() << " matchings are identified. log(nfa)=" << nfa << "." << endl;
return 0;
} else {
matchings.clear();
Minfoall.clear();
cout << "The two images do not match. The matching is not significant: log(nfa)=" << nfa <<
"." << endl;
}
} else {
matchings.clear();
Minfoall.clear();
cout << "The two images do not match. Not enough matches to do epipolar filtering." << endl;
}
} else {
cout << "The two images do not match.\n" << endl;
}
return -1;
}