diff --git a/include/vips/proto.h b/include/vips/proto.h index 65a917a7..9cbb78a2 100644 --- a/include/vips/proto.h +++ b/include/vips/proto.h @@ -665,6 +665,7 @@ int im_remosaic( IMAGE *in, IMAGE *out, const char *old_str, const char *new_str ); int im_align_bands( IMAGE *in, IMAGE *out ); +int im_maxpos_subpel( IMAGE *in, double *x, double *y ); /* inplace */ diff --git a/libsrc/mosaicing/Makefile.am b/libsrc/mosaicing/Makefile.am index 907553d0..7ebc4c2c 100644 --- a/libsrc/mosaicing/Makefile.am +++ b/libsrc/mosaicing/Makefile.am @@ -16,6 +16,7 @@ libmosaicing_la_SOURCES = \ im_lrcalcon.c \ im_lrmerge.c \ im_lrmosaic.c \ + im_maxpos_subpel.c \ im_tbcalcon.c \ im_tbmerge.c \ im_remosaic.c \ diff --git a/libsrc/mosaicing/im_maxpos_subpel.c b/libsrc/mosaicing/im_maxpos_subpel.c new file mode 100644 index 00000000..51c0c5fc --- /dev/null +++ b/libsrc/mosaicing/im_maxpos_subpel.c @@ -0,0 +1,157 @@ +/* This function implements: + * "Extension of Phase Correlation to Subpixel Registration" + * by H. Foroosh, from IEEE trans. Im. Proc. 11(3), 2002. + * + * If the best three matches in the correlation are aranged: + * + * 02 or 01 + * 1 2 + * + * then we return a subpixel match using the ratio of correlations in the + * vertical and horizontal dimension. + * + * ( xs[0], ys[0] ) is the best integer alignment + * ( xs[ use_x ], ys[ use_x ] ) is equal in y and (+/-)1 off in x + * ( xs[ use_y ], ys[ use_y ] ) is equal in x and (+/-)1 off in y + * + * + * Alternatively if the best four matches in the correlation are aranged in + * a square: + * + * 01 or 03 or 02 or 03 + * 32 12 31 21 + * + * then we return a subpixel match weighting with the sum the two on each + * side over the sum of all four, but only if all four of them are very + * close to the best, and the fifth is nowhere near. + * + * This alternative method is not described by Foroosh, but is often the + * case where the match is close to n-and-a-half pixels in both dimensions. + * + * Copyright: 2008, Nottingham Trent University + * + * Author: Tom Vajzovic + * + * Written on: 2008-02-07 + */ + +/* + + This file is part of VIPS. + + VIPS is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + */ + +/* + + These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk + + */ + +#ifdef HAVE_CONFIG_H +#include +#endif /*HAVE_CONFIG_H */ +#include + +#include +#include + +#ifdef WITH_DMALLOC +#include +#endif /*WITH_DMALLOC */ + + +#define MOST_OF( A, B ) ( (A) > 0.9 * (B) ) +#define LITTLE_OF( A, B ) ( (B) < 0.1 * (B) ) + +int im_maxpos_subpel( IMAGE *in, double *x, double *y ){ +#define FUNCTION_NAME "im_maxpos_subpel" + + int xs[5]; + int ys[5]; + double vals[5]; + int xa, ya, xb, yb; + double vxa, vya, vxb, vyb; + + if( im_maxpos_vec( in, xs, ys, vals, 5 )) + return -1; + +#define WRAP_TEST_RETURN() \ + \ + /* wrap around if we have alignment -1 < d <= 0 */ \ + /* (change it to: size - 1 <= d < size ) */ \ + \ + if( ! xa && in-> Xsize - 1 == xb ) \ + xa= in-> Xsize; \ + \ + else if( ! xb && in-> Xsize - 1 == xa ) \ + xb= in-> Xsize; \ + \ + if( ! ya && in-> Ysize - 1 == yb ) \ + ya= in-> Ysize; \ + \ + else if( ! yb && in-> Ysize - 1 == ya ) \ + yb= in-> Ysize; \ + \ + if( 1 == abs( xb - xa ) && 1 == abs( yb - ya )){ \ + *x= ((double)xa) + ((double)( xb - xa )) * ( vxb / ( vxa + vxb )); \ + *y= ((double)ya) + ((double)( yb - ya )) * ( vyb / ( vya + vyb )); \ + return 0; \ + } + +#define TEST3( A, B ) \ + if( xs[0] == xs[A] && ys[0] == ys[B] ){ \ + xa= xs[0]; \ + ya= ys[0]; \ + xb= xs[B]; \ + yb= ys[A]; \ + vxa= vals[0]; \ + vya= vals[0]; \ + vxb= vals[B]; \ + vyb= vals[A]; \ + WRAP_TEST_RETURN() \ + } + + TEST3( 1, 2 ) + TEST3( 2, 1 ) + + if( MOST_OF( vals[1], vals[0] ) && MOST_OF( vals[2], vals[0] ) + && MOST_OF( vals[3], vals[0] ) && LITTLE_OF( vals[4], vals[0] )){ + +#define TEST4( A, B, C, D, E, F, G, H ) \ + if( xs[A] == xs[B] && xs[C] == xs[D] && ys[E] == ys[F] && ys[G] == ys[H] ){ \ + xa= xs[A]; \ + xb= xs[C]; \ + ya= ys[E]; \ + yb= ys[G]; \ + vxa= vals[A] + vals[B]; \ + vxb= vals[C] + vals[D]; \ + vya= vals[E] + vals[F]; \ + vyb= vals[G] + vals[H]; \ + WRAP_TEST_RETURN() \ + } + + TEST4( 0, 3, 1, 2, 0, 1, 2, 3 ) + TEST4( 0, 1, 2, 3, 0, 3, 1, 2 ) + TEST4( 0, 3, 1, 2, 0, 2, 1, 3 ) + TEST4( 0, 2, 1, 3, 0, 3, 1, 2 ) + } + + im_warn( FUNCTION_NAME, "registration performed to nearest pixel only: correlation does not have the expected distribution for sub-pixel registration" ); + *x= (double) xs[0]; + *y= (double) ys[0]; + return 0; +} diff --git a/libsrc/mosaicing/mosaicing_dispatch.c b/libsrc/mosaicing/mosaicing_dispatch.c index a19e612e..161ba4db 100644 --- a/libsrc/mosaicing/mosaicing_dispatch.c +++ b/libsrc/mosaicing/mosaicing_dispatch.c @@ -832,6 +832,25 @@ static im_function align_bands_desc= { align_bands_arg_types }; +static int maxpos_subpel_vec( im_object *argv ){ + return im_maxpos_subpel( (IMAGE*)argv[0], (double*)argv[1], (double*)argv[2] ); +} + +static im_arg_desc maxpos_subpel_arg_types[]= { + IM_INPUT_IMAGE( "im" ), + IM_OUTPUT_DOUBLE( "x" ), + IM_OUTPUT_DOUBLE( "y" ) +}; + +static im_function maxpos_subpel_desc= { + "im_maxpos_subpel", + "subpixel position of maximum of (phase correlation) image", + IM_FN_PIO, + maxpos_subpel_vec, + IM_NUMBER( maxpos_subpel_arg_types ), + maxpos_subpel_arg_types +}; + /* Package up all these functions. */ static im_function *mos_list[] = { @@ -848,6 +867,7 @@ static im_function *mos_list[] = { &lrmosaic1_desc, &match_linear_desc, &match_linear_search_desc, + &maxpos_subpel_desc, &remosaic_desc, &similarity_area_desc, &similarity_desc,