PLplot
5.15.0
lpi.c
Go to the documentation of this file.
1
//--------------------------------------------------------------------------
2
//
3
// File: linear.c
4
//
5
// Created: 04/08/2000
6
//
7
// Author: Pavel Sakov
8
// CSIRO Marine Research
9
//
10
// Purpose: 2D linear interpolation
11
//
12
// Description: `lpi' -- "Linear Point Interpolator" -- is
13
// a structure for conducting linear interpolation on a given
14
// data on a "point-to-point" basis. It interpolates linearly
15
// within each triangle resulted from the Delaunay
16
// triangluation of input data. `lpi' is much
17
// faster than all Natural Neighbours interpolators in `nn'
18
// library.
19
//
20
// Revisions: None
21
//
22
//--------------------------------------------------------------------------
23
24
#include <stdlib.h>
25
#include <stdio.h>
26
#include "
nan.h
"
27
#include "
delaunay.h
"
28
29
typedef
struct
30
{
31
double
w[3];
32
}
lweights
;
33
34
struct
lpi
35
{
36
delaunay
*
d
;
37
lweights
*
weights
;
38
};
39
40
int
delaunay_xytoi
(
delaunay
* d,
point
* p,
int
seed );
41
42
// Builds linear interpolator.
43
//
44
// @param d Delaunay triangulation
45
// @return Linear interpolator
46
//
47
lpi
*
lpi_build
(
delaunay
* d )
48
{
49
int
i;
50
lpi
* l = malloc(
sizeof
(
lpi
) );
51
52
l->
d
= d;
53
l->
weights
= malloc( (
size_t
) d->
ntriangles
* sizeof (
lweights
) );
54
55
for
( i = 0; i < d->
ntriangles
; ++i )
56
{
57
triangle
* t = &d->
triangles
[i];
58
lweights
* lw = &l->
weights
[i];
59
double
x0 = d->
points
[t->
vids
[0]].
x
;
60
double
y0 = d->
points
[t->
vids
[0]].
y
;
61
double
z0 = d->
points
[t->
vids
[0]].
z
;
62
double
x1 = d->
points
[t->
vids
[1]].
x
;
63
double
y1 = d->
points
[t->
vids
[1]].
y
;
64
double
z1 = d->
points
[t->
vids
[1]].
z
;
65
double
x2 = d->
points
[t->
vids
[2]].
x
;
66
double
y2 = d->
points
[t->
vids
[2]].
y
;
67
double
z2 = d->
points
[t->
vids
[2]].
z
;
68
double
x02 = x0 - x2;
69
double
y02 = y0 - y2;
70
double
z02 = z0 - z2;
71
double
x12 = x1 - x2;
72
double
y12 = y1 - y2;
73
double
z12 = z1 - z2;
74
75
if
( y12 != 0.0 )
76
{
77
double
y0212 = y02 / y12;
78
79
lw->
w
[0] = ( z02 - z12 * y0212 ) / ( x02 - x12 * y0212 );
80
lw->
w
[1] = ( z12 - lw->
w
[0] * x12 ) / y12;
81
lw->
w
[2] = ( z2 - lw->
w
[0] * x2 - lw->
w
[1] * y2 );
82
}
83
else
84
{
85
double
x0212 = x02 / x12;
86
87
lw->
w
[1] = ( z02 - z12 * x0212 ) / ( y02 - y12 * x0212 );
88
lw->
w
[0] = ( z12 - lw->
w
[1] * y12 ) / x12;
89
lw->
w
[2] = ( z2 - lw->
w
[0] * x2 - lw->
w
[1] * y2 );
90
}
91
}
92
93
return
l;
94
}
95
96
// Destroys linear interpolator.
97
//
98
// @param l Structure to be destroyed
99
//
100
void
lpi_destroy
(
lpi
* l )
101
{
102
free( l->
weights
);
103
free( l );
104
}
105
106
// Finds linearly interpolated value in a point.
107
//
108
// @param l Linear interpolation
109
// @param p Point to be interpolated (p->x, p->y -- input; p->z -- output)
110
//
111
void
lpi_interpolate_point
(
lpi
* l,
point
* p )
112
{
113
delaunay
* d = l->
d
;
114
int
tid =
delaunay_xytoi
( d, p, d->
first_id
);
115
116
if
( tid >= 0 )
117
{
118
lweights
* lw = &l->
weights
[tid];
119
120
d->
first_id
= tid;
121
p->
z
= p->
x
* lw->
w
[0] + p->
y
* lw->
w
[1] + lw->
w
[2];
122
}
123
else
124
p->
z
=
NaN
;
125
}
126
127
// Linearly interpolates data from one array of points for another array of
128
// points.
129
//
130
// @param nin Number of input points
131
// @param pin Array of input points [pin]
132
// @param nout Number of ouput points
133
// @param pout Array of output points [nout]
134
//
135
void
lpi_interpolate_points
(
int
nin,
point
pin[],
int
nout,
point
pout[] )
136
{
137
delaunay
* d =
delaunay_build
( nin, pin, 0, NULL, 0, NULL );
138
lpi
* l =
lpi_build
( d );
139
int
seed = 0;
140
int
i;
141
142
if
(
nn_verbose
)
143
{
144
fprintf( stderr,
"xytoi:\n"
);
145
for
( i = 0; i < nout; ++i )
146
{
147
point
* p = &pout[i];
148
149
fprintf( stderr,
"(%.7g,%.7g) -> %d\n"
, p->
x
, p->
y
,
delaunay_xytoi
( d, p, seed ) );
150
}
151
}
152
153
for
( i = 0; i < nout; ++i )
154
lpi_interpolate_point
( l, &pout[i] );
155
156
if
(
nn_verbose
)
157
{
158
fprintf( stderr,
"output:\n"
);
159
for
( i = 0; i < nout; ++i )
160
{
161
point
* p = &pout[i];;
162
fprintf( stderr,
" %d:%15.7g %15.7g %15.7g\n"
, i, p->
x
, p->
y
, p->
z
);
163
}
164
}
165
166
lpi_destroy
( l );
167
delaunay_destroy
( d );
168
}
triangle
Definition:
csa.c:55
delaunay::ntriangles
int ntriangles
Definition:
delaunay.h:54
lpi_interpolate_point
void lpi_interpolate_point(lpi *l, point *p)
Definition:
lpi.c:111
delaunay_destroy
void delaunay_destroy(delaunay *d)
Definition:
delaunay.c:578
lpi::d
delaunay * d
Definition:
lpi.c:36
nan.h
delaunay.h
delaunay_build
delaunay * delaunay_build(int np, point points[], int ns, int segments[], int nh, double holes[])
Definition:
delaunay.c:265
delaunay::points
point * points
Definition:
delaunay.h:48
point::z
double z
Definition:
csa.h:33
delaunay::first_id
int first_id
Definition:
delaunay.h:74
lweights::w
double w[3]
Definition:
lpi.c:31
lweights
Definition:
lpi.c:29
NaN
#define NaN
Definition:
csa/nan.h:62
lpi
Definition:
lpi.c:34
point::y
double y
Definition:
csa.h:32
delaunay
Definition:
delaunay.h:45
point
Definition:
csa.h:29
lpi_build
lpi * lpi_build(delaunay *d)
Definition:
lpi.c:47
lpi_destroy
void lpi_destroy(lpi *l)
Definition:
lpi.c:100
lpi_interpolate_points
void lpi_interpolate_points(int nin, point pin[], int nout, point pout[])
Definition:
lpi.c:135
delaunay_xytoi
int delaunay_xytoi(delaunay *d, point *p, int seed)
Definition:
delaunay.c:631
point::x
double x
Definition:
csa.h:31
lpi::weights
lweights * weights
Definition:
lpi.c:37
nn_verbose
int nn_verbose
Definition:
nncommon.c:43
delaunay::triangles
triangle * triangles
Definition:
delaunay.h:55
triangle::vids
int vids[3]
Definition:
delaunay.h:25
lib
nn
lpi.c
Generated on Thu Nov 7 2019 00:00:00 for PLplot by
1.8.16