1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265

#![no_std]
mod predicates;

/**
 * Adaptive exact 2D orientation test.  Robust.
 *
 * Return a positive value if the points `pa`, `pb`, and `pc` occur
 * in counterclockwise order; a negative value if they occur
 * in clockwise order; and zero if they are collinear.  The
 * result is also a rough approximation of twice the signed
 * area of the triangle defined by the three points.
 *
 * The result returned is the determinant of a matrix.
 * This determinant is computed adaptively, in the sense that exact
 * arithmetic is used only to the degree it is needed to ensure that the
 * returned value has the correct sign.  Hence, `orient2d()` is usually quite
 * fast, but will run more slowly when the input points are collinear or
 * nearly so.
 */
#[inline]
pub fn orient2d(pa: [f64; 2], pb: [f64; 2], pc: [f64; 2]) -> f64 {
    unsafe { predicates::orient2d(pa.as_ptr(), pb.as_ptr(), pc.as_ptr()) }
}

/**
 * Adaptive exact 3D orientation test. Robust.
 *
 * Return a positive value if the point `pd` lies below the
 * plane passing through `pa`, `pb`, and `pc`; "below" is defined so
 * that `pa`, `pb`, and `pc` appear in counterclockwise order when
 * viewed from above the plane.  Returns a negative value if
 * pd lies above the plane.  Returns zero if the points are
 * coplanar.  The result is also a rough approximation of six
 * times the signed volume of the tetrahedron defined by the
 * four points.
 *
 * The result returned is the determinant of a matrix.
 * This determinant is computed adaptively, in the sense that exact
 * arithmetic is used only to the degree it is needed to ensure that the
 * returned value has the correct sign.  Hence, orient3d() is usually quite
 * fast, but will run more slowly when the input points are coplanar or
 * nearly so.
 */
#[inline]
pub fn orient3d(pa: [f64; 3], pb: [f64; 3], pc: [f64; 3], pd: [f64; 3]) -> f64 {
    unsafe { predicates::orient3d(pa.as_ptr(), pb.as_ptr(), pc.as_ptr(), pd.as_ptr()) }
}

/**
 * Adaptive exaact 2D incircle test. Robust.
 *
 * Return a positive value if the point `pd` lies inside the
 * circle passing through `pa`, `pb`, and `pc`; a negative value if
 * it lies outside; and zero if the four points are cocircular.
 * The points `pa`, `pb`, and `pc` must be in counterclockwise
 * order, or the sign of the result will be reversed.
 *
 * The result returned is the determinant of a matrix.
 * This determinant is computed adaptively, in the sense that exact
 * arithmetic is used only to the degree it is needed to ensure that the
 * returned value has the correct sign.  Hence, `incircle()` is usually quite
 * fast, but will run more slowly when the input points are cocircular or
 * nearly so.
 */
#[inline]
pub fn incircle(pa: [f64; 2], pb: [f64; 2], pc: [f64; 2], pd: [f64; 2]) -> f64 {
    unsafe { predicates::incircle(pa.as_ptr(), pb.as_ptr(), pc.as_ptr(), pd.as_ptr()) }
}

/**
 * Adaptive exact 3D insphere test. Robust.
 *
 * Return a positive value if the point `pe` lies inside the
 * sphere passing through `pa`, `pb`, `pc`, and `pd`; a negative value
 * if it lies outside; and zero if the five points are
 * cospherical.  The points `pa`, `pb`, `pc`, and `pd` must be ordered
 * so that they have a positive orientation (as defined by
 * `orient3d()`), or the sign of the result will be reversed.
 *
 * The result returned is the determinant of a matrix.
 * this determinant is computed adaptively, in the sense that exact
 * arithmetic is used only to the degree it is needed to ensure that the
 * returned value has the correct sign.  Hence, `insphere()` is usually quite
 * fast, but will run more slowly when the input points are cospherical or
 * nearly so.
 */
#[inline]
pub fn insphere(pa: [f64; 3], pb: [f64; 3], pc: [f64; 3], pd: [f64; 3], pe: [f64; 3]) -> f64 {
    unsafe {
        predicates::insphere(
            pa.as_ptr(),
            pb.as_ptr(),
            pc.as_ptr(),
            pd.as_ptr(),
            pe.as_ptr(),
        )
    }
}

/// Approximate 2D orientation test. Nonrobust version of `orient2d`.
#[inline]
pub fn orient2d_fast(pa: [f64; 2], pb: [f64; 2], pc: [f64; 2]) -> f64 {
    unsafe { predicates::orient2dfast(pa.as_ptr(), pb.as_ptr(), pc.as_ptr()) }
}

/// Approximate 3D orientation test. Nonrobust version of `orient3d`.
#[inline]
pub fn orient3d_fast(pa: [f64; 3], pb: [f64; 3], pc: [f64; 3], pd: [f64; 3]) -> f64 {
    unsafe { predicates::orient3dfast(pa.as_ptr(), pb.as_ptr(), pc.as_ptr(), pd.as_ptr()) }
}

/// Approximate 2D incircle test. Nonrobust version of `incircle`
#[inline]
pub fn incircle_fast(pa: [f64; 2], pb: [f64; 2], pc: [f64; 2], pd: [f64; 2]) -> f64 {
    unsafe { predicates::incirclefast(pa.as_ptr(), pb.as_ptr(), pc.as_ptr(), pd.as_ptr()) }
}

/// Approximate 3D insphere test. Nonrobust version of `insphere`
#[inline]
pub fn insphere_fast(pa: [f64; 3], pb: [f64; 3], pc: [f64; 3], pd: [f64; 3], pe: [f64; 3]) -> f64 {
    unsafe {
        predicates::inspherefast(
            pa.as_ptr(),
            pb.as_ptr(),
            pc.as_ptr(),
            pd.as_ptr(),
            pe.as_ptr(),
        )
    }
}

