libvips/libsrc/mosaicing/im_clinear.c

184 lines
4.5 KiB
C
Raw Normal View History

2007-08-29 18:23:50 +02:00
/* @(#) Function which calculates the coefficients between corresponding
* @(#) points from reference and secondary images (probably from the scanner),
* @(#) previously calculated using the functions im_calcon() and im_chpair()
* @(#) It is assummed that a selection of the best(?) possible points has
* @(#) been already carried out and that those nopoints points are in arrays
* @(#) x1, y1 and x2, y2
* @(#) No IMAGES are involved in this function and the calculated parameters
* @(#) are returned in scale angle deltax and deltay of the TIE_POINTS struct.
* @(#)
* @(#) int im_clinear( points )
* @(#) TIE_POINTS *points;
* @(#)
* @(#) Returns 0 on sucess and -1 on error.
*
* Copyright: 1990, N. Dessipris.
*
* Author: Nicos Dessipris
* Written on: 20/12/1990
* Modified on : 18/04/1991
* 24/1/97 JC
* - tiny mem leak fixed
*/
/*
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 <config.h>
#endif /*HAVE_CONFIG_H*/
#include <vips/intl.h>
#include <stdio.h>
#include <math.h>
#include <vips/vips.h>
#include "mosaic.h"
#ifdef WITH_DMALLOC
#include <dmalloc.h>
#endif /*WITH_DMALLOC*/
int
im__clinear( TIE_POINTS *points )
{
double **mat; /* matrix mar[4][4] */
double *g; /* vector g[1][4] */
double value;
double sx1=0.0, sx1x1=0.0, sy1=0.0, sy1y1=0.0, sx1y1 = 0.0;
double sx2x1=0.0, sx2y1=0.0, sx2=0.0, sy2=0.0, sy2y1=0.0, sy2x1=0.0;
int i, j;
int elms;
double scale, angle, xdelta, ydelta;
int *xref, *yref, *xsec, *ysec;
double *dx, *dy, *dev;
double resx, resy;
xref = &points->x_reference[0];
yref = &points->y_reference[0];
xsec = &points->x_secondary[0];
ysec = &points->y_secondary[0];
dx = &points->dx[0];
dy = &points->dy[0];
dev = &points->deviation[0];
elms = points->nopoints;
if( !(mat = im_dmat_alloc( 0, 3, 0, 3 )) )
return( -1 );
if( !(g = im_dvector( 0, 3 )) ) {
im_free_dmat( mat, 0, 3, 0, 3 );
return( -1 );
}
resx = 0.0;
resy = 0.0;
for( i = 0; i < points->nopoints; i++ ) {
sx1 += xref[i];
sx1x1 += xref[i] * xref[i];
sy1 += yref[i];
sy1y1 += yref[i] * yref[i];
sx1y1 += xref[i] * yref[i];
sx2x1 += xsec[i] * xref[i];
sx2y1 += xsec[i] * yref[i];
sy2y1 += ysec[i] * yref[i];
sy2x1 += ysec[i] * xref[i];
sx2 += xsec[i];
sy2 += ysec[i];
}
resx = fabs( sx1-sx2 )/points->nopoints;
resy = fabs( sy1-sy2 )/points->nopoints;
mat[0][0] = sx1x1 + sy1y1;
mat[0][1] = 0;
mat[0][2] = sx1;
mat[0][3] = sy1;
mat[1][0] = 0;
mat[1][1] = sx1x1 + sy1y1;
mat[1][2] = -sy1;
mat[1][3] = sx1;
mat[2][0] = sx1;
mat[2][1] = -sy1;
mat[2][2] = (double)elms;
mat[2][3] = 0.0;
mat[3][0] = sy1;
mat[3][1] = sx1;
mat[3][2] = 0.0;
mat[3][3] = (double)elms;
g[0] = sx2x1 + sy2y1;
g[1] = -sx2y1 + sy2x1;
g[2] = sx2;
g[3] = sy2;
if( im_invmat( mat, 4 ) ) {
im_free_dmat( mat, 0, 3, 0, 3 );
im_free_dvector( g, 0, 3 );
im_errormsg( "im_clinear: im_invmat failed" );
return( -1 );
}
scale = 0.0; angle = 0.0;
xdelta = 0.0; ydelta = 0.0;
for( j = 0; j < 4; j++ ) {
scale += mat[0][j] * g[j];
angle += mat[1][j] * g[j];
xdelta += mat[2][j] * g[j];
ydelta += mat[3][j] * g[j];
}
/* find the deviation of each point for the estimated variables
* if it greater than 1 then the solution is not good enough
* but this is handled by the main program
*/
for( i = 0; i < points->nopoints; i++ ) {
dx[i] = xsec[i] -
((scale * xref[i]) - (angle * yref[i]) + xdelta);
dy[i] = ysec[i] -
((angle * xref[i]) + (scale * yref[i]) + ydelta);
value = sqrt( dx[i]*dx[i] + dy[i]*dy[i] );
dev[i] = value;
}
points->l_scale = scale;
points->l_angle = angle;
points->l_deltax = xdelta;
points->l_deltay = ydelta;
im_free_dmat( mat, 0, 3, 0, 3 );
im_free_dvector( g, 0, 3 );
return( 0 );
}