libvips/libsrc/arithmetic/im_linreg.c

/* @(#) Function to find perform pixelwise linear regression on an array of 
 * @(#) single band images.
 * @(#)
 * @(#) int im_linreg(
 * @(#)   IMAGE **ins,
 * @(#)   IMAGE *out,
 * @(#)   double *xs
 * @(#) );
 * @(#) 
 *
 * 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
 */

/*

    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

 */


/** HEADERS **/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif /* HAVE_CONFIG_H */
#include <vips/intl.h>

#include <stdlib.h>
#include <math.h>
#include <vips/vips.h>

#ifdef WITH_DMALLOC
#include <dmalloc.h>
#endif /*WITH_DMALLOC */


/** TYPES **/

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 );


/** LOCAL FUNCTION DECLARATIONS **/

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 );


/** EXPORTED FUNCTION DEFINITION **/

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( ! im_isscalar( ins[ 0 ] ) ){
        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-> Bbits= IM_BBITS_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
}


/** LOCAL FUNCTION DECLARATIONS **/

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 );