#[cfg(test)]
mod tests {
    extern crate rand;
    use self::rand::{Rng, SeedableRng, StdRng};
    use super::*;

    const SEED: &[usize] = &[1, 2, 3, 4];

    #[test]
    fn orient2d_test() {
        let a = [0.0, 1.0];
        let b = [2.0, 3.0];
        let c = [4.0, 5.0];
        assert_eq!(orient2d(a, b, c), 0.0);
    }

    #[test]
    fn orient2d_robustness_test() {
        let mut rng: StdRng = SeedableRng::from_seed(SEED);
        let tol = 5.0e-14;

        #[cfg(miri)]
        let n = 99;
        #[cfg(not(miri))]
        let n = 9999999;
        for _ in 0..n {
            let a = [0.5 + tol * rng.gen::<f64>(), 0.5 + tol * rng.gen::<f64>()];
            let b = [12.0, 12.0];
            let c = [24.0, 24.0];
            assert_eq!(orient2d(a, b, c) > 0.0, orient2d(b, c, a) > 0.0);
            assert_eq!(orient2d(b, c, a) > 0.0, orient2d(c, a, b) > 0.0);
            assert_eq!(orient2d(a, b, c) > 0.0, orient2d(b, a, c) < 0.0);
            assert_eq!(orient2d(a, b, c) > 0.0, orient2d(a, c, b) < 0.0);
        }
    }

    #[test]
    fn orient2d_fast_test() {
        let a = [0.0, 1.0];
        let b = [2.0, 3.0];
        let c = [4.0, 5.0];
        assert_eq!(orient2d_fast(a, b, c), 0.0);

        // The fast orientation test should also work when the given points are sufficiently
        // non-colinear.
        let mut rng: StdRng = SeedableRng::from_seed(SEED);
        let tol = 5.0e-10; // will not work with 5.0e-14

        for _ in 0..999 {
            let a = [0.5 + tol * rng.gen::<f64>(), 0.5 + tol * rng.gen::<f64>()];
            let b = [12.0, 12.0];
            let c = [24.0, 24.0];
            assert_eq!(orient2d_fast(a, b, c) > 0.0, orient2d_fast(b, c, a) > 0.0);
            assert_eq!(orient2d_fast(b, c, a) > 0.0, orient2d_fast(c, a, b) > 0.0);
            assert_eq!(orient2d_fast(a, b, c) > 0.0, orient2d_fast(b, a, c) < 0.0);
            assert_eq!(orient2d_fast(a, b, c) > 0.0, orient2d_fast(a, c, b) < 0.0);
        }
    }

    #[test]
    fn orient3d_test() {
        let a = [0.0, 1.0, 6.0];
        let b = [2.0, 3.0, 4.0];
        let c = [4.0, 5.0, 1.0];
        let d = [6.0, 2.0, 5.3];
        assert_eq!(orient3d(a, b, c, d), 10.0);
    }

    #[test]
    fn orient3d_robustness_test() {
        let mut rng: StdRng = SeedableRng::from_seed(SEED);
        let tol = 5.0e-14;

        #[cfg(miri)]
        let n = 99;
        #[cfg(not(miri))]
        let n = 99999;
        for _ in 0..n {
            let a = [
                0.5 + tol * rng.gen::<f64>(),
                0.5 + tol * rng.gen::<f64>(),
                0.5 + tol * rng.gen::<f64>(),
            ];
            let b = [12.0, 12.0, 12.0];
            let c = [24.0, 24.0, 24.0];
            let d = [48.0, 48.0, 48.0];
            assert_eq!(orient3d(a, b, c, d) > 0.0, orient3d(b, c, d, a) > 0.0);
            assert_eq!(orient3d(b, c, d, a) > 0.0, orient3d(c, d, a, b) > 0.0);
            assert_eq!(orient3d(b, a, c, d) < 0.0, orient3d(c, b, d, a) < 0.0);
            assert_eq!(orient3d(b, d, c, a) < 0.0, orient3d(c, d, b, a) < 0.0);
        }
    }

    #[test]
    fn orient3d_fast_test() {
        let a = [0.0, 1.0, 6.0];
        let b = [2.0, 3.0, 4.0];
        let c = [4.0, 5.0, 1.0];
        let d = [6.0, 2.0, 5.3];
        assert!(f64::abs(orient3d_fast(a, b, c, d) - 10.0) < 1.0e-4);

        // The fast orientation test should also work when the given points are sufficiently
        // non-colinear.
        let mut rng: StdRng = SeedableRng::from_seed(SEED);
        let tol = 5.0; // will fail (see below) with all tolerances I tried

        let a = [
            0.5 + tol * rng.gen::<f64>(),
            0.5 + tol * rng.gen::<f64>(),
            0.5 + tol * rng.gen::<f64>(),
        ];
        let b = [12.0, 12.0, 12.0];
        let c = [24.0, 24.0, 24.0];
        let d = [48.0, 48.0, 48.0];
        assert_eq!(
            orient3d_fast(a, b, c, d) > 0.0,
            orient3d_fast(b, c, d, a) > 0.0
        );
        assert_eq!(
            orient3d_fast(b, a, c, d) < 0.0,
            orient3d_fast(c, b, d, a) < 0.0
        );
        // The following orientations are expected to be inconsistent
        // (this is why we need exact predicateexact predicatess)
        assert_ne!(
            orient3d_fast(b, c, d, a) > 0.0,
            orient3d_fast(c, d, a, b) > 0.0
        );
        assert_ne!(
            orient3d_fast(b, d, c, a) < 0.0,
            orient3d_fast(c, d, b, a) < 0.0
        );
    }
}