/* @(#) im_affine() ... affine transform, bi-linear interpolation. * @(#) * @(#) int im_affine(in, out, a, b, c, d, dx, dy, w, h, x, y) * @(#) * @(#) IMAGE *in, *out; * @(#) double a, b, c, d, dx, dy; * @(#) int w, h, x, y; * @(#) * @(#) Forward transform * @(#) X = a * x + b * y + dx * @(#) Y = c * x + d * y + dy * @(#) * @(#) x and y are the coordinates in input image. * @(#) X and Y are the coordinates in output image. * @(#) (0,0) is the upper left corner. * * Copyright N. Dessipris * Written on: 01/11/1991 * Modified on: 12/3/92 JC * - rounding error in interpolation routine fixed * - test for scale=1, angle=0 case fixed * - clipping of output removed: redundant * - various little tidies * - problems remain with scale>20, size<10 * * Re-written on: 20/08/92, J.Ph Laurent * * 21/02/93, JC * - speed-ups * - simplifications * - im_similarity now calculates a window and calls this routine * 6/7/93 JC * - rewritten for partials * - ANSIfied * - now rotates any non-complex type * 3/6/94 JC * - C revised in bug search * 9/6/94 JC * - im_prepare() was preparing too small an area! oops * 22/5/95 JC * - added code to detect all-black output area case - helps lazy ip * 3/7/95 JC * - IM_CODING_LABQ handling moved to here * 31/7/97 JC * - dx/dy sign reversed to be less confusing ... now follows comment at * top ... ax - by + dx etc. * - tiny speed up, replaced the *++ on interpolation with [z] * - im_similarity() moved in here * - args swapped: was whxy, now xywh * - didn't agree with dispatch fns before :( * 3/3/98 JC * - im_demand_hint() added * 20/12/99 JC * - im_affine() made from im_similarity_area() * - transform stuff cleaned up a bit * 14/4/01 JC * - oops, invert_point() had a rounding problem * 23/2/02 JC * - pre-calculate interpolation matricies * - integer interpolation for int8/16 types, double for * int32/float/double * - faster transformation * 15/8/02 JC * - records Xoffset/Yoffset * 14/4/04 * - rounding, clipping and transforming revised, now pixel-perfect (or * better than gimp, anyway) * 22/6/05 * - all revised again, simpler and more reliable now * 30/3/06 * - gah, still an occasional clipping problem * 12/7/06 * - still more tweaking, gah again * 7/10/06 * - set THINSTRIP for no-rotate affines */ /* 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 */ /* #define DEBUG #define DEBUG_GEOMETRY */ #ifdef HAVE_CONFIG_H #include #endif /*HAVE_CONFIG_H*/ #include #include #include #include #include #include #include #include "merge.h" #ifdef WITH_DMALLOC #include #endif /*WITH_DMALLOC*/ /* "fast" floor() ... on my laptop, anyway. */ #define FLOOR( V ) ((V) >= 0 ? (int)(V) : (int)((V) - 1)) /* Precalculate a whole bunch of interpolation matricies. int (used for pel * sizes up to short), and double (for all others). We go to scale + 1, so * we can round-to-nearest safely. FIXME ... should use seperable tables really */ static int im_affine_linear_int [TRANSFORM_SCALE + 1][TRANSFORM_SCALE + 1][4]; static double im_affine_linear_double [TRANSFORM_SCALE + 1][TRANSFORM_SCALE + 1][4]; /* Make sure the interpolation tables are built. */ static void affine_interpol_calc( void ) { static int calced = 0; int x, y; if( calced ) return; for( x = 0; x < TRANSFORM_SCALE + 1; x++ ) for( y = 0; y < TRANSFORM_SCALE + 1; y++ ) { double X, Y, Xd, Yd; double c1, c2, c3, c4; /* Interpolation errors. */ X = (double) x / TRANSFORM_SCALE; Y = (double) y / TRANSFORM_SCALE; Xd = 1.0 - X; Yd = 1.0 - Y; /* Weights. */ c1 = Xd*Yd; c2 = X*Yd; c3 = X*Y; c4 = Xd*Y; im_affine_linear_double[x][y][0] = c1; im_affine_linear_double[x][y][1] = c2; im_affine_linear_double[x][y][2] = c3; im_affine_linear_double[x][y][3] = c4; im_affine_linear_int[x][y][0] = c1 * INTERPOL_SCALE; im_affine_linear_int[x][y][1] = c2 * INTERPOL_SCALE; im_affine_linear_int[x][y][2] = c3 * INTERPOL_SCALE; im_affine_linear_int[x][y][3] = c4 * INTERPOL_SCALE; } calced = 1; } /* Calculate the inverse transformation. */ int im__transform_calc_inverse( Transformation *trn ) { DOUBLEMASK *msk, *msk2; if( !(msk = im_create_dmaskv( "boink", 2, 2, trn->a, trn->b, trn->c, trn->d )) ) return( -1 ); if( !(msk2 = im_matinv( msk, "boink2" )) ) { (void) im_free_dmask( msk ); return( -1 ); } trn->ia = msk2->coeff[0]; trn->ib = msk2->coeff[1]; trn->ic = msk2->coeff[2]; trn->id = msk2->coeff[3]; (void) im_free_dmask( msk ); (void) im_free_dmask( msk2 ); return( 0 ); } /* Init a Transform. */ void im__transform_init( Transformation *trn ) { trn->oarea.left = 0; trn->oarea.top = 0; trn->oarea.width = -1; trn->oarea.height = -1; trn->iarea.left = 0; trn->iarea.top = 0; trn->iarea.width = -1; trn->iarea.height = -1; trn->a = 1.0; /* Identity transform */ trn->b = 0.0; trn->c = 0.0; trn->d = 1.0; trn->dx = 0.0; trn->dy = 0.0; (void) im__transform_calc_inverse( trn ); } /* Test for transform is identity function. */ int im__transform_isidentity( Transformation *trn ) { if( trn->a == 1.0 && trn->b == 0.0 && trn->c == 0.0 && trn->d == 1.0 && trn->dx == 0.0 && trn->dy == 0.0 ) return( 1 ); else return( 0 ); } /* Map a pixel coordinate through the transform. */ void im__transform_forward( Transformation *trn, double x, double y, /* In input space */ double *ox, double *oy ) /* In output space */ { *ox = trn->a * x + trn->b * y + trn->dx; *oy = trn->c * x + trn->d * y + trn->dy; } /* Map a pixel coordinate through the inverse transform. */ void im__transform_inverse( Transformation *trn, double x, double y, /* In output space */ double *ox, double *oy ) /* In input space */ { double mx = x - trn->dx; double my = y - trn->dy; *ox = trn->ia * mx + trn->ib * my; *oy = trn->ic * mx + trn->id * my; } /* Combine two transformations. out can be one of the ins. */ int im__transform_add( Transformation *in1, Transformation *in2, Transformation *out ) { out->a = in1->a * in2->a + in1->c * in2->b; out->b = in1->b * in2->a + in1->d * in2->b; out->c = in1->a * in2->c + in1->c * in2->d; out->d = in1->b * in2->c + in1->d * in2->d; out->dx = in1->dx * in2->a + in1->dy * in2->b + in2->dx; out->dy = in1->dx * in2->c + in1->dy * in2->d + in2->dy; if( im__transform_calc_inverse( out ) ) return( -1 ); return( 0 ); } void im__transform_print( Transformation *trn ) { printf( "im__transform_print:\n" ); printf( " iarea: left=%d, top=%d, width=%d, height=%d\n", trn->iarea.left, trn->iarea.top, trn->iarea.width, trn->iarea.height ); printf( " oarea: left=%d, top=%d, width=%d, height=%d\n", trn->oarea.left, trn->oarea.top, trn->oarea.width, trn->oarea.height ); printf( " mat: a=%g, b=%g, c=%g, d=%g\n", trn->a, trn->b, trn->c, trn->d ); printf( " off: dx=%g, dy=%g\n", trn->dx, trn->dy ); } /* Map a point through the inverse transform. Used for clipping calculations, * so it takes account of iarea and oarea. */ static void invert_point( Transformation *trn, double x, double y, /* In output space */ double *ox, double *oy ) /* In input space */ { double xin = x - trn->oarea.left - trn->dx; double yin = y - trn->oarea.top - trn->dy; /* Find the inverse transform of current (x, y) */ *ox = trn->ia * xin + trn->ib * yin; *oy = trn->ic * xin + trn->id * yin; } /* Given a bounding box for an area in the output image, set the bounding box * for the corresponding pixels in the input image. */ static void invert_rect( Transformation *trn, Rect *in, /* In output space */ Rect *out ) /* In input space */ { double x1, y1; /* Map corners */ double x2, y2; double x3, y3; double x4, y4; double left, right, top, bottom; /* Map input Rect. */ invert_point( trn, in->left, in->top, &x1, &y1 ); invert_point( trn, in->left, IM_RECT_BOTTOM(in), &x2, &y2 ); invert_point( trn, IM_RECT_RIGHT(in), in->top, &x3, &y3 ); invert_point( trn, IM_RECT_RIGHT(in), IM_RECT_BOTTOM(in), &x4, &y4 ); /* Find bounding box for these four corners. */ left = IM_MIN( x1, IM_MIN( x2, IM_MIN( x3, x4 ) ) ); right = IM_MAX( x1, IM_MAX( x2, IM_MAX( x3, x4 ) ) ); top = IM_MIN( y1, IM_MIN( y2, IM_MIN( y3, y4 ) ) ); bottom = IM_MAX( y1, IM_MAX( y2, IM_MAX( y3, y4 ) ) ); /* Set output Rect. */ out->left = left; out->top = top; out->width = right - left + 1; out->height = bottom - top + 1; /* Add a border for interpolation. You'd think +1 would do it, but * we need to allow for rounding clipping as well. FIXME ... will need adjusting when we add bicubic */ im_rect_marginadjust( out, 2 ); } /* Interpolate a section ... int8/16 types. */ #define DO_IPEL(TYPE) { \ TYPE *tq = (TYPE *) q; \ \ int c1 = im_affine_linear_int[xi][yi][0]; \ int c2 = im_affine_linear_int[xi][yi][1]; \ int c3 = im_affine_linear_int[xi][yi][2]; \ int c4 = im_affine_linear_int[xi][yi][3]; \ \ /* p1 points to location (x_int, y_int) \ * p2 " " " (x_int+1, y_int) \ * p4 " " " (x_int+1, y_int+1) \ * p3 " " " (x_int, y_int+1) \ */ \ PEL *p1 = (PEL *) IM_REGION_ADDR( ir, x_int, y_int ); \ PEL *p2 = p1 + ofs2; \ PEL *p3 = p1 + ofs3; \ PEL *p4 = p1 + ofs4; \ TYPE *tp1 = (TYPE *) p1; \ TYPE *tp2 = (TYPE *) p2; \ TYPE *tp3 = (TYPE *) p3; \ TYPE *tp4 = (TYPE *) p4; \ \ /* Interpolate each band. \ */ \ for( z = 0; z < in->Bands; z++ ) \ tq[z] = (c1*tp1[z] + c2*tp2[z] + \ c3*tp3[z] + c4*tp4[z]) >> INTERPOL_SHIFT; \ } /* Interpolate a pel ... int32 and float types. */ #define DO_FPEL(TYPE) { \ TYPE *tq = (TYPE *) q; \ \ double c1 = im_affine_linear_double[xi][yi][0]; \ double c2 = im_affine_linear_double[xi][yi][1]; \ double c3 = im_affine_linear_double[xi][yi][2]; \ double c4 = im_affine_linear_double[xi][yi][3]; \ \ /* p1 points to location (x_int, y_int) \ * p2 " " " (x_int+1, y_int) \ * p4 " " " (x_int+1, y_int+1) \ * p3 " " " (x_int, y_int+1) \ */ \ PEL *p1 = (PEL *) IM_REGION_ADDR( ir, x_int, y_int ); \ PEL *p2 = p1 + ofs2; \ PEL *p3 = p1 + ofs3; \ PEL *p4 = p1 + ofs4; \ TYPE *tp1 = (TYPE *) p1; \ TYPE *tp2 = (TYPE *) p2; \ TYPE *tp3 = (TYPE *) p3; \ TYPE *tp4 = (TYPE *) p4; \ \ /* Interpolate each band. \ */ \ for( z = 0; z < in->Bands; z++ ) \ tq[z] = c1*tp1[z] + c2*tp2[z] + \ c3*tp3[z] + c4*tp4[z]; \ } static int affine_gen( REGION *or, void *seq, void *a, void *b ) { REGION *ir = (REGION *) seq; IMAGE *in = (IMAGE *) a; Transformation *trn = (Transformation *) b; /* Output area for this call. */ Rect *r = &or->valid; int le = r->left; int ri = IM_RECT_RIGHT(r); int to = r->top; int bo = IM_RECT_BOTTOM(r); Rect *iarea = &trn->iarea; Rect *oarea = &trn->oarea; int ps = IM_IMAGE_SIZEOF_PEL( in ); int x, y, z; /* Interpolation variables. */ int ofs2, ofs3, ofs4; /* Clipping Rects. */ Rect image, need, clipped; /* Find the area of the input image we need. */ image.left = 0; image.top = 0; image.width = in->Xsize; image.height = in->Ysize; invert_rect( trn, r, &need ); im_rect_intersectrect( &need, &image, &clipped ); /* Outside input image? All black. */ if( im_rect_isempty( &clipped ) ) { im__black_region( or ); return( 0 ); } /* We do need some pixels from the input image to make our output - * ask for them. */ if( im_prepare( ir, &clipped ) ) return( -1 ); #ifdef DEBUG printf( "affine: preparing left=%d, top=%d, width=%d, height=%d\n", clipped.left, clipped.top, clipped.width, clipped.height ); #endif /*DEBUG*/ /* Calculate pel offsets. */ ofs2 = IM_IMAGE_SIZEOF_PEL( in ); ofs3 = ofs2 + IM_REGION_LSKIP( ir ); ofs4 = IM_REGION_LSKIP( ir ); /* Resample! */ for( y = to; y < bo; y++ ) { /* Continuous cods in output space. */ double oy = y - oarea->top - trn->dy; double ox; /* Input clipping rectangle. */ int ile = iarea->left; int ito = iarea->top; int iri = iarea->left + iarea->width; int ibo = iarea->top + iarea->height; /* Derivative of matrix. */ double dx = trn->ia; double dy = trn->ic; /* Continuous cods in input space. */ double ix, iy; PEL *q; q = (PEL *) IM_REGION_ADDR( or, le, y ); ox = le - oarea->left - trn->dx; ix = trn->ia * ox + trn->ib * oy; iy = trn->ic * ox + trn->id * oy; /* Offset ix/iy input by iarea.left/top ... so we skip the * image edges we added for interpolation. */ ix += iarea->left; iy += iarea->top; for( x = le; x < ri; x++ ) { int fx, fy; fx = FLOOR( ix ); fy = FLOOR( iy ); /* Clipping! Use >= for right/bottom, since IPOL needs * to see one pixel more each way. */ if( fx < ile || fx >= iri || fy < ito || fy >= ibo ) { for( z = 0; z < ps; z++ ) q[z] = 0; } else { double sx, sy; int x_int, y_int; int xi, yi; /* Subtract 0.5 to centre the bilinear. FIXME ... need to adjust for bicubic. */ sx = ix - 0.5; sy = iy - 0.5; /* Now go to scaled int. */ sx *= TRANSFORM_SCALE; sy *= TRANSFORM_SCALE; x_int = FLOOR( sx ); y_int = FLOOR( sy ); /* Get index into interpolation table and * unscaled integer position. */ xi = x_int & (TRANSFORM_SCALE - 1); yi = y_int & (TRANSFORM_SCALE - 1); x_int = x_int >> TRANSFORM_SHIFT; y_int = y_int >> TRANSFORM_SHIFT; /* Interpolate for each input type. */ switch( in->BandFmt ) { case IM_BANDFMT_UCHAR: DO_IPEL( unsigned char ); break; case IM_BANDFMT_CHAR: DO_IPEL( char ); break; case IM_BANDFMT_USHORT: DO_IPEL( unsigned short ); break; case IM_BANDFMT_SHORT: DO_IPEL( short ); break; case IM_BANDFMT_UINT: DO_FPEL( unsigned int ); break; case IM_BANDFMT_INT: DO_FPEL( int ); break; case IM_BANDFMT_FLOAT: DO_FPEL( float ); break; case IM_BANDFMT_DOUBLE: DO_FPEL( double ); break; default: error_exit( "im_affine: panic!"); /*NOTREACHED*/ } } ix += dx; iy += dy; q += ps; } } return( 0 ); } static int affine( IMAGE *in, IMAGE *out, Transformation *trn ) { Transformation *trn2; double edge; if( im_iscomplex( in ) ) { im_errormsg( "im_affine: complex input not supported" ); return( -1 ); } /* We output at (0,0), so displace output by that amount -ve to get * output at (ox,oy). Alter our copy of trn. */ if( !(trn2 = IM_NEW( out, Transformation )) ) return( -1 ); *trn2 = *trn; trn2->oarea.left = -trn->oarea.left; trn2->oarea.top = -trn->oarea.top; if( im__transform_calc_inverse( trn2 ) ) return( -1 ); /* Make output image. */ if( im_piocheck( in, out ) ) return( -1 ); if( im_cp_desc( out, in ) ) return( -1 ); out->Xsize = trn2->oarea.width; out->Ysize = trn2->oarea.height; /* Normally SMALLTILE ... except if this is a size up/down affine. */ if( trn->b == 0.0 && trn->c == 0.0 ) { if( im_demand_hint( out, IM_FATSTRIP, in, NULL ) ) return( -1 ); } else { if( im_demand_hint( out, IM_SMALLTILE, in, NULL ) ) return( -1 ); } /* Check for coordinate overflow ... we want to be able to hold the * output space inside INT_MAX / TRANSFORM_SCALE. */ edge = INT_MAX / TRANSFORM_SCALE; if( trn2->oarea.left < -edge || trn2->oarea.top < -edge || IM_RECT_RIGHT( &trn2->oarea ) > edge || IM_RECT_BOTTOM( &trn2->oarea ) > edge ) { im_errormsg( "im_affine: output coordinates out of range" ); return( -1 ); } /* Generate! */ if( im_generate( out, im_start_one, affine_gen, im_stop_one, in, trn2 ) ) return( -1 ); return( 0 ); } /* As above, but do IM_CODING_LABQ too. And embed the input. */ int im__affine( IMAGE *in, IMAGE *out, Transformation *trn ) { IMAGE *t3 = im_open_local( out, "im_affine:3", "p" ); Transformation trn2; #ifdef DEBUG_GEOMETRY printf( "im__affine: %s\n", in->filename ); im__transform_print( trn ); #endif /*DEBUG_GEOMETRY*/ /* Add new pixels around the input so we can interpolate at the edges. * Bilinear needs 0.5 pixels on all edges. FIXME ... will need to fiddle with this when we add bicubic */ if( !t3 || im_embed( in, t3, 1, 1, 1, in->Xsize + 2, in->Ysize + 2 ) ) return( -1 ); /* Set iarea so we know what part of the input we can take. */ trn2 = *trn; trn2.iarea.left += 1; trn2.iarea.top += 1; affine_interpol_calc(); if( in->Coding == IM_CODING_LABQ ) { IMAGE *t1 = im_open_local( out, "im_affine:1", "p" ); IMAGE *t2 = im_open_local( out, "im_affine:2", "p" ); if( !t1 || !t2 || im_LabQ2LabS( t3, t1 ) || affine( t1, t2, &trn2 ) || im_LabS2LabQ( t2, out ) ) return( -1 ); } else if( in->Coding == IM_CODING_NONE ) { if( affine( t3, out, &trn2 ) ) return( -1 ); } else { im_errormsg( "im_affine: unknown coding type" ); return( -1 ); } /* Finally: can now set Xoffset/Yoffset. */ out->Xoffset = trn->dx - trn->oarea.left; out->Yoffset = trn->dy - trn->oarea.top; return( 0 ); } int im_affine( IMAGE *in, IMAGE *out, double a, double b, double c, double d, double dx, double dy, int ox, int oy, int ow, int oh ) { Transformation trn; trn.oarea.left = ox; trn.oarea.top = oy; trn.oarea.width = ow; trn.oarea.height = oh; trn.iarea.left = 0; trn.iarea.top = 0; trn.iarea.width = in->Xsize; trn.iarea.height = in->Ysize; trn.a = a; trn.b = b; trn.c = c; trn.d = d; trn.dx = dx; trn.dy = dy; return( im__affine( in, out, &trn ) ); }