PLplot 5.15.0
Loading...
Searching...
No Matches
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
29typedef struct
30{
31 double w[3];
32} lweights;
33
34struct lpi
35{
38};
39
40int delaunay_xytoi( delaunay* d, point* p, int seed );
41
42// Builds linear interpolator.
43//
44// @param d Delaunay triangulation
45// @return Linear interpolator
46//
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//
100void 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//
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//
135void 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}
#define NaN
Definition csa/nan.h:62
int delaunay_xytoi(delaunay *d, point *p, int id)
Definition delaunay.c:631
delaunay * delaunay_build(int np, point points[], int ns, int segments[], int nh, double holes[])
Definition delaunay.c:265
void delaunay_destroy(delaunay *d)
Definition delaunay.c:578
void lpi_destroy(lpi *l)
Definition lpi.c:100
void lpi_interpolate_points(int nin, point pin[], int nout, point pout[])
Definition lpi.c:135
void lpi_interpolate_point(lpi *l, point *p)
Definition lpi.c:111
int delaunay_xytoi(delaunay *d, point *p, int seed)
Definition delaunay.c:631
lpi * lpi_build(delaunay *d)
Definition lpi.c:47
int nn_verbose
Definition nncommon.c:43
point * points
Definition delaunay.h:48
int ntriangles
Definition delaunay.h:54
int first_id
Definition delaunay.h:74
triangle * triangles
Definition delaunay.h:55
Definition lpi.c:35
delaunay * d
Definition lpi.c:36
lweights * weights
Definition lpi.c:37
Definition lpi.c:30
double w[3]
Definition lpi.c:31
Definition csa.h:30
double y
Definition csa.h:32
double x
Definition csa.h:31
double z
Definition csa.h:33
Definition csa.c:56
int vids[3]
Definition delaunay.h:25