/* im_linreg.c * * Copyright: 2006, The Nottingham Trent University * * Author: Tom Vajzovic * * Written on: 2006-12-26 * * 12/5/09 * - make x_anal() static, fix some signed/unsigned warnings * 3/9/09 * - gtkdoc comment */ /* 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 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 #include #include typedef struct { unsigned int n; double *xs; double *difs; double mean; double nsig2; double err_term; } x_set; #define LINREG_SEQ( TYPE ) typedef struct { \ REGION **regs; \ TYPE **ptrs; \ size_t *skips; \ } linreg_seq_ ## TYPE LINREG_SEQ( gint8 ); LINREG_SEQ( guint8 ); LINREG_SEQ( gint16 ); LINREG_SEQ( guint16 ); LINREG_SEQ( gint32 ); LINREG_SEQ( guint32 ); LINREG_SEQ( float ); LINREG_SEQ( double ); static x_set *x_anal( IMAGE *im, double *xs, unsigned int n ); #define LINREG_START_DECL( TYPE ) static void * linreg_start_ ## TYPE( IMAGE *, void *, void * ); #define LINREG_GEN_DECL( TYPE ) static int linreg_gen_ ## TYPE( REGION *, void *, void *, void * ); #define LINREG_STOP_DECL( TYPE ) static int linreg_stop_ ## TYPE( void *, void *, void * ); #define INCR_ALL_DECL( TYPE ) static void incr_all_ ## TYPE( TYPE **ptrs, unsigned int n ) #define SKIP_ALL_DECL( TYPE ) static void skip_all_ ## TYPE( TYPE **ptrs, size_t *skips, unsigned int n ) LINREG_START_DECL( gint8 ); LINREG_START_DECL( guint8 ); LINREG_START_DECL( gint16 ); LINREG_START_DECL( guint16 ); LINREG_START_DECL( gint32 ); LINREG_START_DECL( guint32 ); LINREG_START_DECL( float ); LINREG_START_DECL( double ); LINREG_GEN_DECL( gint8 ); LINREG_GEN_DECL( guint8 ); LINREG_GEN_DECL( gint16 ); LINREG_GEN_DECL( guint16 ); LINREG_GEN_DECL( gint32 ); LINREG_GEN_DECL( guint32 ); LINREG_GEN_DECL( float ); LINREG_GEN_DECL( double ); LINREG_STOP_DECL( gint8 ); LINREG_STOP_DECL( guint8 ); LINREG_STOP_DECL( gint16 ); LINREG_STOP_DECL( guint16 ); LINREG_STOP_DECL( gint32 ); LINREG_STOP_DECL( guint32 ); LINREG_STOP_DECL( float ); LINREG_STOP_DECL( double ); INCR_ALL_DECL( gint8 ); INCR_ALL_DECL( guint8 ); INCR_ALL_DECL( gint16 ); INCR_ALL_DECL( guint16 ); INCR_ALL_DECL( gint32 ); INCR_ALL_DECL( guint32 ); INCR_ALL_DECL( float ); INCR_ALL_DECL( double ); SKIP_ALL_DECL( gint8 ); SKIP_ALL_DECL( guint8 ); SKIP_ALL_DECL( gint16 ); SKIP_ALL_DECL( guint16 ); SKIP_ALL_DECL( gint32 ); SKIP_ALL_DECL( guint32 ); SKIP_ALL_DECL( float ); SKIP_ALL_DECL( double ); /** * im_linreg: * @ins: NULL-terminated array of input images * @out: results of analysis * @xs: X position of each image (pixel value is Y) * * Function to find perform pixelwise linear regression on an array of * single band images. The output is a seven-band douuble image * * TODO: figure out how this works and fix up these docs! */ int im_linreg( IMAGE **ins, IMAGE *out, double *xs ){ #define FUNCTION_NAME "im_linreg" int n; x_set *x_vals; if( im_poutcheck( out ) ) return( -1 ); for( n= 0; ins[ n ]; ++n ){ /* if( ! isfinite( xs[ n ] ) ){ im_error( FUNCTION_NAME, "invalid argument" ); return( -1 ); } */ if( im_pincheck( ins[ n ] ) ) return( -1 ); if( 1 != ins[ n ]-> Bands ){ im_error( FUNCTION_NAME, "image is not single band" ); return( -1 ); } if( ins[ n ]-> Coding ){ im_error( FUNCTION_NAME, "image is not uncoded" ); return( -1 ); } if( n ){ if( ins[ n ]-> BandFmt != ins[ 0 ]-> BandFmt ){ im_error( FUNCTION_NAME, "image band formats differ" ); return( -1 ); } } else { if( vips_band_format_iscomplex( ins[ 0 ]->BandFmt ) ){ im_error( FUNCTION_NAME, "image has non-scalar band format" ); return( -1 ); } } if( n && ( ins[ n ]-> Xsize != ins[ 0 ]-> Xsize || ins[ n ]-> Ysize != ins[ 0 ]-> Ysize ) ){ im_error( FUNCTION_NAME, "image sizes differ" ); return( -1 ); } } if( n < 3 ){ im_error( FUNCTION_NAME, "not enough input images" ); return( -1 ); } if( im_cp_desc_array( out, ins ) ) return( -1 ); out-> Bands= 7; out-> BandFmt= IM_BANDFMT_DOUBLE; out-> Type= 0; if( im_demand_hint_array( out, IM_THINSTRIP, ins ) ) return( -1 ); x_vals= x_anal( out, xs, n ); if( ! x_vals ) return( -1 ); switch( ins[ 0 ]-> BandFmt ){ #define LINREG_RET( TYPE ) return im_generate( out, linreg_start_ ## TYPE, linreg_gen_ ## TYPE, linreg_stop_ ## TYPE, ins, x_vals ) case IM_BANDFMT_CHAR: LINREG_RET( gint8 ); case IM_BANDFMT_UCHAR: LINREG_RET( guint8 ); case IM_BANDFMT_SHORT: LINREG_RET( gint16 ); case IM_BANDFMT_USHORT: LINREG_RET( guint16 ); case IM_BANDFMT_INT: LINREG_RET( gint32 ); case IM_BANDFMT_UINT: LINREG_RET( guint32 ); case IM_BANDFMT_FLOAT: LINREG_RET( float ); case IM_BANDFMT_DOUBLE: LINREG_RET( double ); default: /* keep -Wall happy */ return( -1 ); } #undef FUNCTION_NAME } static x_set *x_anal( IMAGE *im, double *xs, unsigned int n ){ unsigned int i; x_set *x_vals= IM_NEW( im, x_set ); if( ! x_vals ) return( NULL ); x_vals-> xs= IM_ARRAY( im, 2 * n, double ); if( ! x_vals-> xs ) return( NULL ); x_vals-> difs= x_vals-> xs + n; x_vals-> n= n; x_vals-> mean= 0.0; for( i= 0; i < n; ++i ){ x_vals-> xs[ i ]= xs[ i ]; x_vals-> mean+= xs[ i ]; } x_vals-> mean/= n; x_vals-> nsig2= 0.0; for( i= 0; i < n; ++i ){ x_vals-> difs[ i ]= xs[ i ] - x_vals-> mean; x_vals-> nsig2+= x_vals-> difs[ i ] * x_vals-> difs[ i ]; } x_vals-> err_term= ( 1.0 / (double) n ) + ( ( x_vals-> mean * x_vals-> mean ) / x_vals-> nsig2 ); return( x_vals ); } #define LINREG_START_DEFN( TYPE ) static void *linreg_start_ ## TYPE( IMAGE *out, void *a, void *b ){ \ IMAGE **ins= (IMAGE **) a; \ x_set *x_vals= (x_set *) b; \ linreg_seq_ ## TYPE *seq= IM_NEW( out, linreg_seq_ ## TYPE ); \ \ if( ! seq ) \ return NULL; \ \ seq-> regs= im_start_many( NULL, ins, NULL ); \ seq-> ptrs= IM_ARRAY( out, x_vals-> n, TYPE* ); \ seq-> skips= IM_ARRAY( out, x_vals-> n, size_t ); \ \ if( ! seq-> ptrs || ! seq-> regs || ! seq-> skips ){ \ linreg_stop_ ## TYPE( seq, NULL, NULL ); \ return NULL; \ } \ return (void *) seq; \ } #define N ( (double) n ) #define y(a) ( (double) (* seq-> ptrs[(a)] ) ) #define x(a) ( (double) ( x_vals-> xs[(a)] ) ) #define xd(a) ( (double) ( x_vals-> difs[(a)] ) ) #define Sxd2 ( x_vals-> nsig2 ) #define mean_x ( x_vals-> mean ) #define mean_y ( out[0] ) #define dev_y ( out[1] ) #define y_x0 ( out[2] ) #define d_y_x0 ( out[3] ) #define dy_dx ( out[4] ) #define d_dy_dx ( out[5] ) #define R ( out[6] ) #define LINREG_GEN_DEFN( TYPE ) static int linreg_gen_ ## TYPE( REGION *to_make, void *vseq, void *unrequired, void *b ){ \ linreg_seq_ ## TYPE *seq= (linreg_seq_ ## TYPE *) vseq; \ x_set *x_vals= (x_set *) b; \ unsigned int n= x_vals-> n; \ double *out= (double*) IM_REGION_ADDR_TOPLEFT( to_make ); \ size_t out_skip= IM_REGION_LSKIP( to_make ) / sizeof( double ); \ double *out_end= out + out_skip * to_make-> valid. height; \ double *out_stop; \ size_t out_n= IM_REGION_N_ELEMENTS( to_make ); \ unsigned int i; \ \ out_skip-= out_n; \ \ if( im_prepare_many( seq-> regs, & to_make-> valid ) ) \ return -1; \ \ for( i= 0; i < n; ++i ){ \ seq-> ptrs[ i ]= (TYPE*) IM_REGION_ADDR( seq-> regs[ i ], to_make-> valid. left, to_make-> valid. top ); \ seq-> skips[ i ]= ( IM_REGION_LSKIP( seq-> regs[ i ] ) / sizeof( TYPE ) ) - IM_REGION_N_ELEMENTS( seq-> regs[ i ] ); \ } \ \ for( ; out < out_end; out+= out_skip, skip_all_ ## TYPE( seq-> ptrs, seq-> skips, n ) ) \ for( out_stop= out + out_n; out < out_stop; out+= 7, incr_all_ ## TYPE( seq-> ptrs, n ) ){ \ double Sy= 0.0; \ double Sxd_y= 0.0; \ double Syd2= 0.0; \ double Sxd_yd= 0.0; \ double Se2= 0.0; \ \ for( i= 0; i < n; ++i ){ \ Sy+= y(i); \ Sxd_y+= xd(i) * y(i); \ } \ mean_y= Sy / N; \ dy_dx= Sxd_y / Sxd2; \ y_x0= mean_y - dy_dx * mean_x; \ \ for( i= 0; i < n; ++i ){ \ double yd= y(i) - mean_y; \ double e= y(i) - dy_dx * x(i) - y_x0; \ Syd2+= yd * yd; \ Sxd_yd+= xd(i) * yd; \ Se2+= e * e; \ } \ dev_y= sqrt( Syd2 / N ); \ Se2/= ( N - 2.0 ); \ d_dy_dx= sqrt( Se2 / Sxd2 ); \ d_y_x0= sqrt( Se2 * x_vals-> err_term ); \ R= Sxd_yd / sqrt( Sxd2 * Syd2 ); \ } \ return 0; \ } #define LINREG_STOP_DEFN( TYPE ) static int linreg_stop_ ## TYPE( void *vseq, void *a, void *b ){ \ linreg_seq_ ## TYPE *seq = (linreg_seq_ ## TYPE *) vseq; \ if( seq-> regs ) \ im_stop_many( seq-> regs, NULL, NULL ); \ return 0; \ } #define INCR_ALL_DEFN( TYPE ) static void incr_all_ ## TYPE( TYPE **ptrs, unsigned int n ){ \ TYPE **stop= ptrs + n; \ for( ; ptrs < stop; ++ptrs ) \ ++*ptrs; \ } #define SKIP_ALL_DEFN( TYPE ) static void skip_all_ ## TYPE( TYPE **ptrs, size_t *skips, unsigned int n ){ \ TYPE **stop= ptrs + n; \ for( ; ptrs < stop; ++ptrs, ++skips ) \ *ptrs+= *skips; \ } LINREG_START_DEFN( gint8 ); LINREG_START_DEFN( guint8 ); LINREG_START_DEFN( gint16 ); LINREG_START_DEFN( guint16 ); LINREG_START_DEFN( gint32 ); LINREG_START_DEFN( guint32 ); LINREG_START_DEFN( float ); LINREG_START_DEFN( double ); LINREG_GEN_DEFN( gint8 ); LINREG_GEN_DEFN( guint8 ); LINREG_GEN_DEFN( gint16 ); LINREG_GEN_DEFN( guint16 ); LINREG_GEN_DEFN( gint32 ); LINREG_GEN_DEFN( guint32 ); LINREG_GEN_DEFN( float ); LINREG_GEN_DEFN( double ); LINREG_STOP_DEFN( gint8 ); LINREG_STOP_DEFN( guint8 ); LINREG_STOP_DEFN( gint16 ); LINREG_STOP_DEFN( guint16 ); LINREG_STOP_DEFN( gint32 ); LINREG_STOP_DEFN( guint32 ); LINREG_STOP_DEFN( float ); LINREG_STOP_DEFN( double ); INCR_ALL_DEFN( gint8 ); INCR_ALL_DEFN( guint8 ); INCR_ALL_DEFN( gint16 ); INCR_ALL_DEFN( guint16 ); INCR_ALL_DEFN( gint32 ); INCR_ALL_DEFN( guint32 ); INCR_ALL_DEFN( float ); INCR_ALL_DEFN( double ); SKIP_ALL_DEFN( gint8 ); SKIP_ALL_DEFN( guint8 ); SKIP_ALL_DEFN( gint16 ); SKIP_ALL_DEFN( guint16 ); SKIP_ALL_DEFN( gint32 ); SKIP_ALL_DEFN( guint32 ); SKIP_ALL_DEFN( float ); SKIP_ALL_DEFN( double );