1 /* find position of maximum, subpixel estimation
2 *
3 * Copyright: 2008, Nottingham Trent University
4 * Author: Tom Vajzovic
5 * Written on: 2008-02-07
6 *
7 * 25/1/11
8 * - gtk-doc
9 */
10
11 /*
12
13 This file is part of VIPS.
14
15 VIPS is free software; you can redistribute it and/or modify
16 it under the terms of the GNU Lesser General Public License as published by
17 the Free Software Foundation; either version 2 of the License, or
18 (at your option) any later version.
19
20 This program is distributed in the hope that it will be useful,
21 but WITHOUT ANY WARRANTY; without even the implied warranty of
22 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 GNU Lesser General Public License for more details.
24
25 You should have received a copy of the GNU Lesser General Public License
26 along with this program; if not, write to the Free Software
27 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
28 02110-1301 USA
29
30 */
31
32 /*
33
34 These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
35
36 */
37
38 #ifdef HAVE_CONFIG_H
39 #include <config.h>
40 #endif /*HAVE_CONFIG_H */
41 #include <vips/intl.h>
42
43 #include <stdlib.h>
44 #include <vips/vips.h>
45 #include <vips/vips7compat.h>
46
47
48 #define MOST_OF( A, B ) ( (A) > 0.9 * (B) )
49 #define LITTLE_OF( A, B ) ( (B) < 0.1 * (B) )
50
51 /**
52 * im_maxpos_subpel:
53 * @in: input image
54 * @x: output position of maximum
55 * @y: output position of maximum
56 *
57 * This function implements:
58 *
59 * "Extension of Phase Correlation to Subpixel Registration"
60 * by H. Foroosh, from IEEE trans. Im. Proc. 11(3), 2002.
61 *
62 * If the best three matches in the correlation are aranged:
63 *
64 * 02 or 01
65 * 1 2
66 *
67 * then we return a subpixel match using the ratio of correlations in the
68 * vertical and horizontal dimension.
69 *
70 * ( xs[0], ys[0] ) is the best integer alignment
71 * ( xs[ use_x ], ys[ use_x ] ) is equal in y and (+/-)1 off in x
72 * ( xs[ use_y ], ys[ use_y ] ) is equal in x and (+/-)1 off in y
73 *
74 * Alternatively if the best four matches in the correlation are aranged in
75 * a square:
76 *
77 * 01 or 03 or 02 or 03
78 * 32 12 31 21
79 *
80 * then we return a subpixel match weighting with the sum the two on each
81 * side over the sum of all four, but only if all four of them are very
82 * close to the best, and the fifth is nowhere near.
83 *
84 * This alternative method is not described by Foroosh, but is often the
85 * case where the match is close to n-and-a-half pixels in both dimensions.
86 *
87 * See also: im_maxpos(), im_min(), im_stats().
88 *
89 * Returns: 0 on success, -1 on error
90 */
91
im_maxpos_subpel(IMAGE * in,double * x,double * y)92 int im_maxpos_subpel( IMAGE *in, double *x, double *y ){
93 #define FUNCTION_NAME "im_maxpos_subpel"
94
95 int xs[5];
96 int ys[5];
97 double vals[5];
98 int xa, ya, xb, yb;
99 double vxa, vya, vxb, vyb;
100
101 if( im_maxpos_vec( in, xs, ys, vals, 5 ))
102 return -1;
103
104 #define WRAP_TEST_RETURN() \
105 \
106 /* wrap around if we have alignment -1 < d <= 0 */ \
107 /* (change it to: size - 1 <= d < size ) */ \
108 \
109 if( ! xa && in-> Xsize - 1 == xb ) \
110 xa= in-> Xsize; \
111 \
112 else if( ! xb && in-> Xsize - 1 == xa ) \
113 xb= in-> Xsize; \
114 \
115 if( ! ya && in-> Ysize - 1 == yb ) \
116 ya= in-> Ysize; \
117 \
118 else if( ! yb && in-> Ysize - 1 == ya ) \
119 yb= in-> Ysize; \
120 \
121 if( 1 == abs( xb - xa ) && 1 == abs( yb - ya )){ \
122 *x= ((double)xa) + ((double)( xb - xa )) * ( vxb / ( vxa + vxb )); \
123 *y= ((double)ya) + ((double)( yb - ya )) * ( vyb / ( vya + vyb )); \
124 return 0; \
125 }
126
127 #define TEST3( A, B ) \
128 if( xs[0] == xs[A] && ys[0] == ys[B] ){ \
129 xa= xs[0]; \
130 ya= ys[0]; \
131 xb= xs[B]; \
132 yb= ys[A]; \
133 vxa= vals[0]; \
134 vya= vals[0]; \
135 vxb= vals[B]; \
136 vyb= vals[A]; \
137 WRAP_TEST_RETURN() \
138 }
139
140 TEST3( 1, 2 )
141 TEST3( 2, 1 )
142
143 if( MOST_OF( vals[1], vals[0] ) && MOST_OF( vals[2], vals[0] )
144 && MOST_OF( vals[3], vals[0] ) && LITTLE_OF( vals[4], vals[0] )){
145
146 #define TEST4( A, B, C, D, E, F, G, H ) \
147 if( xs[A] == xs[B] && xs[C] == xs[D] && ys[E] == ys[F] && ys[G] == ys[H] ){ \
148 xa= xs[A]; \
149 xb= xs[C]; \
150 ya= ys[E]; \
151 yb= ys[G]; \
152 vxa= vals[A] + vals[B]; \
153 vxb= vals[C] + vals[D]; \
154 vya= vals[E] + vals[F]; \
155 vyb= vals[G] + vals[H]; \
156 WRAP_TEST_RETURN() \
157 }
158
159 TEST4( 0, 3, 1, 2, 0, 1, 2, 3 )
160 TEST4( 0, 1, 2, 3, 0, 3, 1, 2 )
161 TEST4( 0, 3, 1, 2, 0, 2, 1, 3 )
162 TEST4( 0, 2, 1, 3, 0, 3, 1, 2 )
163 }
164
165 im_warn( FUNCTION_NAME, "registration performed to nearest pixel only: correlation does not have the expected distribution for sub-pixel registration" );
166 *x= (double) xs[0];
167 *y= (double) ys[0];
168 return 0;
169 }
170