PLplot  5.10.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
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 {
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 //
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 //
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 }