📄 shepard.cpp
字号:
///////////////////////////////////////////////////////////
// //
// SAGA //
// //
// System for Automated Geoscientific Analyses //
// //
// Module Library: //
// Grid_Gridding //
// //
//-------------------------------------------------------//
// //
// Shepard.cpp //
// //
// Copyright (C) 2003 by //
// Andre Ringeler //
// //
//-------------------------------------------------------//
// //
// This file is part of 'SAGA - System for Automated //
// Geoscientific Analyses'. SAGA is free software; you //
// can redistribute it and/or modify it under the terms //
// of the GNU General Public License as published by the //
// Free Software Foundation; version 2 of the License. //
// //
// SAGA 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 General Public //
// License for more details. //
// //
// You should have received a copy of the GNU 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. //
// //
//-------------------------------------------------------//
// //
// e-mail: aringel@gwdg.de //
// //
// contact: Andre Ringeler //
// Institute of Geography //
// University of Goettingen //
// Goldschmidtstr. 5 //
// 37077 Goettingen //
// Germany //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
#include <math.h>
#include <stdlib.h>
#include "Shepard.h"
double _missing_;
CShepard2d::CShepard2d (void)
{
m_cells = NULL;
m_next = NULL;
m_rsq = NULL;
m_a = NULL;
_missing_ = -9999.99;
}
CShepard2d::~CShepard2d (void)
{
if (m_cells)
{
free (m_cells);
m_cells = NULL;
}
if (m_next)
{
free (m_next);
m_next = NULL;
}
if (m_rsq)
{
free (m_rsq);
m_rsq = NULL;
}
if (m_a)
{
free (m_a);
m_a = NULL;
}
}
void CShepard2d::Set_Missing(double missing)
{
_missing_ = missing;
}
int CShepard2d::Interpolate (double *X, double * Y, double * F, int N_Points, int Quadratic_Neighbors, int Weighting_Neighbors )
{
int nr, lmax;
int status;
if (N_Points < 6)
return -1;
lmax = min(40, N_Points - 1);
if (Quadratic_Neighbors < 5 || Quadratic_Neighbors > lmax)
return -1;
if (Weighting_Neighbors < 1 || Weighting_Neighbors > lmax)
return -1;
nr = (int)(sqrt(N_Points / 3.0));
if (nr < 1)
nr = 1;
this->CShepard2d::~CShepard2d();
m_cells = (int *) malloc(nr * nr * sizeof(int));
m_next = (int *) malloc(N_Points * sizeof(int));
m_rsq = (double *) malloc(N_Points * sizeof(double));
m_a = (double *) malloc(N_Points * 5 * sizeof(double));
m_x = X;
m_y = Y;
m_f = F;
m_nPoints = N_Points;
m_nr = nr;
qshep2_(&N_Points, X, Y, F, &Quadratic_Neighbors, &Weighting_Neighbors, &nr, m_cells, m_next, &xmin, &ymin, &dx, &dy, &rmax, m_rsq, m_a, &status);
return status;
}
void CShepard2d::GetValue(double px, double py, double &q)
{
if ( ! m_a )
q = _missing_;
else
q = qs2val_(&px, &py, &m_nPoints , m_x, m_y, m_f, &m_nr, m_cells, m_next, &xmin, &ymin, &dx, &dy, &rmax, m_rsq, m_a);
}
/*/////////////////////////////////////////////////////////////////////////////////////
The folowing source ist an automatic translation
with f2c of Module 660 in TOMS
QSHEP2D: Fortran routines implementing the quadratic Shepard method for
bivariate interpolation of scattered data. (See R. J. Renka, ACM TOMS 14 (1988)
pp. 149-150.).
Classes : E2b . Interpolation of scattered, non-gridded multivariate data
Type : Fortran software in TOMS collection.
Access : Some uses prohibited. Portable.
Precision: Single.
Details : Fullsource
Sites : (1) NETLIB
//////////////////////////////////////////////////////////////////////////////////////*/
struct
{
double y;
}
stcom_;
#define stcom_1 stcom_
/* Table of constant values */
static int c__1 = 1;
int qshep2_(int *n, double *x, double *y, double *f, int *
nq, int *nw, int *nr, int *lcell, int *lnext, double *
xmin, double *ymin, double *dx, double *dy, double *rmax, double *rsq, double *a,
int *ier)
{
static double rtol = 1e-5f;
static double dtol = .01f;
static double sf = 1.f;
int lcell_dim1, lcell_offset, i__1, i__2, i__3;
double r__1, r__2;
// double sqrt(double);
static double b[36] /* was[6][6] */, c__;
static int i__, j, k;
static double s, t;
static int ib;
static double fk, av;
static int nn, np;
static double rq, xk, rs, yk;
static int ip1, jp1;
static double ddx, ddy;
static int neq, lnp, nnq, nnr, nnw;
static double xmn, sum, ymn, rws;
static int irm1;
static double dmin__;
static int ierr, lmax;
static double avsq;
static int irow, npts[40];
static double rsmx, rsold;
extern int getnp2_(double *, double *, double *, double *,
int *, int *, int *, double *, double *, double *, double *,
int *, double *), store2_(int *, double *, double *, int *,
int *, int *, double *, double *, double *, double *, int *),
setup2_(double *, double *, double *, double *, double *, double *, double *,
double *, double *, double *), givens_(double *, double *, double *, double *),
rotate_(int *, double *, double *, double *, double *);
static int nqwmax;
/****************************************************************/
/* */
/* ROBERT RENKA */
/* UNIV. OF NORTH TEXAS */
/*(817) 565 - 2767 */
/* 01/08/90 */
/* */
/* THIS SUBROUTINE COMPUTES A SET OF PARAMETERS A AND RSQ */
/* DEFINING A SMOOTH(ONCE CONTINUOUSLY DIFFERENTIABLE) BI- */
/* VARIATE FUNCTION Q(X, Y) WHICH INTERPOLATES DATA VALUES F */
/* AT SCATTERED NODES(X, Y). THE INTERPOLANT Q MAY BE EVAL- */
/* UATED AT AN ARBITRARY POINT BY FUNCTION QS2VAL, AND ITS */
/* FIRST DERIVATIVES ARE COMPUTED BY SUBROUTINE QS2GRD. */
/* THE INTERPOLATION SCHEME IS A MODIFIED QUADRATIC SHEPARD */
/* METHOD -- */
/* */
/* Q =(W(1)*Q(1) + W(2)*Q(2)+..+W(N)*Q(N))/(W(1) + W(2)+..+W(N))*/
/* */
/* FOR BIVARIATE FUNCTIONS W(K) AND Q(K). THE NODAL FUNC- */
/* TIONS ARE GIVEN BY */
/* */
/* Q(K)(X, Y) = A(1, K)*(X - X(K))**2 */
/* + A(2, K)*(X - X(K))*(Y - Y(K)) */
/* + A(3, K)*(Y - Y(K))**2 + A(4, K)*(X - X(K)) */
/* + A(5, K)*(Y - Y(K)) + F(K) . */
/* */
/* THUS, Q(K) IS A QUADRATIC FUNCTION WHICH INTERPOLATES THE */
/* DATA VALUE AT NODE K. ITS COEFFICIENTS A(, K) ARE OBTAINED */
/* BY A WEIGHTED LEAST SQUARES FIT TO THE CLOSEST NQ DATA */
/* POINTS WITH WEIGHTS SIMILAR TO W(K). NOTE THAT THE RADIUS */
/* OF INFLUENCE FOR THE LEAST SQUARES FIT IS FIXED FOR EACH */
/* K, BUT VARIES WITH K. */
/* THE WEIGHTS ARE TAKEN TO BE */
/* */
/* W(K)(X, Y) =((R(K) - D(K))+ / R(K)*D(K))**2 */
/* */
/* WHERE(R(K) - D(K))+ = 0 IF R(K) .LE. D(K) AND D(K)(X, Y) IS */
/* THE EUCLIDEAN DISTANCE BETWEEN(X, Y) AND(X(K), Y(K)). THE */
/* RADIUS OF INFLUENCE R(K) VARIES WITH K AND IS CHOSEN SO */
/* THAT NW NODES ARE WITHIN THE RADIUS. NOTE THAT W(K) IS */
/* NOT DEFINED AT NODE(X(K), Y(K)), BUT Q(X, Y) HAS LIMIT F(K) */
/* AS(X, Y) APPROACHES(X(K), Y(K)). */
/* */
/* ON INPUT -- */
/* */
/* N = NUMBER OF NODES AND ASSOCIATED DATA VALUES. */
/* N .GE. 6. */
/* */
/* X, Y = ARRAYS OF LENGTH N CONTAINING THE CARTESIAN */
/* COORDINATES OF THE NODES. */
/* */
/* F = ARRAY OF LENGTH N CONTAINING THE DATA VALUES */
/* IN ONE - TO - ONE CORRESPONDENCE WITH THE NODES. */
/* */
/* NQ = NUMBER OF DATA POINTS TO BE USED IN THE LEAST */
/* SQUARES FIT FOR COEFFICIENTS DEFINING THE NODAL */
/* FUNCTIONS Q(K). A HIGHLY RECOMMENDED VALUE IS */
/* NQ = 13. 5 .LE. NQ .LE. MIN(40, N - 1). */
/* */
/* NW = NUMBER OF NODES WITHIN(AND DEFINING) THE RADII */
/* OF INFLUENCE R(K) WHICH ENTER INTO THE WEIGHTS */
/* W(K). FOR N SUFFICIENTLY LARGE, A RECOMMENDED */
/* VALUE IS NW = 19. 1 .LE. NW .LE. MIN(40, N - 1). */
/* */
/* NR = NUMBER OF ROWS AND COLUMNS IN THE CELL GRID DE- */
/* FINED IN SUBROUTINE STORE2. A RECTANGLE CON- */
/* TAINING THE NODES IS PARTITIONED INTO CELLS IN */
/* ORDER TO INCREASE SEARCH EFFICIENCY. NR = */
/* SQRT(N/3) IS RECOMMENDED. NR .GE. 1. */
/* */
/* THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. */
/* */
/* LCELL = ARRAY OF LENGTH .GE. NR**2. */
/* */
/* LNEXT = ARRAY OF LENGTH .GE. N. */
/* */
/* RSQ = ARRAY OF LENGTH .GE. N. */
/* */
/* A = ARRAY OF LENGTH .GE. 5N. */
/* */
/* ON OUTPUT -- */
/* */
/* LCELL = NR BY NR ARRAY OF NODAL INDICES ASSOCIATED */
/* WITH CELLS. REFER TO STORE2. */
/* */
/* LNEXT = ARRAY OF LENGTH N CONTAINING NEXT - NODE INDI- */
/* CES. REFER TO STORE2. */
/* */
/* XMIN, YMIN, DX, DY = MINIMUM NODAL COORDINATES AND CELL */
/* DIMENSIONS. REFER TO STORE2. */
/* */
/* RMAX = SQUARE ROOT OF THE LARGEST ELEMENT IN RSQ -- */
/* MAXIMUM RADIUS R(K). */
/* */
/* RSQ = ARRAY CONTAINING THE SQUARES OF THE RADII R(K) */
/* WHICH ENTER INTO THE WEIGHTS W(K). */
/* */
/* A = 5 BY N ARRAY CONTAINING THE COEFFICIENTS FOR */
/* QUADRATIC NODAL FUNCTION Q(K) IN COLUMN K. */
/* */
/* NOTE THAT THE ABOVE OUTPUT PARAMETERS ARE NOT DEFINED */
/* UNLESS IER = 0. */
/* */
/* IER = ERROR INDICATOR -- */
/* IER = 0 IF NO ERRORS WERE ENCOUNTERED. */
/* IER = 1 IF N, NQ, NW, OR NR IS OUT OF RANGE. */
/* IER = 2 IF DUPLICATE NODES WERE ENCOUNTERED. */
/* IER = 3 IF ALL NODES ARE COLLINEAR. */
/* */
/* MODULES REQUIRED BY QSHEP2 -- GETNP2, GIVENS, ROTATE, */
/* SETUP2, STORE2 */
/* */
/* INTRINSIC FUNCTIONS CALLED BY QSHEP2 -- ABS, AMIN1, FLOAT, */
/* MAX0, MIN0, SQRT */
/* */
/****************************************************************/
a -= 6;
--rsq;
--lnext;
--f;
--y;
--x;
lcell_dim1 = *nr;
lcell_offset = 1 + lcell_dim1;
lcell -= lcell_offset;
/* LOCAL PARAMETERS -- */
/* */
/* AV = ROOT - MEAN - SQUARE DISTANCE BETWEEN K AND THE */
/* NODES IN THE LEAST SQUARES SYSTEM(UNLESS */
/* ADDITIONAL NODES ARE INTRODUCED FOR STABIL- */
/* ITY). THE FIRST 3 COLUMNS OF THE MATRIX */
/* ARE SCALED BY 1/AVSQ, THE LAST 2 BY 1/AV */
/* AVSQ = AV*AV */
/* B = TRANSPOSE OF THE AUGMENTED REGRESSION MATRIX */
/* C = FIRST COMPONENT OF THE PLANE ROTATION USED TO */
/* ZERO THE LOWER TRIANGLE OF B**T -- COMPUTED */
/* BY SUBROUTINE GIVENS */
/* DDX, DDY = LOCAL VARIABLES FOR DX AND DY */
/* DMIN = MINIMUM OF THE MAGNITUDES OF THE DIAGONAL */
/* ELEMENTS OF THE REGRESSION MATRIX AFTER */
/* ZEROS ARE INTRODUCED BELOW THE DIAGONAL */
/* DTOL = TOLERANCE FOR DETECTING AN ILL - CONDITIONED */
/* SYSTEM. THE SYSTEM IS ACCEPTED WHEN DMIN */
/* .GE. DTOL */
/* FK = DATA VALUE AT NODE K -- F(K) */
/* I = INDEX FOR A, B, AND NPTS */
/* IB = DO - LOOP INDEX FOR BACK SOLVE */
/* IERR = ERROR FLAG FOR THE CALL TO STORE2 */
/* IP1 = I + 1 */
/* IRM1 = IROW - 1 */
/* IROW = ROW INDEX FOR B */
/* J = INDEX FOR A AND B */
/* JP1 = J + 1 */
/* K = NODAL FUNCTION INDEX AND COLUMN INDEX FOR A */
/* LMAX = MAXIMUM NUMBER OF NPTS ELEMENTS(MUST BE CON- */
/* SISTENT WITH THE DIMENSION STATEMENT ABOVE) */
/* LNP = CURRENT LENGTH OF NPTS */
/* NEQ = NUMBER OF EQUATIONS IN THE LEAST SQUARES FIT */
/* NN, NNQ, NNR = LOCAL COPIES OF N, NQ, AND NR */
/* NNW = LOCAL COPY OF NW */
/* NP = NPTS ELEMENT */
/* NPTS = ARRAY CONTAINING THE INDICES OF A SEQUENCE OF */
/* NODES TO BE USED IN THE LEAST SQUARES FIT */
/* OR TO COMPUTE RSQ. THE NODES ARE ORDERED */
/* BY DISTANCE FROM K AND THE LAST ELEMENT */
/*(USUALLY INDEXED BY LNP) IS USED ONLY TO */
/* DETERMINE RQ, OR RSQ(K) IF NW .GT. NQ */
/* NQWMAX = MAX(NQ, NW) */
/* RQ = RADIUS OF INFLUENCE WHICH ENTERS INTO THE */
/* WEIGHTS FOR Q(K)(SEE SUBROUTINE SETUP2) */
/* RS = SQUARED DISTANCE BETWEEN K AND NPTS(LNP) -- */
/* USED TO COMPUTE RQ AND RSQ(K) */
/* RSMX = MAXIMUM RSQ ELEMENT ENCOUNTERED */
/* RSOLD = SQUARED DISTANCE BETWEEN K AND NPTS(LNP - 1) -- */
/* USED TO COMPUTE A RELATIVE CHANGE IN RS */
/* BETWEEN SUCCEEDING NPTS ELEMENTS */
/* RTOL = TOLERANCE FOR DETECTING A SUFFICIENTLY LARGE */
/* RELATIVE CHANGE IN RS. IF THE CHANGE IS */
/* NOT GREATER THAN RTOL, THE NODES ARE */
/* TREATED AS BEING THE SAME DISTANCE FROM K */
/* RWS = CURRENT VALUE OF RSQ(K) */
/* S = SECOND COMPONENT OF THE PLANE GIVENS ROTATION */
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -