startin/interpolation/
mod.rs

1use crate::StartinError;
2use crate::Triangulation;
3use kdbush::KDBush;
4
5use crate::geom;
6
7pub trait Interpolant {
8    fn interpolate(
9        &self,
10        dt: &mut Triangulation,
11        locations: &Vec<[f64; 2]>,
12    ) -> Vec<Result<f64, StartinError>>;
13}
14
15pub fn interpolate(
16    interpolant: &impl Interpolant,
17    dt: &mut Triangulation,
18    locs: &Vec<[f64; 2]>,
19) -> Vec<Result<f64, StartinError>> {
20    interpolant.interpolate(dt, locs)
21}
22
23/// Estimation of z-value with interpolation: IDW
24/// (this function doesn't use the TIN at all, added here for
25/// convenience and teaching purposes)
26pub struct IDW {
27    pub radius: f64,
28    pub power: f64,
29}
30impl Interpolant for IDW {
31    fn interpolate(
32        &self,
33        dt: &mut Triangulation,
34        locs: &Vec<[f64; 2]>,
35    ) -> Vec<Result<f64, StartinError>> {
36        //-- build a kd-tree
37        let mut allpts: Vec<(f64, f64)> = Vec::new();
38        for i in 0..dt.stars.len() {
39            allpts.push((dt.stars[i].pt[0], dt.stars[i].pt[1]));
40        }
41        let index = KDBush::create(allpts, kdbush::DEFAULT_NODE_SIZE);
42        //-- perform interpolations
43        let mut re: Vec<Result<f64, StartinError>> = Vec::new();
44        for p in locs {
45            let mut ns: Vec<usize> = Vec::new();
46            index.within(p[0], p[1], self.radius, |id| ns.push(id));
47            if ns.is_empty() {
48                re.push(Err(StartinError::SearchCircleEmpty));
49            } else {
50                let mut weights: Vec<f64> = Vec::new();
51                let mut exisiting = false;
52                let mut value: f64 = std::f64::MAX;
53                for each in &ns {
54                    let d = geom::distance2d(p, &dt.stars[*each].pt);
55                    if d <= dt.get_snap_tolerance() {
56                        exisiting = true;
57                        value = dt.stars[*each].pt[2];
58                        break;
59                    }
60                    weights.push(d.powf(-self.power));
61                }
62                if exisiting {
63                    re.push(Ok(value));
64                } else {
65                    let mut z = 0_f64;
66                    for (i, w) in weights.iter().enumerate() {
67                        z += dt.stars[ns[i]].pt[2] * w;
68                    }
69                    re.push(Ok(z / weights.iter().sum::<f64>()));
70                }
71            }
72        }
73        re
74    }
75}
76
77/// Estimation of z-value with interpolation: Laplace interpolation
78///
79/// Details about Laplace: <http://dilbert.engr.ucdavis.edu/~suku/nem/index.html>, which
80/// is a variation of nni with distances instead of stolen areas, which yields a much
81/// faster implementation.
82pub struct Laplace {}
83impl Interpolant for Laplace {
84    fn interpolate(
85        &self,
86        dt: &mut Triangulation,
87        locs: &Vec<[f64; 2]>,
88    ) -> Vec<Result<f64, StartinError>> {
89        let mut re: Vec<Result<f64, StartinError>> = Vec::new();
90        for p in locs {
91            //-- cannot interpolate if no TIN
92            if !dt.is_init {
93                re.push(Err(StartinError::EmptyTriangulation));
94                continue;
95            }
96            //-- no extrapolation
97            let loc = dt.locate(p[0], p[1]);
98            match loc {
99                Ok(_tr) => {
100                    match dt.insert_one_pt_interpol(p[0], p[1]) {
101                        Ok(pi) => {
102                            //-- no extrapolation
103                            if dt.is_vertex_convex_hull(pi) {
104                                //-- interpolation point was added on boundary of CH
105                                //-- nothing to be done, Voronoi cell is unbounded
106                                let _rr = dt.remove(pi);
107                                re.push(Err(StartinError::OutsideConvexHull));
108                            } else {
109                                let l = &dt.stars[pi].link;
110                                let mut centres: Vec<Vec<f64>> = Vec::new();
111                                for (i, v) in l.iter().enumerate() {
112                                    let j = l.next_index(i);
113                                    centres.push(geom::circle_centre(
114                                        &dt.stars[pi].pt,
115                                        &dt.stars[*v].pt,
116                                        &dt.stars[l[j]].pt,
117                                    ));
118                                }
119                                let mut weights: Vec<f64> = Vec::new();
120                                for (i, v) in l.iter().enumerate() {
121                                    // fetch 2 voronoi centres
122                                    let e =
123                                        geom::distance2d(&centres[i], &centres[l.prev_index(i)]);
124                                    let w = geom::distance2d(&dt.stars[pi].pt, &dt.stars[*v].pt);
125                                    weights.push(e / w);
126                                }
127                                let mut z: f64 = 0.0;
128                                for (i, v) in l.iter().enumerate() {
129                                    z += weights[i] * dt.stars[*v].pt[2];
130                                }
131                                let sumweights: f64 = weights.iter().sum();
132                                //-- delete the interpolation location point
133                                let _rr = dt.remove(pi);
134                                re.push(Ok(z / sumweights));
135                            }
136                        }
137                        Err((pi, _updated)) => {
138                            re.push(Ok(dt.stars[pi].pt[2]));
139                        }
140                    }
141                }
142                Err(_e) => re.push(Err(StartinError::OutsideConvexHull)),
143            }
144        }
145        re
146    }
147}
148
149/// Estimation of z-value with interpolation: nearest/closest neighbour
150pub struct NN {}
151impl Interpolant for NN {
152    fn interpolate(
153        &self,
154        dt: &mut Triangulation,
155        locs: &Vec<[f64; 2]>,
156    ) -> Vec<Result<f64, StartinError>> {
157        let mut re: Vec<Result<f64, StartinError>> = Vec::new();
158        for p in locs {
159            //-- cannot interpolation if no TIN
160            if !dt.is_init {
161                re.push(Err(StartinError::EmptyTriangulation));
162                continue;
163            }
164            //-- TODO: should interpolate_nn() extrapolate?
165            match dt.closest_point(p[0], p[1]) {
166                Ok(vi) => re.push(Ok(dt.stars[vi].pt[2])),
167                Err(why) => re.push(Err(why)),
168            }
169        }
170        re
171    }
172}
173
174/// Estimation of z-value with interpolation: linear in TIN
175pub struct TIN {}
176impl Interpolant for TIN {
177    fn interpolate(
178        &self,
179        dt: &mut Triangulation,
180        locs: &Vec<[f64; 2]>,
181    ) -> Vec<Result<f64, StartinError>> {
182        let mut re: Vec<Result<f64, StartinError>> = Vec::new();
183        for p in locs {
184            //-- cannot interpolate if no TIN
185            if !dt.is_init {
186                re.push(Err(StartinError::EmptyTriangulation));
187                continue;
188            }
189            //-- no extrapolation
190            let loc = dt.locate(p[0], p[1]);
191            match loc {
192                Ok(tr) => {
193                    let q: [f64; 3] = [p[0], p[1], 0.0];
194                    let a0: f64 =
195                        geom::area2d_triangle(&q, &dt.stars[tr.v[1]].pt, &dt.stars[tr.v[2]].pt);
196                    let a1: f64 =
197                        geom::area2d_triangle(&q, &dt.stars[tr.v[2]].pt, &dt.stars[tr.v[0]].pt);
198                    let a2: f64 =
199                        geom::area2d_triangle(&q, &dt.stars[tr.v[0]].pt, &dt.stars[tr.v[1]].pt);
200                    let mut total = 0.;
201                    total += dt.stars[tr.v[0]].pt[2] * a0;
202                    total += dt.stars[tr.v[1]].pt[2] * a1;
203                    total += dt.stars[tr.v[2]].pt[2] * a2;
204                    re.push(Ok(total / (a0 + a1 + a2)));
205                }
206                Err(_e) => re.push(Err(StartinError::OutsideConvexHull)),
207            }
208        }
209        re
210    }
211}
212
213/// Estimation of z-value with interpolation: natural neighbour interpolation (nni),
214/// also called Sibson's interpolation
215pub struct NNI {
216    pub precompute: bool,
217}
218impl Interpolant for NNI {
219    fn interpolate(
220        &self,
221        dt: &mut Triangulation,
222        locs: &Vec<[f64; 2]>,
223    ) -> Vec<Result<f64, StartinError>> {
224        //-- store temporarily all the Voronoi cells areas
225        let mut vorareas: Vec<f64> = Vec::new();
226        if self.precompute {
227            vorareas.reserve_exact(dt.stars.len());
228            vorareas.push(0.);
229            for vi in 1..dt.stars.len() {
230                if !dt.stars[vi].is_deleted() {
231                    vorareas.push(dt.voronoi_cell_area(vi, true).unwrap());
232                } else {
233                    vorareas.push(0.);
234                }
235            }
236        }
237        let mut re: Vec<Result<f64, StartinError>> = Vec::new();
238        for p in locs {
239            //-- cannot interpolate if no TIN
240            if !dt.is_init {
241                re.push(Err(StartinError::EmptyTriangulation));
242                continue;
243            }
244            //-- no extrapolation
245            let loc = dt.locate(p[0], p[1]);
246            match loc {
247                Ok(_tr) => {
248                    match dt.insert_one_pt_interpol(p[0], p[1]) {
249                        Ok(pi) => {
250                            //-- no extrapolation
251                            if dt.is_vertex_convex_hull(pi) {
252                                //-- interpolation point was added on boundary of CH
253                                //-- nothing to be done, Voronoi cell is unbounded
254                                let _rr = dt.remove(pi);
255                                re.push(Err(StartinError::OutsideConvexHull));
256                            } else {
257                                let nns = dt.adjacent_vertices_to_vertex(pi).unwrap();
258                                let mut weights: Vec<f64> = Vec::new();
259                                for nn in &nns {
260                                    let a = dt.voronoi_cell_area(*nn, true).unwrap();
261                                    weights.push(a);
262                                }
263                                let newarea = dt.voronoi_cell_area(pi, true).unwrap();
264                                let _rr = dt.remove(pi);
265                                for (i, nn) in nns.iter().enumerate() {
266                                    if self.precompute {
267                                        weights[i] = vorareas[*nn] - weights[i];
268                                    } else {
269                                        //-- TODO : is it faster to save them?!
270                                        weights[i] =
271                                            dt.voronoi_cell_area(*nn, true).unwrap() - weights[i];
272                                    }
273                                }
274                                let mut z: f64 = 0.0;
275                                for (i, nn) in nns.iter().enumerate() {
276                                    z += weights[i] * dt.stars[*nn].pt[2];
277                                }
278                                re.push(Ok(z / newarea));
279                            }
280                        }
281                        Err(e) => re.push(Ok(dt.stars[e.0].pt[2])),
282                    }
283                }
284                Err(_e) => re.push(Err(StartinError::OutsideConvexHull)),
285            }
286        }
287        re
288    }
289}
startin/interpolation/mod.rs

startin/interpolation/
mod.rs
