# include # include # include # include # include "polygon_integrals.h" /******************************************************************************/ int i4_max ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MAX returns the maximum of two I4's. Licensing: This code is distributed under the MIT license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, are two integers to be compared. Output, int I4_MAX, the larger of I1 and I2. */ { int value; if ( i2 < i1 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ int i4_min ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MIN returns the smaller of two I4's. Licensing: This code is distributed under the MIT license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, two integers to be compared. Output, int I4_MIN, the smaller of I1 and I2. */ { int value; if ( i1 < i2 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ double moment ( int n, double x[], double y[], int p, int q ) /******************************************************************************/ /* Purpose: MOMENT computes an unnormalized moment of a polygon. Discussion: Nu(P,Q) = Integral ( x, y in polygon ) x^p y^q dx dy Licensing: This code is distributed under the MIT license. Modified: 03 October 2012 Author: John Burkardt Reference: Carsten Steger, On the calculation of arbitrary moments of polygons, Technical Report FGBV-96-05, Forschungsgruppe Bildverstehen, Informatik IX, Technische Universitaet Muenchen, October 1996. Parameters: Input, int N, the number of vertices of the polygon. Input, double X[N], Y[N], the vertex coordinates. Input, int P, Q, the indices of the moment. Output, double MOMENT, the unnormalized moment Nu(P,Q). */ { int i; int k; int l; double nu_pq; double s_pq; double xi; double xj; double yi; double yj; nu_pq = 0.0; xj = x[n-1]; yj = y[n-1]; for ( i = 0; i < n; i++ ) { xi = x[i]; yi = y[i]; s_pq = 0.0; for ( k = 0; k <= p; k++ ) { for ( l = 0; l <= q; l++ ) { s_pq = s_pq + r8_choose ( k + l, l ) * r8_choose ( p + q - k - l, q - l ) * pow ( xi, k ) * pow ( xj, p - k ) * pow ( yi, l ) * pow ( yj, q - l ); } } nu_pq = nu_pq + ( xj * yi - xi * yj ) * s_pq; xj = xi; yj = yi; } nu_pq = nu_pq / ( double ) ( p + q + 2 ) / ( double ) ( p + q + 1 ) / r8_choose ( p + q, p ); return nu_pq; } /******************************************************************************/ double moment_central ( int n, double x[], double y[], int p, int q ) /******************************************************************************/ /* Purpose: MOMENT_CENTRAL computes central moments of a polygon. Discussion: The central moment Mu(P,Q) is defined by Mu(P,Q) = Integral ( polygon ) (x-Alpha(1,0))^p (y-Alpha(0,1))^q dx dy / Area ( polygon ) where Alpha(1,0) = Integral ( polygon ) x dx dy / Area ( polygon ) Alpha(0,1) = Integral ( polygon ) y dx dy / Area ( polygon ) Licensing: This code is distributed under the MIT license. Modified: 03 October 2012 Author: John Burkardt Reference: Carsten Steger, On the calculation of arbitrary moments of polygons, Technical Report FGBV-96-05, Forschungsgruppe Bildverstehen, Informatik IX, Technische Universitaet Muenchen, October 1996. Parameters: Input, int N, the number of vertices of the polygon. Input, double X[N], Y[N], the vertex coordinates. Input, int P, Q, the indices of the moment. Output, double MOMENT_CENTRAL, the unnormalized moment Mu(P,Q). */ { double alpha_01; double alpha_10; double alpha_ij; int i; int j; double mu_pq; alpha_10 = moment_normalized ( n, x, y, 1, 0 ); alpha_01 = moment_normalized ( n, x, y, 0, 1 ); mu_pq = 0.0; for ( i = 0; i <= p; i++ ) { for ( j = 0; j <= q; j++ ) { alpha_ij = moment_normalized ( n, x, y, i, j ); mu_pq = mu_pq + r8_mop ( p + q - i - j ) * r8_choose ( p, i ) * r8_choose ( q, j ) * pow ( alpha_10, p - i ) * pow ( alpha_01, q - j ) * alpha_ij; } } return mu_pq; } /******************************************************************************/ double moment_normalized ( int n, double x[], double y[], int p, int q ) /******************************************************************************/ /* Purpose: MOMENT_NORMALIZED computes a normalized moment of a polygon. Discussion: Alpha(P,Q) = Integral ( x, y in polygon ) x^p y^q dx dy / Area ( polygon ) Licensing: This code is distributed under the MIT license. Modified: 03 October 2012 Author: John Burkardt Reference: Carsten Steger, On the calculation of arbitrary moments of polygons, Technical Report FGBV-96-05, Forschungsgruppe Bildverstehen, Informatik IX, Technische Universitaet Muenchen, October 1996. Parameters: Input, int N, the number of vertices of the polygon. Input, double X[N], Y[N], the vertex coordinates. Input, int P, Q, the indices of the moment. Output, double MOMENT_NORMALIZED, the normalized moment Alpha(P,Q). */ { double alpha_pq; double nu_00; double nu_pq; nu_pq = moment ( n, x, y, p, q ); nu_00 = moment ( n, x, y, 0, 0 ); alpha_pq = nu_pq / nu_00; return alpha_pq; } /******************************************************************************/ double r8_choose ( int n, int k ) /******************************************************************************/ /* Purpose: R8_CHOOSE computes the binomial coefficient C(N,K) as an R8. Discussion: The value is calculated in such a way as to avoid overflow and roundoff. The calculation is done in R8 arithmetic. The formula used is: C(N,K) = N / ( K * (N-K) ) Licensing: This code is distributed under the MIT license. Modified: 01 July 2008 Author: John Burkardt Reference: ML Wolfson, HV Wright, Algorithm 160: Combinatorial of M Things Taken N at a Time, Communications of the ACM, Volume 6, Number 4, April 1963, page 161. Parameters: Input, int N, K, the values of N and K. Output, double R8_CHOOSE, the number of combinations of N things taken K at a time. */ { int i; int mn; int mx; double value; mn = i4_min ( k, n - k ); if ( mn < 0 ) { value = 0.0; } else if ( mn == 0 ) { value = 1.0; } else { mx = i4_max ( k, n - k ); value = ( double ) ( mx + 1 ); for ( i = 2; i <= mn; i++ ) { value = ( value * ( double ) ( mx + i ) ) / ( double ) i; } } return value; } /******************************************************************************/ double r8_mop ( int i ) /******************************************************************************/ /* Purpose: R8_MOP returns the I-th power of -1 as an R8 value. Discussion: An R8 is an double value. Licensing: This code is distributed under the MIT license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input, int I, the power of -1. Output, double R8_MOP, the I-th power of -1. */ { double value; if ( ( i % 2 ) == 0 ) { value = + 1.0; } else { value = - 1.0; } return value; } /******************************************************************************/ void timestamp ( void ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }