PLplot
5.10.0
|
00001 //-------------------------------------------------------------------------- 00002 // 00003 // File: linear.c 00004 // 00005 // Created: 04/08/2000 00006 // 00007 // Author: Pavel Sakov 00008 // CSIRO Marine Research 00009 // 00010 // Purpose: 2D linear interpolation 00011 // 00012 // Description: `lpi' -- "Linear Point Interpolator" -- is 00013 // a structure for conducting linear interpolation on a given 00014 // data on a "point-to-point" basis. It interpolates linearly 00015 // within each triangle resulted from the Delaunay 00016 // triangluation of input data. `lpi' is much 00017 // faster than all Natural Neighbours interpolators in `nn' 00018 // library. 00019 // 00020 // Revisions: None 00021 // 00022 //-------------------------------------------------------------------------- 00023 00024 #include <stdlib.h> 00025 #include <stdio.h> 00026 #include "nan.h" 00027 #include "delaunay.h" 00028 00029 typedef struct 00030 { 00031 double w[3]; 00032 } lweights; 00033 00034 struct lpi 00035 { 00036 delaunay* d; 00037 lweights* weights; 00038 }; 00039 00040 int delaunay_xytoi( delaunay* d, point* p, int seed ); 00041 00042 // Builds linear interpolator. 00043 // 00044 // @param d Delaunay triangulation 00045 // @return Linear interpolator 00046 // 00047 lpi* lpi_build( delaunay* d ) 00048 { 00049 int i; 00050 lpi * l = malloc( sizeof ( lpi ) ); 00051 00052 l->d = d; 00053 l->weights = malloc( (size_t) d->ntriangles * sizeof ( lweights ) ); 00054 00055 for ( i = 0; i < d->ntriangles; ++i ) 00056 { 00057 triangle* t = &d->triangles[i]; 00058 lweights* lw = &l->weights[i]; 00059 double x0 = d->points[t->vids[0]].x; 00060 double y0 = d->points[t->vids[0]].y; 00061 double z0 = d->points[t->vids[0]].z; 00062 double x1 = d->points[t->vids[1]].x; 00063 double y1 = d->points[t->vids[1]].y; 00064 double z1 = d->points[t->vids[1]].z; 00065 double x2 = d->points[t->vids[2]].x; 00066 double y2 = d->points[t->vids[2]].y; 00067 double z2 = d->points[t->vids[2]].z; 00068 double x02 = x0 - x2; 00069 double y02 = y0 - y2; 00070 double z02 = z0 - z2; 00071 double x12 = x1 - x2; 00072 double y12 = y1 - y2; 00073 double z12 = z1 - z2; 00074 00075 if ( y12 != 0.0 ) 00076 { 00077 double y0212 = y02 / y12; 00078 00079 lw->w[0] = ( z02 - z12 * y0212 ) / ( x02 - x12 * y0212 ); 00080 lw->w[1] = ( z12 - lw->w[0] * x12 ) / y12; 00081 lw->w[2] = ( z2 - lw->w[0] * x2 - lw->w[1] * y2 ); 00082 } 00083 else 00084 { 00085 double x0212 = x02 / x12; 00086 00087 lw->w[1] = ( z02 - z12 * x0212 ) / ( y02 - y12 * x0212 ); 00088 lw->w[0] = ( z12 - lw->w[1] * y12 ) / x12; 00089 lw->w[2] = ( z2 - lw->w[0] * x2 - lw->w[1] * y2 ); 00090 } 00091 } 00092 00093 return l; 00094 } 00095 00096 // Destroys linear interpolator. 00097 // 00098 // @param l Structure to be destroyed 00099 // 00100 void lpi_destroy( lpi* l ) 00101 { 00102 free( l->weights ); 00103 free( l ); 00104 } 00105 00106 // Finds linearly interpolated value in a point. 00107 // 00108 // @param l Linear interpolation 00109 // @param p Point to be interpolated (p->x, p->y -- input; p->z -- output) 00110 // 00111 void lpi_interpolate_point( lpi* l, point* p ) 00112 { 00113 delaunay* d = l->d; 00114 int tid = delaunay_xytoi( d, p, d->first_id ); 00115 00116 if ( tid >= 0 ) 00117 { 00118 lweights* lw = &l->weights[tid]; 00119 00120 d->first_id = tid; 00121 p->z = p->x * lw->w[0] + p->y * lw->w[1] + lw->w[2]; 00122 } 00123 else 00124 p->z = NaN; 00125 } 00126 00127 // Linearly interpolates data from one array of points for another array of 00128 // points. 00129 // 00130 // @param nin Number of input points 00131 // @param pin Array of input points [pin] 00132 // @param nout Number of ouput points 00133 // @param pout Array of output points [nout] 00134 // 00135 void lpi_interpolate_points( int nin, point pin[], int nout, point pout[] ) 00136 { 00137 delaunay* d = delaunay_build( nin, pin, 0, NULL, 0, NULL ); 00138 lpi * l = lpi_build( d ); 00139 int seed = 0; 00140 int i; 00141 00142 if ( nn_verbose ) 00143 { 00144 fprintf( stderr, "xytoi:\n" ); 00145 for ( i = 0; i < nout; ++i ) 00146 { 00147 point* p = &pout[i]; 00148 00149 fprintf( stderr, "(%.7g,%.7g) -> %d\n", p->x, p->y, delaunay_xytoi( d, p, seed ) ); 00150 } 00151 } 00152 00153 for ( i = 0; i < nout; ++i ) 00154 lpi_interpolate_point( l, &pout[i] ); 00155 00156 if ( nn_verbose ) 00157 { 00158 fprintf( stderr, "output:\n" ); 00159 for ( i = 0; i < nout; ++i ) 00160 { 00161 point* p = &pout[i];; 00162 fprintf( stderr, " %d:%15.7g %15.7g %15.7g\n", i, p->x, p->y, p->z ); 00163 } 00164 } 00165 00166 lpi_destroy( l ); 00167 delaunay_destroy( d ); 00168 }