diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-cubes.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-cubes.rkt index 15ab082e62..11b1d12457 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-cubes.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-cubes.rkt @@ -17,20 +17,20 @@ ;; edge vertexes -(define-syntax-rule (edge-1-2 d d1 d2) (vector (unsafe-solve-t d d1 d2) 0.0 0.0)) -(define-syntax-rule (edge-2-3 d d2 d3) (vector 1.0 (unsafe-solve-t d d2 d3) 0.0)) -(define-syntax-rule (edge-3-4 d d3 d4) (vector (unsafe-solve-t d d4 d3) 1.0 0.0)) -(define-syntax-rule (edge-1-4 d d1 d4) (vector 0.0 (unsafe-solve-t d d1 d4) 0.0)) +(define-syntax-rule (edge-1-2 d d1 d2) (vector (solve-t d d1 d2) 0 0)) +(define-syntax-rule (edge-2-3 d d2 d3) (vector 1 (solve-t d d2 d3) 0)) +(define-syntax-rule (edge-3-4 d d3 d4) (vector (solve-t d d4 d3) 1 0)) +(define-syntax-rule (edge-1-4 d d1 d4) (vector 0 (solve-t d d1 d4) 0)) -(define-syntax-rule (edge-5-6 d d5 d6) (vector (unsafe-solve-t d d5 d6) 0.0 1.0)) -(define-syntax-rule (edge-6-7 d d6 d7) (vector 1.0 (unsafe-solve-t d d6 d7) 1.0)) -(define-syntax-rule (edge-7-8 d d7 d8) (vector (unsafe-solve-t d d7 d8) 1.0 1.0)) -(define-syntax-rule (edge-5-8 d d5 d8) (vector 0.0 (unsafe-solve-t d d5 d8) 1.0)) +(define-syntax-rule (edge-5-6 d d5 d6) (vector (solve-t d d5 d6) 0 1)) +(define-syntax-rule (edge-6-7 d d6 d7) (vector 1 (solve-t d d6 d7) 1)) +(define-syntax-rule (edge-7-8 d d7 d8) (vector (solve-t d d7 d8) 1 1)) +(define-syntax-rule (edge-5-8 d d5 d8) (vector 0 (solve-t d d5 d8) 1)) -(define-syntax-rule (edge-1-5 d d1 d5) (vector 0.0 0.0 (unsafe-solve-t d d1 d5))) -(define-syntax-rule (edge-2-6 d d2 d6) (vector 1.0 0.0 (unsafe-solve-t d d2 d6))) -(define-syntax-rule (edge-3-7 d d3 d7) (vector 1.0 1.0 (unsafe-solve-t d d3 d7))) -(define-syntax-rule (edge-4-8 d d4 d8) (vector 0.0 1.0 (unsafe-solve-t d d4 d8))) +(define-syntax-rule (edge-1-5 d d1 d5) (vector 0 0 (solve-t d d1 d5))) +(define-syntax-rule (edge-2-6 d d2 d6) (vector 1 0 (solve-t d d2 d6))) +(define-syntax-rule (edge-3-7 d d3 d7) (vector 1 1 (solve-t d d3 d7))) +(define-syntax-rule (edge-4-8 d d4 d8) (vector 0 1 (solve-t d d4 d8))) #| Cube vertex numbers: @@ -72,7 +72,7 @@ Cube vertex numbers: (define known-cube-0000-1111 known-cube-1111-0000) (define ((make-known-cube-1010-0000 test?) d d1 d2 d3 d4 d5 d6 d7 d8) - (define da (unsafe-flavg4 d1 d2 d3 d4)) + (define da (/ (+ d1 d2 d3 d4) 4)) (cond [(test? da d) (list (list (edge-1-2 d d1 d2) (edge-2-3 d d2 d3) @@ -83,8 +83,8 @@ Cube vertex numbers: (list (list (edge-1-2 d d1 d2) (edge-1-5 d d1 d5) (edge-1-4 d d1 d4)) (list (edge-2-3 d d2 d3) (edge-3-4 d d3 d4) (edge-3-7 d d3 d7)))])) -(define known-cube-1010-0000 (make-known-cube-1010-0000 unsafe-fl>=)) -(define known-cube-0101-1111 (make-known-cube-1010-0000 unsafe-fl<)) +(define known-cube-1010-0000 (make-known-cube-1010-0000 >=)) +(define known-cube-0101-1111 (make-known-cube-1010-0000 <)) (define (known-cube-1000-0010 d d1 d2 d3 d4 d5 d6 d7 d8) (list (list (edge-1-2 d d1 d2) (edge-1-5 d d1 d5) (edge-1-4 d d1 d4)) @@ -93,7 +93,7 @@ Cube vertex numbers: (define known-cube-0111-1101 known-cube-1000-0010) (define ((make-known-cube-1100-0010 test?) d d1 d2 d3 d4 d5 d6 d7 d8) - (define da (unsafe-flavg4 d2 d6 d7 d3)) + (define da (/ (+ d2 d6 d7 d3) 4)) (cond [(test? da d) (list (list (edge-1-5 d d1 d5) (edge-2-6 d d2 d6) @@ -106,12 +106,12 @@ Cube vertex numbers: (edge-2-3 d d2 d3) (edge-1-4 d d1 d4)) (list (edge-6-7 d d6 d7) (edge-3-7 d d3 d7) (edge-7-8 d d7 d8)))])) -(define known-cube-1100-0010 (make-known-cube-1100-0010 unsafe-fl>=)) -(define known-cube-0011-1101 (make-known-cube-1100-0010 unsafe-fl<)) +(define known-cube-1100-0010 (make-known-cube-1100-0010 >=)) +(define known-cube-0011-1101 (make-known-cube-1100-0010 <)) (define ((make-known-cube-1100-0011 test?) d d1 d2 d3 d4 d5 d6 d7 d8) - (define da (unsafe-flavg4 d1 d5 d8 d4)) - (define db (unsafe-flavg4 d2 d6 d7 d3)) + (define da (/ (+ d1 d5 d8 d4) 4)) + (define db (/ (+ d2 d6 d7 d3) 4)) (cond [(and (test? da d) (test? db d)) (list (list (edge-1-5 d d1 d5) (edge-5-8 d d5 d8) @@ -123,7 +123,7 @@ Cube vertex numbers: (v+ (edge-4-8 d d4 d8) (edge-5-8 d d5 d8))) (v+ (v+ (edge-6-7 d d6 d7) (edge-2-6 d d2 d6)) (v+ (edge-2-3 d d2 d3) (edge-3-7 d d3 d7)))) - 0.125)) + 1/8)) (list (list ec (edge-5-8 d d5 d8) (edge-6-7 d d6 d7)) (list ec (edge-1-5 d d1 d5) (edge-2-6 d d2 d6)) (list ec (edge-4-8 d d4 d8) (edge-3-7 d d3 d7)) @@ -137,7 +137,7 @@ Cube vertex numbers: (v+ (edge-4-8 d d4 d8) (edge-5-8 d d5 d8))) (v+ (v+ (edge-6-7 d d6 d7) (edge-2-6 d d2 d6)) (v+ (edge-2-3 d d2 d3) (edge-3-7 d d3 d7)))) - 0.125)) + 1/8)) (list (list ec (edge-5-8 d d5 d8) (edge-6-7 d d6 d7)) (list ec (edge-1-5 d d1 d5) (edge-2-6 d d2 d6)) (list ec (edge-4-8 d d4 d8) (edge-3-7 d d3 d7)) @@ -152,8 +152,8 @@ Cube vertex numbers: (list (edge-5-8 d d5 d8) (edge-6-7 d d6 d7) (edge-3-7 d d3 d7) (edge-4-8 d d4 d8)))])) -(define known-cube-1100-0011 (make-known-cube-1100-0011 unsafe-fl>=)) -(define known-cube-0011-1100 (make-known-cube-1100-0011 unsafe-fl<)) +(define known-cube-1100-0011 (make-known-cube-1100-0011 >=)) +(define known-cube-0011-1100 (make-known-cube-1100-0011 <)) #| Cube vertex numbers: @@ -169,12 +169,12 @@ Cube vertex numbers: (define ((make-known-cube-1010-0101 test?) d d1 d2 d3 d4 d5 d6 d7 d8) - (define da (unsafe-flavg4 d1 d2 d3 d4)) - (define db (unsafe-flavg4 d1 d5 d8 d4)) - (define dc (unsafe-flavg4 d3 d4 d8 d7)) - (define dd (unsafe-flavg4 d1 d2 d6 d5)) - (define de (unsafe-flavg4 d2 d3 d7 d6)) - (define df (unsafe-flavg4 d5 d6 d7 d8)) + (define da (/ (+ d1 d2 d3 d4) 4)) + (define db (/ (+ d1 d5 d8 d4) 4)) + (define dc (/ (+ d3 d4 d8 d7) 4)) + (define dd (/ (+ d1 d2 d6 d5) 4)) + (define de (/ (+ d2 d3 d7 d6) 4)) + (define df (/ (+ d5 d6 d7 d8) 4)) (append (list (list (edge-1-5 d d1 d5) (edge-1-2 d d1 d2) (edge-1-4 d d1 d4)) (list (edge-2-3 d d2 d3) (edge-3-7 d d3 d7) (edge-3-4 d d3 d4)) @@ -205,26 +205,26 @@ Cube vertex numbers: (edge-7-8 d d7 d8) (edge-5-8 d d5 d8))) empty))) -(define known-cube-1010-0101 (make-known-cube-1010-0101 unsafe-fl>=)) -(define known-cube-0101-1010 (make-known-cube-1010-0101 unsafe-fl<)) +(define known-cube-1010-0101 (make-known-cube-1010-0101 >=)) +(define known-cube-0101-1010 (make-known-cube-1010-0101 <)) -(define ((make-known-cube-1110-0001 unsafe-fl>=) d d1 d2 d3 d4 d5 d6 d7 d8) - (define da (unsafe-flavg4 d1 d5 d8 d4)) - (define db (unsafe-flavg4 d7 d8 d4 d3)) +(define ((make-known-cube-1110-0001 >=) d d1 d2 d3 d4 d5 d6 d7 d8) + (define da (/ (+ d1 d5 d8 d4) 4)) + (define db (/ (+ d7 d8 d4 d3) 4)) (cond - [(and (da . unsafe-fl>= . d) (db . unsafe-fl>= . d)) + [(and (da . >= . d) (db . >= . d)) (list (list (edge-1-5 d d1 d5) (edge-2-6 d d2 d6) (edge-3-7 d d3 d7)) (list (edge-1-5 d d1 d5) (edge-3-7 d d3 d7) (edge-7-8 d d7 d8) (edge-5-8 d d5 d8)) (list (edge-1-4 d d1 d4) (edge-3-4 d d3 d4) (edge-4-8 d d4 d8)))] - [(da . unsafe-fl>= . d) - (define ec (v* (v+ (edge-1-5 d d1 d5) (edge-3-7 d d3 d7)) 0.5)) + [(da . >= . d) + (define ec (v* (v+ (edge-1-5 d d1 d5) (edge-3-7 d d3 d7)) 1/2)) (list (list (edge-1-5 d d1 d5) (edge-2-6 d d2 d6) (edge-3-7 d d3 d7)) (list ec (edge-3-7 d d3 d7) (edge-3-4 d d3 d4) (edge-1-4 d d1 d4)) (list ec (edge-1-5 d d1 d5) (edge-5-8 d d5 d8) (edge-7-8 d d7 d8)) (list ec (edge-1-4 d d1 d4) (edge-4-8 d d4 d8) (edge-7-8 d d7 d8)))] - [(db . unsafe-fl>= . d) - (define ec (v* (v+ (edge-1-5 d d1 d5) (edge-3-7 d d3 d7)) 0.5)) + [(db . >= . d) + (define ec (v* (v+ (edge-1-5 d d1 d5) (edge-3-7 d d3 d7)) 1/2)) (list (list (edge-1-5 d d1 d5) (edge-2-6 d d2 d6) (edge-3-7 d d3 d7)) (list ec (edge-1-5 d d1 d5) (edge-1-4 d d1 d4) (edge-3-4 d d3 d4)) (list ec (edge-3-7 d d3 d7) (edge-7-8 d d7 d8) (edge-5-8 d d5 d8)) @@ -235,8 +235,8 @@ Cube vertex numbers: (edge-3-4 d d3 d4) (edge-1-4 d d1 d4)) (list (edge-7-8 d d7 d8) (edge-5-8 d d5 d8) (edge-4-8 d d4 d8)))])) -(define known-cube-1110-0001 (make-known-cube-1110-0001 unsafe-fl>=)) -(define known-cube-0001-1110 (make-known-cube-1110-0001 unsafe-fl<)) +(define known-cube-1110-0001 (make-known-cube-1110-0001 >=)) +(define known-cube-0001-1110 (make-known-cube-1110-0001 <)) (define (known-cube-1110-0100 d d1 d2 d3 d4 d5 d6 d7 d8) (list (list (edge-1-5 d d1 d5) (edge-5-6 d d5 d6) (edge-6-7 d d6 d7) @@ -251,9 +251,9 @@ Cube vertex numbers: (define known-cube-0001-1101 known-cube-1110-0010) (define ((make-known-cube-1010-0001 test?) d d1 d2 d3 d4 d5 d6 d7 d8) - (define da (unsafe-flavg4 d1 d2 d3 d4)) - (define db (unsafe-flavg4 d1 d5 d8 d4)) - (define dc (unsafe-flavg4 d3 d4 d8 d7)) + (define da (/ (+ d1 d2 d3 d4) 4)) + (define db (/ (+ d1 d5 d8 d4) 4)) + (define dc (/ (+ d3 d4 d8 d7) 4)) (append (list (list (edge-1-5 d d1 d5) (edge-1-2 d d1 d2) (edge-1-4 d d1 d4)) (list (edge-2-3 d d2 d3) (edge-3-7 d d3 d7) (edge-3-4 d d3 d4)) @@ -271,8 +271,8 @@ Cube vertex numbers: (edge-7-8 d d7 d8) (edge-4-8 d d4 d8))) empty))) -(define known-cube-1010-0001 (make-known-cube-1010-0001 unsafe-fl>=)) -(define known-cube-0101-1110 (make-known-cube-1010-0001 unsafe-fl<)) +(define known-cube-1010-0001 (make-known-cube-1010-0001 >=)) +(define known-cube-0101-1110 (make-known-cube-1010-0001 <)) (define-for-syntax known-cubes '((0 0 0 0 0 0 0 0) @@ -308,7 +308,7 @@ Cube vertex numbers: (define (mirror-vec-d v) (match-define (vector x y d) v) - (vector x y (unsafe-fl- 1.0 d))) + (vector x y (- 1 d))) (define ((mirror-cube-d f) d d1 d2 d3 d4 d5 d6 d7 d8) (map (λ (poly) (map mirror-vec-d poly)) @@ -316,7 +316,7 @@ Cube vertex numbers: (define (mirror-vec-y v) (match-define (vector x y d) v) - (vector x (unsafe-fl- 1.0 y) d)) + (vector x (- 1 y) d)) (define ((mirror-cube-y f) d d1 d2 d3 d4 d5 d6 d7 d8) (map (λ (poly) (map mirror-vec-y poly)) @@ -324,7 +324,7 @@ Cube vertex numbers: (define (mirror-vec-x v) (match-define (vector x y d) v) - (vector (unsafe-fl- 1.0 x) y d)) + (vector (- 1 x) y d)) (define ((mirror-cube-x f) d d1 d2 d3 d4 d5 d6 d7 d8) (map (λ (poly) (map mirror-vec-x poly)) @@ -334,7 +334,7 @@ Cube vertex numbers: (define (unrotate-vec-d v) (match-define (vector x y d) v) - (vector (unsafe-fl- 1.0 y) x d)) + (vector (- 1 y) x d)) (define ((rotate-cube-d f) d d1 d2 d3 d4 d5 d6 d7 d8) (map (λ (poly) (map unrotate-vec-d poly)) @@ -342,7 +342,7 @@ Cube vertex numbers: (define (unrotate-vec-y v) (match-define (vector x y d) v) - (vector d y (unsafe-fl- 1.0 x))) + (vector d y (- 1 x))) (define ((rotate-cube-y f) d d1 d2 d3 d4 d5 d6 d7 d8) (map (λ (poly) (map unrotate-vec-y poly)) @@ -350,7 +350,7 @@ Cube vertex numbers: (define (unrotate-vec-x v) (match-define (vector x y d) v) - (vector x (unsafe-fl- 1.0 d) y)) + (vector x (- 1 d) y)) (define ((rotate-cube-x f) d d1 d2 d3 d4 d5 d6 d7 d8) (map (λ (poly) (map unrotate-vec-x poly)) @@ -360,7 +360,7 @@ Cube vertex numbers: (define (rotate-vec-d v) (match-define (vector x y d) v) - (vector y (unsafe-fl- 1.0 x) d)) + (vector y (- 1 x) d)) (define ((unrotate-cube-d f) d d1 d2 d3 d4 d5 d6 d7 d8) (map (λ (poly) (map rotate-vec-d poly)) @@ -368,7 +368,7 @@ Cube vertex numbers: (define (rotate-vec-y v) (match-define (vector x y d) v) - (vector (unsafe-fl- 1.0 d) y x)) + (vector (- 1 d) y x)) (define ((unrotate-cube-y f) d d1 d2 d3 d4 d5 d6 d7 d8) (map (λ (poly) (map rotate-vec-y poly)) @@ -376,7 +376,7 @@ Cube vertex numbers: (define (rotate-vec-x v) (match-define (vector x y d) v) - (vector x d (unsafe-fl- 1.0 y))) + (vector x d (- 1 y))) (define ((unrotate-cube-x f) d d1 d2 d3 d4 d5 d6 d7 d8) (map (λ (poly) (map rotate-vec-x poly)) @@ -511,15 +511,15 @@ Cube vertex numbers: (define-syntax-rule (add-digit j idx) (unsafe-fx+ (unsafe-fx* idx 2) j)) -(define (unsafe-heights->cube-polys d d1 d2 d3 d4 d5 d6 d7 d8) - (define j1 (if (d1 . unsafe-fl< . d) 0 1)) - (define j2 (if (d2 . unsafe-fl< . d) 0 1)) - (define j3 (if (d3 . unsafe-fl< . d) 0 1)) - (define j4 (if (d4 . unsafe-fl< . d) 0 1)) - (define j5 (if (d5 . unsafe-fl< . d) 0 1)) - (define j6 (if (d6 . unsafe-fl< . d) 0 1)) - (define j7 (if (d7 . unsafe-fl< . d) 0 1)) - (define j8 (if (d8 . unsafe-fl< . d) 0 1)) +(define (do-heights->cube-polys d d1 d2 d3 d4 d5 d6 d7 d8) + (define j1 (if (d1 . < . d) 0 1)) + (define j2 (if (d2 . < . d) 0 1)) + (define j3 (if (d3 . < . d) 0 1)) + (define j4 (if (d4 . < . d) 0 1)) + (define j5 (if (d5 . < . d) 0 1)) + (define j6 (if (d6 . < . d) 0 1)) + (define j7 (if (d7 . < . d) 0 1)) + (define j8 (if (d8 . < . d) 0 1)) (define facet-num (add-digit j1 (add-digit @@ -536,37 +536,12 @@ Cube vertex numbers: [d1 real?] [d2 real?] [d3 real?] [d4 real?] [d5 real?] [d6 real?] [d7 real?] [d8 real?] ) (listof (listof (vector/c real? real? real?))) - (cond [(all inexact-real? xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8) - (define polys (unsafe-heights->cube-polys d d1 d2 d3 d4 d5 d6 d7 d8)) + (cond [(= d d1 d2 d3 d4 d5 d6 d7 d8) empty] + [else + (define polys (do-heights->cube-polys d d1 d2 d3 d4 d5 d6 d7 d8)) (for/list ([poly (in-list polys)]) (for/list ([uvw (in-list poly)]) (match-define (vector u v w) uvw) - (vector (unsafe-unsolve-t xa xb u) - (unsafe-unsolve-t ya yb v) - (unsafe-unsolve-t za zb w))))] - [(find-failure-index real? xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8) - => (λ (i) (raise-type-error 'heights->polys "real number" - i xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))] - [(= d d1 d2 d3 d4 d5 d6 d7 d8) empty] - [else - (let-map - (d d1 d2 d3 d4 d5 d6 d7 d8) inexact->exact - (define d-min (min d d1 d2 d3 d4 d5 d6 d7 d8)) - (define d-max (max d d1 d2 d3 d4 d5 d6 d7 d8)) - (define d-scale (- d-max d-min)) - (define polys - (unsafe-heights->cube-polys (exact->inexact (/ (- d d-min) d-scale)) - (exact->inexact (/ (- d1 d-min) d-scale)) - (exact->inexact (/ (- d2 d-min) d-scale)) - (exact->inexact (/ (- d3 d-min) d-scale)) - (exact->inexact (/ (- d4 d-min) d-scale)) - (exact->inexact (/ (- d5 d-min) d-scale)) - (exact->inexact (/ (- d6 d-min) d-scale)) - (exact->inexact (/ (- d7 d-min) d-scale)) - (exact->inexact (/ (- d8 d-min) d-scale)))) - (for/list ([poly (in-list polys)]) - (for/list ([uvw (in-list poly)]) - (match-define (vector u v w) uvw) - (vector (unsolve-t xa xb u) - (unsolve-t ya yb v) - (unsolve-t za zb w)))))])) + (vector (unsolve-t xa xb u) + (unsolve-t ya yb v) + (unsolve-t za zb w))))])) diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-squares.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-squares.rkt index 4bed76be26..09598c2826 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-squares.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-squares.rkt @@ -2,6 +2,7 @@ (require racket/flonum racket/fixnum racket/list racket/match racket/unsafe/ops racket/contract unstable/latent-contract/defthing + (only-in math/flonum fl) (for-syntax racket/base) "math.rkt" "utils.rkt" @@ -29,18 +30,6 @@ first argument is the contour value; the rest are z coordinates arranged as above. |# -(define-syntax-rule (unrotate-vec v) - (match-let ([(vector x y z) v]) - (vector (unsafe-fl- 1.0 y) x z))) - -(define-syntax-rule (mirror-x-vec v) - (match-let ([(vector x y z) v]) - (vector (fl- 1.0 x) y z))) - -(define-syntax-rule (mirror-y-vec v) - (match-let ([(vector x y z) v]) - (vector x (fl- 1.0 y) z))) - ;; ============================================================================= ;; Contour lines @@ -59,23 +48,23 @@ above. ;(: lines1000 (FType Lines)) (define-syntax-rule (lines1000 z z1 z2 z3 z4) - (list (vector (unsafe-solve-t z z1 z2) 0.0 - 0.0 (unsafe-solve-t z z1 z4)))) + (list (vector (solve-t z z1 z2) 0 + 0 (solve-t z z1 z4)))) ;(: lines0100 (FType Lines)) (define-syntax-rule (lines0100 z z1 z2 z3 z4) - (list (vector (unsafe-solve-t z z1 z2) 0.0 - 1.0 (unsafe-solve-t z z2 z3)))) + (list (vector (solve-t z z1 z2) 0 + 1 (solve-t z z2 z3)))) ;(: lines0010 (FType Lines)) (define-syntax-rule (lines0010 z z1 z2 z3 z4) - (list (vector 1.0 (unsafe-solve-t z z2 z3) - (unsafe-solve-t z z4 z3) 1.0))) + (list (vector 1 (solve-t z z2 z3) + (solve-t z z4 z3) 1))) ;(: lines0001 (FType Lines)) (define-syntax-rule (lines0001 z z1 z2 z3 z4) - (list (vector 0.0 (unsafe-solve-t z z1 z4) - (unsafe-solve-t z z4 z3) 1.0))) + (list (vector 0 (solve-t z z1 z4) + (solve-t z z4 z3) 1))) (define-syntax-rule (lines0111 z z1 z2 z3 z4) (lines1000 z z1 z2 z3 z4)) (define-syntax-rule (lines1011 z z1 z2 z3 z4) (lines0100 z z1 z2 z3 z4)) @@ -87,13 +76,13 @@ above. ;(: lines1100 (FType Lines)) (define-syntax-rule (lines1100 z z1 z2 z3 z4) - (list (vector 0.0 (unsafe-solve-t z z1 z4) - 1.0 (unsafe-solve-t z z2 z3)))) + (list (vector 0 (solve-t z z1 z4) + 1 (solve-t z z2 z3)))) ;(: lines0110 (FType Lines)) (define-syntax-rule (lines0110 z z1 z2 z3 z4) - (list (vector (unsafe-solve-t z z1 z2) 0.0 - (unsafe-solve-t z z4 z3) 1.0))) + (list (vector (solve-t z z1 z2) 0 + (solve-t z z4 z3) 1))) (define-syntax-rule (lines0011 z z1 z2 z3 z4) (lines1100 z z1 z2 z3 z4)) (define-syntax-rule (lines1001 z z1 z2 z3 z4) (lines0110 z z1 z2 z3 z4)) @@ -104,25 +93,27 @@ above. ;(: lines-opposite ((Float Float -> Boolean) -> (FType Lines))) (define-syntax-rule (lines-opposite test? z z1 z2 z3 z4) ; disambiguate using average of corners as guess for center value - (let ([z5 (unsafe-flavg4 z1 z2 z3 z4)]) + (let ([z5 (/ (+ z1 z2 z3 z4) 4)]) (if (test? z5 z) - (list (vector (unsafe-solve-t z z1 z2) 0.0 - 1.0 (unsafe-solve-t z z2 z3)) - (vector 0.0 (unsafe-solve-t z z1 z4) - (unsafe-solve-t z z4 z3) 1.0)) - (list (vector (unsafe-solve-t z z1 z2) 0.0 - 0.0 (unsafe-solve-t z z1 z4)) - (vector 1.0 (unsafe-solve-t z z2 z3) - (unsafe-solve-t z z4 z3) 1.0))))) + (list (vector (solve-t z z1 z2) 0 + 1 (solve-t z z2 z3)) + (vector 0 (solve-t z z1 z4) + (solve-t z z4 z3) 1)) + (list (vector (solve-t z z1 z2) 0 + 0 (solve-t z z1 z4)) + (vector 1 (solve-t z z2 z3) + (solve-t z z4 z3) 1))))) -(define-syntax-rule (lines1010 z z1 z2 z3 z4) (lines-opposite unsafe-fl>= z z1 z2 z3 z4)) -(define-syntax-rule (lines0101 z z1 z2 z3 z4) (lines-opposite unsafe-fl< z z1 z2 z3 z4)) +(define-syntax-rule (lines1010 z z1 z2 z3 z4) (lines-opposite >= z z1 z2 z3 z4)) +(define-syntax-rule (lines0101 z z1 z2 z3 z4) (lines-opposite < z z1 z2 z3 z4)) -(define (unsafe-heights->lines z z1 z2 z3 z4) - (define p1 (z1 . unsafe-fl>= . z)) - (define p2 (z2 . unsafe-fl>= . z)) - (define p3 (z3 . unsafe-fl>= . z)) - (define p4 (z4 . unsafe-fl>= . z)) +(define (do-heights->lines z z1 z2 z3 z4) + ;(printf "~v ~v ~v ~v ~v~n" z z1 z2 z3 z4) + (define p1 (z1 . >= . z)) + (define p2 (z2 . >= . z)) + (define p3 (z3 . >= . z)) + (define p4 (z4 . >= . z)) + ;(printf "~v ~v ~v ~v~n~n" p1 p2 p3 p4) (if p1 (if p2 (if p3 @@ -158,46 +149,45 @@ above. (defproc (heights->lines [xa real?] [xb real?] [ya real?] [yb real?] [z real?] [z1 real?] [z2 real?] [z3 real?] [z4 real?] ) (listof (list/c (vector/c real? real? real?) (vector/c real? real? real?))) - (cond [(all inexact-real? xa xb ya yb z z1 z2 z3 z4) - (define lines (unsafe-heights->lines z z1 z2 z3 z4)) - (for/list ([line (in-list lines)]) - (match-define (vector u1 v1 u2 v2) line) - (list (vector (unsafe-unsolve-t xa xb u1) (unsafe-unsolve-t ya yb v1) z) - (vector (unsafe-unsolve-t xa xb u2) (unsafe-unsolve-t ya yb v2) z)))] - [(find-failure-index real? xa xb ya yb z z1 z2 z3 z4) - => (λ (i) (raise-type-error 'heights->liens "real number" i xa xb ya yb z z1 z2 z3 z4))] + (cond [(and (z . < . z1) (z . < . z2) (z . < . z3) (z . < . z4)) empty] + [(and (z . > . z1) (z . > . z2) (z . > . z3) (z . > . z4)) empty] [(= z z1 z2 z3 z4) empty] [else - (let-map - (z z1 z2 z3 z4) inexact->exact - (define z-min (min z z1 z2 z3 z4)) - (define z-max (max z z1 z2 z3 z4)) - (define z-scale (- z-max z-min)) - (define lines - (unsafe-heights->lines (exact->inexact (/ (- z z-min) z-scale)) - (exact->inexact (/ (- z1 z-min) z-scale)) - (exact->inexact (/ (- z2 z-min) z-scale)) - (exact->inexact (/ (- z3 z-min) z-scale)) - (exact->inexact (/ (- z4 z-min) z-scale)))) - (for/list ([line (in-list lines)]) - (match-define (vector u1 v1 u2 v2) line) - (list (vector (unsolve-t xa xb u1) (unsolve-t ya yb v1) z) - (vector (unsolve-t xa xb u2) (unsolve-t ya yb v2) z))))])) + (define lines (do-heights->lines z z1 z2 z3 z4)) + (for/list ([line (in-list lines)]) + (match-define (vector u1 v1 u2 v2) line) + (list (vector (unsolve-t xa xb u1) (unsolve-t ya yb v1) z) + (vector (unsolve-t xa xb u2) (unsolve-t ya yb v2) z)))])) ;; ============================================================================= ;; Isoband marching squares: polygonizes contour between two isoline values -(define ((rotate-facet f) za zb z1 z2 z3 z4) - (map (λ (poly) (map (λ (v) (unrotate-vec v)) poly)) - (f za zb z2 z3 z4 z1))) +(define-syntax-rule (unrotate-vec v) + (match-let ([(vector x y z) v]) + (vector (- 1 y) x z))) -(define ((mirror-x-facet f) za zb z1 z2 z3 z4) - (map (λ (poly) (map (λ (v) (mirror-x-vec v)) poly)) - (f za zb z2 z1 z4 z3))) +(define-syntax-rule (unrotate-vec2 v) + (match-let ([(vector x y z) v]) + (vector (- 1 x) (- 1 y) z))) -(define ((mirror-y-facet f) za zb z1 z2 z3 z4) - (map (λ (poly) (map (λ (v) (mirror-y-vec v)) poly)) - (f za zb z4 z3 z2 z1))) +(define-syntax-rule (unrotate-vec3 v) + (match-let ([(vector x y z) v]) + (vector y (- 1 x) z))) + +(define-syntax-rule (rotate-facet f) + (λ (za zb z1 z2 z3 z4) + (map (λ (vs) (map (λ (v) (unrotate-vec v)) vs)) + (f za zb z2 z3 z4 z1)))) + +(define-syntax-rule (rotate-facet2 f) + (λ (za zb z1 z2 z3 z4) + (map (λ (vs) (map (λ (v) (unrotate-vec2 v)) vs)) + (f za zb z3 z4 z1 z2)))) + +(define-syntax-rule (rotate-facet3 f) + (λ (za zb z1 z2 z3 z4) + (map (λ (vs) (map (λ (v) (unrotate-vec3 v)) vs)) + (f za zb z4 z1 z2 z3)))) ;; ----------------------------------------------------------------------------- ;; all corners same @@ -210,233 +200,233 @@ above. ;; single triangle (define (polys1000 za zb z1 z2 z3 z4) - (list (list (vector 0.0 0.0 z1) - (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector 0.0 (unsafe-solve-t za z1 z4) za)))) + (list (list (vector 0 0 z1) + (vector (solve-t za z1 z2) 0 za) + (vector 0 (solve-t za z1 z4) za)))) (define polys0100 (rotate-facet polys1000)) -(define polys0010 (rotate-facet polys0100)) -(define polys0001 (rotate-facet polys0010)) +(define polys0010 (rotate-facet2 polys1000)) +(define polys0001 (rotate-facet3 polys1000)) (define (polys1222 za zb z1 z2 z3 z4) - (list (list (vector 0.0 0.0 z1) - (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb)))) + (list (list (vector 0 0 z1) + (vector (solve-t zb z1 z2) 0 zb) + (vector 0 (solve-t zb z1 z4) zb)))) (define polys2122 (rotate-facet polys1222)) -(define polys2212 (rotate-facet polys2122)) -(define polys2221 (rotate-facet polys2212)) +(define polys2212 (rotate-facet2 polys1222)) +(define polys2221 (rotate-facet3 polys1222)) ;; ----------------------------------------------------------------------------- ;; single trapezoid (define (polys2000 za zb z1 z2 z3 z4) - (list (list (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector 0.0 (unsafe-solve-t za z1 z4) za) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb)))) + (list (list (vector (solve-t zb z1 z2) 0 zb) + (vector (solve-t za z1 z2) 0 za) + (vector 0 (solve-t za z1 z4) za) + (vector 0 (solve-t zb z1 z4) zb)))) (define polys0200 (rotate-facet polys2000)) -(define polys0020 (rotate-facet polys0200)) -(define polys0002 (rotate-facet polys0020)) +(define polys0020 (rotate-facet2 polys2000)) +(define polys0002 (rotate-facet3 polys2000)) (define (polys0222 za zb z1 z2 z3 z4) - (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb) - (vector 0.0 (unsafe-solve-t za z1 z4) za)))) + (list (list (vector (solve-t za z1 z2) 0 za) + (vector (solve-t zb z1 z2) 0 zb) + (vector 0 (solve-t zb z1 z4) zb) + (vector 0 (solve-t za z1 z4) za)))) (define polys2022 (rotate-facet polys0222)) -(define polys2202 (rotate-facet polys2022)) -(define polys2220 (rotate-facet polys2202)) +(define polys2202 (rotate-facet2 polys0222)) +(define polys2220 (rotate-facet3 polys0222)) ;; ----------------------------------------------------------------------------- ;; single rectangle (define (polys1100 za zb z1 z2 z3 z4) - (list (list (vector 0.0 0.0 z1) - (vector 1.0 0.0 z2) - (vector 1.0 (unsafe-solve-t za z2 z3) za) - (vector 0.0 (unsafe-solve-t za z1 z4) za)))) + (list (list (vector 0 0 z1) + (vector 1 0 z2) + (vector 1 (solve-t za z2 z3) za) + (vector 0 (solve-t za z1 z4) za)))) (define polys0110 (rotate-facet polys1100)) -(define polys0011 (rotate-facet polys0110)) -(define polys1001 (rotate-facet polys0011)) +(define polys0011 (rotate-facet2 polys1100)) +(define polys1001 (rotate-facet3 polys1100)) (define (polys1122 za zb z1 z2 z3 z4) - (list (list (vector 0.0 0.0 z1) - (vector 1.0 0.0 z2) - (vector 1.0 (unsafe-solve-t zb z2 z3) zb) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb)))) + (list (list (vector 0 0 z1) + (vector 1 0 z2) + (vector 1 (solve-t zb z2 z3) zb) + (vector 0 (solve-t zb z1 z4) zb)))) (define polys2112 (rotate-facet polys1122)) -(define polys2211 (rotate-facet polys2112)) -(define polys1221 (rotate-facet polys2211)) +(define polys2211 (rotate-facet2 polys1122)) +(define polys1221 (rotate-facet3 polys1122)) (define (polys0022 za zb z1 z2 z3 z4) - (list (list (vector 0.0 (unsafe-solve-t za z1 z4) za) - (vector 1.0 (unsafe-solve-t za z2 z3) za) - (vector 1.0 (unsafe-solve-t zb z2 z3) zb) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb)))) + (list (list (vector 0 (solve-t za z1 z4) za) + (vector 1 (solve-t za z2 z3) za) + (vector 1 (solve-t zb z2 z3) zb) + (vector 0 (solve-t zb z1 z4) zb)))) (define polys2002 (rotate-facet polys0022)) -(define polys2200 (rotate-facet polys2002)) -(define polys0220 (rotate-facet polys2200)) +(define polys2200 (rotate-facet2 polys0022)) +(define polys0220 (rotate-facet3 polys0022)) ;; ----------------------------------------------------------------------------- ;; single pentagon (define (polys0111 za zb z1 z2 z3 z4) - (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector 1.0 0.0 z2) - (vector 1.0 1.0 z3) - (vector 0.0 1.0 z4) - (vector 0.0 (unsafe-solve-t za z1 z4) za)))) + (list (list (vector (solve-t za z1 z2) 0 za) + (vector 1 0 z2) + (vector 1 1 z3) + (vector 0 1 z4) + (vector 0 (solve-t za z1 z4) za)))) (define polys1011 (rotate-facet polys0111)) -(define polys1101 (rotate-facet polys1011)) -(define polys1110 (rotate-facet polys1101)) +(define polys1101 (rotate-facet2 polys0111)) +(define polys1110 (rotate-facet3 polys0111)) (define (polys2111 za zb z1 z2 z3 z4) - (list (list (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector 1.0 0.0 z2) - (vector 1.0 1.0 z3) - (vector 0.0 1.0 z4) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb)))) + (list (list (vector (solve-t zb z1 z2) 0 zb) + (vector 1 0 z2) + (vector 1 1 z3) + (vector 0 1 z4) + (vector 0 (solve-t zb z1 z4) zb)))) (define polys1211 (rotate-facet polys2111)) -(define polys1121 (rotate-facet polys1211)) -(define polys1112 (rotate-facet polys1121)) +(define polys1121 (rotate-facet2 polys2111)) +(define polys1112 (rotate-facet3 polys2111)) (define (polys1002 za zb z1 z2 z3 z4) - (list (list (vector 0.0 0.0 z1) - (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector (unsafe-solve-t za z4 z3) 1.0 za) - (vector (unsafe-solve-t zb z4 z3) 1.0 zb) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb)))) + (list (list (vector 0 0 z1) + (vector (solve-t za z1 z2) 0 za) + (vector (solve-t za z4 z3) 1 za) + (vector (solve-t zb z4 z3) 1 zb) + (vector 0 (solve-t zb z1 z4) zb)))) (define polys2100 (rotate-facet polys1002)) -(define polys0210 (rotate-facet polys2100)) -(define polys0021 (rotate-facet polys0210)) +(define polys0210 (rotate-facet2 polys1002)) +(define polys0021 (rotate-facet3 polys1002)) (define (polys1220 za zb z1 z2 z3 z4) - (list (list (vector 0.0 0.0 z1) - (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector (unsafe-solve-t zb z4 z3) 1.0 zb) - (vector (unsafe-solve-t za z4 z3) 1.0 za) - (vector 0.0 (unsafe-solve-t za z1 z4) za)))) + (list (list (vector 0 0 z1) + (vector (solve-t zb z1 z2) 0 zb) + (vector (solve-t zb z4 z3) 1 zb) + (vector (solve-t za z4 z3) 1 za) + (vector 0 (solve-t za z1 z4) za)))) (define polys0122 (rotate-facet polys1220)) -(define polys2012 (rotate-facet polys0122)) -(define polys2201 (rotate-facet polys2012)) +(define polys2012 (rotate-facet2 polys1220)) +(define polys2201 (rotate-facet3 polys1220)) (define (polys1200 za zb z1 z2 z3 z4) - (list (list (vector 0.0 0.0 z1) - (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector 1.0 (unsafe-solve-t zb z2 z3) zb) - (vector 1.0 (unsafe-solve-t za z2 z3) za) - (vector 0.0 (unsafe-solve-t za z1 z4) za)))) + (list (list (vector 0 0 z1) + (vector (solve-t zb z1 z2) 0 zb) + (vector 1 (solve-t zb z2 z3) zb) + (vector 1 (solve-t za z2 z3) za) + (vector 0 (solve-t za z1 z4) za)))) (define polys0120 (rotate-facet polys1200)) -(define polys0012 (rotate-facet polys0120)) -(define polys2001 (rotate-facet polys0012)) +(define polys0012 (rotate-facet2 polys1200)) +(define polys2001 (rotate-facet3 polys1200)) (define (polys1022 za zb z1 z2 z3 z4) - (list (list (vector 0.0 0.0 z1) - (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector 1.0 (unsafe-solve-t za z2 z3) za) - (vector 1.0 (unsafe-solve-t zb z2 z3) zb) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb)))) + (list (list (vector 0 0 z1) + (vector (solve-t za z1 z2) 0 za) + (vector 1 (solve-t za z2 z3) za) + (vector 1 (solve-t zb z2 z3) zb) + (vector 0 (solve-t zb z1 z4) zb)))) (define polys2102 (rotate-facet polys1022)) -(define polys2210 (rotate-facet polys2102)) -(define polys0221 (rotate-facet polys2210)) +(define polys2210 (rotate-facet2 polys1022)) +(define polys0221 (rotate-facet3 polys1022)) ;; ----------------------------------------------------------------------------- ;; single hexagon (define (polys0112 za zb z1 z2 z3 z4) - (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector 1.0 0.0 z2) - (vector 1.0 1.0 z3) - (vector (unsafe-solve-t zb z4 z3) 1.0 zb) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb) - (vector 0.0 (unsafe-solve-t za z1 z4) za)))) + (list (list (vector (solve-t za z1 z2) 0 za) + (vector 1 0 z2) + (vector 1 1 z3) + (vector (solve-t zb z4 z3) 1 zb) + (vector 0 (solve-t zb z1 z4) zb) + (vector 0 (solve-t za z1 z4) za)))) (define polys2011 (rotate-facet polys0112)) -(define polys1201 (rotate-facet polys2011)) -(define polys1120 (rotate-facet polys1201)) +(define polys1201 (rotate-facet2 polys0112)) +(define polys1120 (rotate-facet3 polys0112)) (define (polys2110 za zb z1 z2 z3 z4) - (list (list (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector 1.0 0.0 z2) - (vector 1.0 1.0 z3) - (vector (unsafe-solve-t za z4 z3) 1.0 za) - (vector 0.0 (unsafe-solve-t za z1 z4) za) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb)))) + (list (list (vector (solve-t zb z1 z2) 0 zb) + (vector 1 0 z2) + (vector 1 1 z3) + (vector (solve-t za z4 z3) 1 za) + (vector 0 (solve-t za z1 z4) za) + (vector 0 (solve-t zb z1 z4) zb)))) (define polys0211 (rotate-facet polys2110)) -(define polys1021 (rotate-facet polys0211)) -(define polys1102 (rotate-facet polys1021)) +(define polys1021 (rotate-facet2 polys2110)) +(define polys1102 (rotate-facet3 polys2110)) (define (polys0121 za zb z1 z2 z3 z4) - (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector 1.0 0.0 z2) - (vector 1.0 (unsafe-solve-t zb z2 z3) zb) - (vector (unsafe-solve-t zb z4 z3) 1.0 zb) - (vector 0.0 1.0 z4) - (vector 0.0 (unsafe-solve-t za z1 z4) za)))) + (list (list (vector (solve-t za z1 z2) 0 za) + (vector 1 0 z2) + (vector 1 (solve-t zb z2 z3) zb) + (vector (solve-t zb z4 z3) 1 zb) + (vector 0 1 z4) + (vector 0 (solve-t za z1 z4) za)))) (define polys1012 (rotate-facet polys0121)) -(define polys2101 (rotate-facet polys1012)) -(define polys1210 (rotate-facet polys2101)) +(define polys2101 (rotate-facet2 polys0121)) +(define polys1210 (rotate-facet3 polys0121)) ;; ----------------------------------------------------------------------------- ;; 6-sided saddle (define (polys10100 za zb z1 z2 z3 z4) - (list (list (vector 0.0 0.0 z1) - (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector 0.0 (unsafe-solve-t za z1 z4) za)) - (list (vector 1.0 1.0 z3) - (vector (unsafe-solve-t za z4 z3) 1.0 za) - (vector 1.0 (unsafe-solve-t za z2 z3) za)))) + (list (list (vector 0 0 z1) + (vector (solve-t za z1 z2) 0 za) + (vector 0 (solve-t za z1 z4) za)) + (list (vector 1 1 z3) + (vector (solve-t za z4 z3) 1 za) + (vector 1 (solve-t za z2 z3) za)))) (define (polys10101 za zb z1 z2 z3 z4) - (list (list (vector 0.0 0.0 z1) - (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector 1.0 (unsafe-solve-t za z2 z3) za) - (vector 1.0 1.0 z3) - (vector (unsafe-solve-t za z4 z3) 1.0 za) - (vector 0.0 (unsafe-solve-t za z1 z4) za)))) + (list (list (vector 0 0 z1) + (vector (solve-t za z1 z2) 0 za) + (vector 1 (solve-t za z2 z3) za) + (vector 1 1 z3) + (vector (solve-t za z4 z3) 1 za) + (vector 0 (solve-t za z1 z4) za)))) (define (polys1010 za zb z1 z2 z3 z4) - (define z5 (unsafe-flavg4 z1 z2 z3 z4)) - (cond [(z5 . unsafe-fl< . za) (polys10100 za zb z1 z2 z3 z4)] + (define z5 (/ (+ z1 z2 z3 z4) 4)) + (cond [(z5 . < . za) (polys10100 za zb z1 z2 z3 z4)] ; (z5 . >= . zb) is impossible [else (polys10101 za zb z1 z2 z3 z4)])) (define polys0101 (rotate-facet polys1010)) (define (polys1212-2 za zb z1 z2 z3 z4) - (list (list (vector 0.0 0.0 z1) - (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb)) - (list (vector 1.0 1.0 z3) - (vector (unsafe-solve-t zb z4 z3) 1.0 zb) - (vector 1.0 (unsafe-solve-t zb z2 z3) zb)))) + (list (list (vector 0 0 z1) + (vector (solve-t zb z1 z2) 0 zb) + (vector 0 (solve-t zb z1 z4) zb)) + (list (vector 1 1 z3) + (vector (solve-t zb z4 z3) 1 zb) + (vector 1 (solve-t zb z2 z3) zb)))) (define (polys1212-1 za zb z1 z2 z3 z4) - (list (list (vector 0.0 0.0 z1) - (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector 1.0 (unsafe-solve-t zb z2 z3) zb) - (vector 1.0 1.0 z3) - (vector (unsafe-solve-t zb z4 z3) 1.0 zb) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb)))) + (list (list (vector 0 0 z1) + (vector (solve-t zb z1 z2) 0 zb) + (vector 1 (solve-t zb z2 z3) zb) + (vector 1 1 z3) + (vector (solve-t zb z4 z3) 1 zb) + (vector 0 (solve-t zb z1 z4) zb)))) (define (polys1212 za zb z1 z2 z3 z4) - (define z5 (unsafe-flavg4 z1 z2 z3 z4)) - (cond [(z5 . unsafe-fl>= . zb) (polys1212-2 za zb z1 z2 z3 z4)] + (define z5 (/ (+ z1 z2 z3 z4) 4)) + (cond [(z5 . >= . zb) (polys1212-2 za zb z1 z2 z3 z4)] ; (z5 . < . za) is impossible [else (polys1212-1 za zb z1 z2 z3 z4)])) @@ -446,98 +436,98 @@ above. ;; 7-sided saddle (define (polys0212-1 za zb z1 z2 z3 z4) - (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector 1.0 (unsafe-solve-t zb z2 z3) zb) - (vector 1.0 1.0 z3) - (vector (unsafe-solve-t zb z4 z3) 1.0 zb) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb) - (vector 0.0 (unsafe-solve-t za z1 z4) za)))) + (list (list (vector (solve-t za z1 z2) 0 za) + (vector (solve-t zb z1 z2) 0 zb) + (vector 1 (solve-t zb z2 z3) zb) + (vector 1 1 z3) + (vector (solve-t zb z4 z3) 1 zb) + (vector 0 (solve-t zb z1 z4) zb) + (vector 0 (solve-t za z1 z4) za)))) (define (polys0212-2 za zb z1 z2 z3 z4) - (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb) - (vector 0.0 (unsafe-solve-t za z1 z4) za)) - (list (vector 1.0 (unsafe-solve-t zb z2 z3) zb) - (vector 1.0 1.0 z3) - (vector (unsafe-solve-t zb z4 z3) 1.0 zb)))) + (list (list (vector (solve-t za z1 z2) 0 za) + (vector (solve-t zb z1 z2) 0 zb) + (vector 0 (solve-t zb z1 z4) zb) + (vector 0 (solve-t za z1 z4) za)) + (list (vector 1 (solve-t zb z2 z3) zb) + (vector 1 1 z3) + (vector (solve-t zb z4 z3) 1 zb)))) (define (polys0212 za zb z1 z2 z3 z4) - (define z5 (unsafe-flavg4 z1 z2 z3 z4)) - (cond [(z5 . unsafe-fl< . zb) (polys0212-1 za zb z1 z2 z3 z4)] + (define z5 (/ (+ z1 z2 z3 z4) 4)) + (cond [(z5 . < . zb) (polys0212-1 za zb z1 z2 z3 z4)] ; handling (z5 . < . za) separately results in a non-convex polygon [else (polys0212-2 za zb z1 z2 z3 z4)])) (define polys2021 (rotate-facet polys0212)) -(define polys1202 (rotate-facet polys2021)) -(define polys2120 (rotate-facet polys1202)) +(define polys1202 (rotate-facet2 polys0212)) +(define polys2120 (rotate-facet3 polys0212)) (define (polys2010-1 za zb z1 z2 z3 z4) - (list (list (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector 1.0 (unsafe-solve-t za z2 z3) za) - (vector 1.0 1.0 z3) - (vector (unsafe-solve-t za z4 z3) 1.0 za) - (vector 0.0 (unsafe-solve-t za z1 z4) za) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb)))) + (list (list (vector (solve-t zb z1 z2) 0 zb) + (vector (solve-t za z1 z2) 0 za) + (vector 1 (solve-t za z2 z3) za) + (vector 1 1 z3) + (vector (solve-t za z4 z3) 1 za) + (vector 0 (solve-t za z1 z4) za) + (vector 0 (solve-t zb z1 z4) zb)))) (define (polys2010-0 za zb z1 z2 z3 z4) - (list (list (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector 0.0 (unsafe-solve-t za z1 z4) za) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb)) - (list (vector 1.0 (unsafe-solve-t za z2 z3) za) - (vector 1.0 1.0 z3) - (vector (unsafe-solve-t za z4 z3) 1.0 za)))) + (list (list (vector (solve-t zb z1 z2) 0 zb) + (vector (solve-t za z1 z2) 0 za) + (vector 0 (solve-t za z1 z4) za) + (vector 0 (solve-t zb z1 z4) zb)) + (list (vector 1 (solve-t za z2 z3) za) + (vector 1 1 z3) + (vector (solve-t za z4 z3) 1 za)))) (define (polys2010 za zb z1 z2 z3 z4) - (define z5 (unsafe-flavg4 z1 z2 z3 z4)) - (cond [(z5 . unsafe-fl>= . za) (polys2010-1 za zb z1 z2 z3 z4)] + (define z5 (/ (+ z1 z2 z3 z4) 4)) + (cond [(z5 . >= . za) (polys2010-1 za zb z1 z2 z3 z4)] ; handling (z5 . >= . zb) separately results in a non-convex polygon [else (polys2010-0 za zb z1 z2 z3 z4)])) (define polys0201 (rotate-facet polys2010)) -(define polys1020 (rotate-facet polys0201)) -(define polys0102 (rotate-facet polys1020)) +(define polys1020 (rotate-facet2 polys2010)) +(define polys0102 (rotate-facet3 polys2010)) ;; ----------------------------------------------------------------------------- ;; 8-sided saddle (define (polys0202-0 za zb z1 z2 z3 z4) - (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector 1.0 (unsafe-solve-t zb z2 z3) zb) - (vector 1.0 (unsafe-solve-t za z2 z3) za)) - (list (vector 0.0 (unsafe-solve-t za z1 z4) za) - (vector (unsafe-solve-t za z4 z3) 1.0 za) - (vector (unsafe-solve-t zb z4 z3) 1.0 zb) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb)))) + (list (list (vector (solve-t za z1 z2) 0 za) + (vector (solve-t zb z1 z2) 0 zb) + (vector 1 (solve-t zb z2 z3) zb) + (vector 1 (solve-t za z2 z3) za)) + (list (vector 0 (solve-t za z1 z4) za) + (vector (solve-t za z4 z3) 1 za) + (vector (solve-t zb z4 z3) 1 zb) + (vector 0 (solve-t zb z1 z4) zb)))) (define (polys0202-1 za zb z1 z2 z3 z4) - (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector 1.0 (unsafe-solve-t zb z2 z3) zb) - (vector 1.0 (unsafe-solve-t za z2 z3) za) - (vector (unsafe-solve-t za z4 z3) 1.0 za) - (vector (unsafe-solve-t zb z4 z3) 1.0 zb) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb) - (vector 0.0 (unsafe-solve-t za z1 z4) za)))) + (list (list (vector (solve-t za z1 z2) 0 za) + (vector (solve-t zb z1 z2) 0 zb) + (vector 1 (solve-t zb z2 z3) zb) + (vector 1 (solve-t za z2 z3) za) + (vector (solve-t za z4 z3) 1 za) + (vector (solve-t zb z4 z3) 1 zb) + (vector 0 (solve-t zb z1 z4) zb) + (vector 0 (solve-t za z1 z4) za)))) (define (polys0202-2 za zb z1 z2 z3 z4) - (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za) - (vector (unsafe-solve-t zb z1 z2) 0.0 zb) - (vector 0.0 (unsafe-solve-t zb z1 z4) zb) - (vector 0.0 (unsafe-solve-t za z1 z4) za)) - (list (vector 1.0 (unsafe-solve-t zb z2 z3) zb) - (vector 1.0 (unsafe-solve-t za z2 z3) za) - (vector (unsafe-solve-t za z4 z3) 1.0 za) - (vector (unsafe-solve-t zb z4 z3) 1.0 zb)))) + (list (list (vector (solve-t za z1 z2) 0 za) + (vector (solve-t zb z1 z2) 0 zb) + (vector 0 (solve-t zb z1 z4) zb) + (vector 0 (solve-t za z1 z4) za)) + (list (vector 1 (solve-t zb z2 z3) zb) + (vector 1 (solve-t za z2 z3) za) + (vector (solve-t za z4 z3) 1 za) + (vector (solve-t zb z4 z3) 1 zb)))) (define (polys0202 za zb z1 z2 z3 z4) - (define z5 (unsafe-flavg4 z1 z2 z3 z4)) - (cond [(z5 . unsafe-fl< . za) (polys0202-0 za zb z1 z2 z3 z4)] - [(z5 . unsafe-fl< . zb) (polys0202-1 za zb z1 z2 z3 z4)] + (define z5 (/ (+ z1 z2 z3 z4) 4)) + (cond [(z5 . < . za) (polys0202-0 za zb z1 z2 z3 z4)] + [(z5 . < . zb) (polys0202-1 za zb z1 z2 z3 z4)] [else (polys0202-2 za zb z1 z2 z3 z4)])) (define polys2020 (rotate-facet polys0202)) @@ -583,11 +573,11 @@ above. polys2210 polys2211 polys2212 polys2220 polys2221 polys2222)) -(define (unsafe-heights->polys za zb z1 z2 z3 z4) - (define t1 (if (z1 . unsafe-fl< . za) 0 (if (z1 . unsafe-fl<= . zb) 1 2))) - (define t2 (if (z2 . unsafe-fl< . za) 0 (if (z2 . unsafe-fl<= . zb) 1 2))) - (define t3 (if (z3 . unsafe-fl< . za) 0 (if (z3 . unsafe-fl<= . zb) 1 2))) - (define t4 (if (z4 . unsafe-fl< . za) 0 (if (z4 . unsafe-fl<= . zb) 1 2))) +(define (do-heights->polys za zb z1 z2 z3 z4) + (define t1 (if (z1 . < . za) 0 (if (z1 . <= . zb) 1 2))) + (define t2 (if (z2 . < . za) 0 (if (z2 . <= . zb) 1 2))) + (define t3 (if (z3 . < . za) 0 (if (z3 . <= . zb) 1 2))) + (define t4 (if (z4 . < . za) 0 (if (z4 . <= . zb) 1 2))) (define facet-num (unsafe-fx+ (unsafe-fx+ (unsafe-fx+ (unsafe-fx* (unsafe-fx* (unsafe-fx* t1 3) 3) 3) (unsafe-fx* (unsafe-fx* t2 3) 3)) @@ -600,36 +590,16 @@ above. [za real?] [zb real?] [z1 real?] [z2 real?] [z3 real?] [z4 real?] ) (listof (listof (vector/c real? real? real?))) - (cond [(all inexact-real? xa xb ya yb za zb z1 z2 z3 z4) - (define polys (unsafe-heights->polys za zb z1 z2 z3 z4)) - (for/list ([poly (in-list polys)]) - (cond [(eq? poly 'full) (list (vector xa ya z1) (vector xb ya z2) - (vector xb yb z3) (vector xa yb z4))] - [else (for/list ([uv (in-list poly)]) - (match-define (vector u v z) uv) - (vector (unsafe-unsolve-t xa xb u) (unsafe-unsolve-t ya yb v) z))]))] - [(find-failure-index real? xa xb ya yb za zb z1 z2 z3 z4) - => (λ (i) (raise-type-error 'heights->polys "real number" i xa xb ya yb za zb z1 z2 z3 z4))] - [(= za zb z1 z2 z3 z4) (list (list (vector xa ya z1) (vector xb ya z2) - (vector xb yb z3) (vector xa yb z4)))] - [else - (let-map - (za zb z1 z2 z3 z4) inexact->exact - (define z-min (min za zb z1 z2 z3 z4)) - (define z-max (max za zb z1 z2 z3 z4)) - (define z-scale (- z-max z-min)) - (define polys - (unsafe-heights->polys (exact->inexact (/ (- za z-min) z-scale)) - (exact->inexact (/ (- zb z-min) z-scale)) - (exact->inexact (/ (- z1 z-min) z-scale)) - (exact->inexact (/ (- z2 z-min) z-scale)) - (exact->inexact (/ (- z3 z-min) z-scale)) - (exact->inexact (/ (- z4 z-min) z-scale)))) - (for/list ([poly (in-list polys)]) - (cond [(eq? poly 'full) (list (vector xa ya z1) (vector xb ya z2) - (vector xb yb z3) (vector xa yb z4))] - [else (for/list ([uv (in-list poly)]) - (match-define (vector u v z) uv) - (vector (unsolve-t xa xb u) - (unsolve-t ya yb v) - (+ z-min (* (inexact->exact z) z-scale))))])))])) + (cond + [(and (zb . < . z1) (zb . < . z2) (zb . < . z3) (zb . < . z4)) empty] + [(and (za . > . z1) (za . > . z2) (za . > . z3) (za . > . z4)) empty] + [(= za zb z1 z2 z3 z4) (list (list (vector xa ya z1) (vector xb ya z2) + (vector xb yb z3) (vector xa yb z4)))] + [else + (define polys (do-heights->polys za zb z1 z2 z3 z4)) + (for/list ([poly (in-list polys)]) + (cond [(eq? poly 'full) (list (vector xa ya z1) (vector xb ya z2) + (vector xb yb z3) (vector xa yb z4))] + [else (for/list ([uvz (in-list poly)]) + (match-define (vector u v z) uvz) + (vector (unsolve-t xa xb u) (unsolve-t ya yb v) z))]))])) diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-utils.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-utils.rkt index 6aeace039a..0b02734e9e 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-utils.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-utils.rkt @@ -10,6 +10,9 @@ ;; if z = zb, this returns 1.0 ;; if z = (za + zb) / 2, this returns 0.5 ;; Intuitively, regard a use (solve-t z za zb) as "the point between za and zb". +(define-syntax-rule (solve-t z za zb) + (/ (- z za) (- zb za))) + (define-syntax-rule (unsafe-solve-t z za zb) (unsafe-fl/ (unsafe-fl- z za) (unsafe-fl- zb za))) @@ -17,8 +20,9 @@ (unsafe-fl+ (unsafe-fl* t zb) (unsafe-fl* (unsafe-fl- 1.0 t) za))) (define-syntax-rule (unsolve-t za zb t) - (let ([t (inexact->exact t)]) - (+ (* t zb) (* (- 1 t) za)))) + (cond [(eq? t 0) za] + [(eq? t 1) zb] + [else (+ (* t zb) (* (- 1 t) za))])) (define-syntax-rule (unsafe-flavg4 z1 z2 z3 z4) (unsafe-fl* 0.25 (unsafe-fl+ (unsafe-fl+ (unsafe-fl+ z1 z2) z3) z4))) diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/common/math.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/common/math.rkt index b1ccb7eef9..1ef5248931 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/common/math.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/common/math.rkt @@ -42,29 +42,6 @@ [xs (cond [(not (andmap flonum? xs)) (raise-type-error 'fldistance "flonums" xs)] [else (unsafe-flsqrt (flsum (λ (x) (unsafe-fl* x x)) xs))])])) -(defproc (flonum-ok-for-range? [x-min rational?] [x-max rational?] - [size exact-positive-integer?]) boolean? - (let/ec return - (let ([x-min (inexact->exact (min x-min x-max))] - [x-max (inexact->exact (max x-min x-max))]) - (define step-size (/ (- x-max x-min) size)) - - (define inexact-x-min (exact->inexact x-min)) - (unless (rational? inexact-x-min) (return #f)) - - (define inexact-x-max (exact->inexact x-max)) - (unless (rational? inexact-x-max) (return #f)) - - (define inexact-x-max-prev (flprev inexact-x-max)) - (unless (rational? inexact-x-max-prev) (return #f)) - - (define inexact-x-min-next (flnext inexact-x-min)) - (unless (rational? inexact-x-min-next) (return #f)) - - (define max-diff (- x-max (inexact->exact inexact-x-max-prev))) - (define min-diff (- (inexact->exact inexact-x-min-next) x-min)) - (and (max-diff . < . step-size) (min-diff . < . step-size))))) - ;; =================================================================================================== ;; Reals diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/common/sample.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/common/sample.rkt index d6d406cf2e..29771dac12 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/common/sample.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/common/sample.rkt @@ -142,11 +142,15 @@ (hash-ref! memo (vector outer-ivl num tx) (λ () (define xs (sample-points outer-ivl inner-ivl num tx)) - (define ys (map* f xs)) - (define rys (filter rational? ys)) - (define-values (y-min y-max) - (cond [(empty? rys) (values #f #f)] - [else (values (apply min* rys) (apply max* rys))])) + (define y-min #f) + (define y-max #f) + (define ys (for/list ([x (in-list xs)]) + (define y (f x)) + (cond [(rational? y) + (unless (and y-min (y . >= . y-min)) (set! y-min y)) + (unless (and y-max (y . <= . y-max)) (set! y-max y)) + (inexact->exact y)] + [else y]))) (sample xs ys (maybe-inexact->exact y-min) (maybe-inexact->exact y-max))))))) @@ -173,11 +177,12 @@ (define z-max #f) (define zss (for/vector #:length (length ys) ([y (in-list ys)]) (for/vector #:length (length xs) ([x (in-list xs)]) - (let ([z (f x y)]) - (when (rational? z) - (unless (and z-min (z . >= . z-min)) (set! z-min z)) - (unless (and z-max (z . <= . z-max)) (set! z-max z))) - z)))) + (define z (f x y)) + (cond [(rational? z) + (unless (and z-min (z . >= . z-min)) (set! z-min z)) + (unless (and z-max (z . <= . z-max)) (set! z-max z)) + (inexact->exact z)] + [else z])))) (2d-sample xs ys zss (maybe-inexact->exact z-min) (maybe-inexact->exact z-max))))))) @@ -208,11 +213,12 @@ (define dsss (for/vector #:length (length zs) ([z (in-list zs)]) (for/vector #:length (length ys) ([y (in-list ys)]) (for/vector #:length (length xs) ([x (in-list xs)]) - (let ([d (f x y z)]) - (when (rational? d) - (unless (and d-min (d . >= . d-min)) (set! d-min d)) - (unless (and d-max (d . <= . d-max)) (set! d-max d))) - d))))) + (define d (f x y z)) + (cond [(rational? d) + (unless (and d-min (d . >= . d-min)) (set! d-min d)) + (unless (and d-max (d . <= . d-max)) (set! d-max d)) + (inexact->exact d)] + [else d]))))) (3d-sample xs ys zs dsss (maybe-inexact->exact d-min) (maybe-inexact->exact d-max))))))) @@ -262,45 +268,3 @@ (values yb ds01 ds11)) (values zb dss1)))) -(defproc (sample-exact->inexact [s sample?]) sample? - (match-define (sample xs ys y-min y-max) s) - (sample (map exact->inexact xs) (map exact->inexact ys) - y-min y-max)) - -(defproc (2d-sample-exact->inexact [s 2d-sample?]) 2d-sample? - (match-define (2d-sample xs ys zss z-min z-max) s) - (2d-sample (map exact->inexact xs) (map exact->inexact ys) - (for/vector #:length (vector-length zss) ([zs (in-vector zss)]) - (for/vector #:length (vector-length zs) ([z (in-vector zs)]) - (exact->inexact z))) - z-min z-max)) - -(defproc (3d-sample-exact->inexact [s 3d-sample?]) 3d-sample? - (match-define (3d-sample xs ys zs dsss d-min d-max) s) - (3d-sample (map exact->inexact xs) (map exact->inexact ys) (map exact->inexact zs) - (for/vector #:length (vector-length dsss) ([dss (in-vector dsss)]) - (for/vector #:length (vector-length dss) ([ds (in-vector dss)]) - (for/vector #:length (vector-length ds) ([d (in-vector ds)]) - (exact->inexact d)))) - d-min d-max)) - -(defproc (flonum-ok-for-2d? [x-min rational?] [x-max rational?] - [y-min rational?] [y-max rational?]) boolean? - (and (flonum-ok-for-range? x-min x-max 10000) - (flonum-ok-for-range? y-min y-max 10000))) - -(defproc (flonum-ok-for-3d? [x-min rational?] [x-max rational?] - [y-min rational?] [y-max rational?] - [z-min rational?] [z-max rational?]) boolean? - (and (flonum-ok-for-range? x-min x-max 10000) - (flonum-ok-for-range? y-min y-max 10000) - (flonum-ok-for-range? z-min z-max 10000))) - -(defproc (flonum-ok-for-4d? [x-min rational?] [x-max rational?] - [y-min rational?] [y-max rational?] - [z-min rational?] [z-max rational?] - [d-min rational?] [d-max rational?]) boolean? - (and (flonum-ok-for-range? x-min x-max 10000) - (flonum-ok-for-range? y-min y-max 10000) - (flonum-ok-for-range? z-min z-max 10000) - (flonum-ok-for-range? d-min d-max 10000))) diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/contracted/math.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/contracted/math.rkt index 8ecafffeef..4cf0305255 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/contracted/math.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/contracted/math.rkt @@ -5,7 +5,7 @@ (require "../common/math.rkt") (provide equal?* ;; Flonums - flblend flsum fldistance (activate-contract-out flonum-ok-for-range?) + flblend flsum fldistance ;; Reals maybe-inexact->exact min* max* blend atan2 sum real-modulo distance diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/contracted/sample.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/contracted/sample.rkt index 98b86146d0..34c466cdb5 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/contracted/sample.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/contracted/sample.rkt @@ -24,12 +24,6 @@ make-function->sampler make-2d-function->sampler make-3d-function->sampler - sample-exact->inexact - 2d-sample-exact->inexact - 3d-sample-exact->inexact - flonum-ok-for-2d? - flonum-ok-for-3d? - flonum-ok-for-4d? ) (contract-out (struct mapped-function ([f (any/c . -> . any/c)] [fmap ((listof any/c) . -> . (listof any/c))]))) diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/clip.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/clip.rkt index 8117af9e79..9ab9a7c737 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/clip.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/clip.rkt @@ -1,57 +1,89 @@ #lang racket/base -;; Small library for clipping points, rectangles, lines and polygons against axial planes. +(require (for-syntax racket/base) + racket/match + racket/list) -(require racket/match racket/list - "../common/utils.rkt") - -(provide point-in-bounds? clip-line clip-lines clip-polygon) +(provide point-in-bounds? + point-inside-plane? + clip-line/bounds + clip-line/plane + clip-lines/bounds + clip-lines/plane + clip-polygon/bounds + clip-polygon/plane) ;; =================================================================================================== -;; Point clipping +;; Basic plane intersection and distance tests + +;; Need to use this instead of `*' to make axial plane tests as fast as hand-written because the +;; optimizer doesn't always figure out (* 0 e) = 0, etc. +(define-syntax (times stx) + (syntax-case stx () + [(_ 0 e) (syntax/loc stx 0)] + [(_ 1 e) (syntax/loc stx e)] + [(_ 0.0 e) (syntax/loc stx 0.0)] + [(_ 1.0 e) (syntax/loc stx (real->double-flonum e))] + [(_ . es) (syntax/loc stx (* . es))])) + +(define-syntax-rule (plane-line-intersect a b c x1-stx y1-stx x2-stx y2-stx) + (let ([x1 x1-stx] [y1 y1-stx] [x2 x2-stx] [y2 y2-stx]) + (let ([dot1 (+ (times a x1) (times b y1))] + [dot2 (+ (times a x2) (times b y2))]) + (let ([denom (- dot1 dot2)]) + (if (zero? denom) + (values x2 y2) + (let ([t (/ (+ dot1 c) denom)]) + (values (+ x1 (* t (- x2 x1))) (+ y1 (* t (- y2 y1)))))))))) + +(define-syntax-rule (plane-point-dist a b c x y) + (+ (times a x) (times b y) c)) + +;; =================================================================================================== +;; Points (define (point-in-bounds? v x-min x-max y-min y-max) (match-define (vector x y) v) (and (<= x-min x x-max) (<= y-min y y-max))) +(define (point-inside-plane? v plane) + (match-define (vector x y) v) + (match-define (vector a b c) plane) + (>= (plane-point-dist a b c x y) 0)) + ;; =================================================================================================== ;; Line clipping -(define (clip-line-x start-in-bounds? x x1 y1 x2 y2) - (let-map - (x x1 y1 x2 y2) inexact->exact - (define t (/ (- x x1) (- x2 x1))) - (cond [start-in-bounds? (values x1 y1 x (+ y1 (* t (- y2 y1))))] - [else (values x (+ y1 (* t (- y2 y1))) x2 y2)]))) +(define (clip-line/plane v1 v2 plane) + (match-define (vector x1 y1) v1) + (match-define (vector x2 y2) v2) + (match-define (vector a b c) plane) + (define v1? (>= (plane-point-dist a b c x1 y1) 0)) + (define v2? (>= (plane-point-dist a b c x2 y2) 0)) + (cond [(and v1? v2?) (values v1 v2)] + [(not (or v1? v2?)) (values #f #f)] + [else + (define-values (x y) (plane-line-intersect a b c x1 y1 x2 y2)) + (if v1? (values v1 (vector x y)) (values (vector x y) v2))])) -(define (clip-line-y start-in-bounds? y x1 y1 x2 y2) - (let-map - (y x1 y1 x2 y2) inexact->exact - (define t (/ (- y y1) (- y2 y1))) - (cond [start-in-bounds? (values x1 y1 (+ x1 (* t (- x2 x1))) y)] - [else (values (+ x1 (* t (- x2 x1))) y x2 y2)]))) +;; --------------------------------------------------------------------------------------------------- +;; Clipping inside axial bounding box -(define (clip-line-x-min x-min x1 y1 x2 y2) - (cond [(and (x1 . >= . x-min) (x2 . < . x-min)) (clip-line-x #t x-min x1 y1 x2 y2)] - [(and (x2 . >= . x-min) (x1 . < . x-min)) (clip-line-x #f x-min x1 y1 x2 y2)] - [else (values x1 y1 x2 y2)])) +(define-syntax-rule (make-clip-line/axis a b gte?) + (λ (val x1 y1 x2 y2) + (define v1? (gte? (plane-point-dist a b (- val) x1 y1) 0)) + (define v2? (gte? (plane-point-dist a b (- val) x2 y2) 0)) + (cond [(or (and v1? v2?) (not (or v1? v2?))) (values x1 y1 x2 y2)] + [else + (define-values (x y) (plane-line-intersect a b (- val) x1 y1 x2 y2)) + (if v1? (values x1 y1 x y) (values x y x2 y2))]))) -(define (clip-line-x-max x-max x1 y1 x2 y2) - (cond [(and (x1 . <= . x-max) (x2 . > . x-max)) (clip-line-x #t x-max x1 y1 x2 y2)] - [(and (x2 . <= . x-max) (x1 . > . x-max)) (clip-line-x #f x-max x1 y1 x2 y2)] - [else (values x1 y1 x2 y2)])) +(define clip-line-x-min (make-clip-line/axis 1 0 >=)) +(define clip-line-x-max (make-clip-line/axis 1 0 <=)) +(define clip-line-y-min (make-clip-line/axis 0 1 >=)) +(define clip-line-y-max (make-clip-line/axis 0 1 <=)) -(define (clip-line-y-min y-min x1 y1 x2 y2) - (cond [(and (y1 . >= . y-min) (y2 . < . y-min)) (clip-line-y #t y-min x1 y1 x2 y2)] - [(and (y2 . >= . y-min) (y1 . < . y-min)) (clip-line-y #f y-min x1 y1 x2 y2)] - [else (values x1 y1 x2 y2)])) - -(define (clip-line-y-max y-max x1 y1 x2 y2) - (cond [(and (y1 . <= . y-max) (y2 . > . y-max)) (clip-line-y #t y-max x1 y1 x2 y2)] - [(and (y2 . <= . y-max) (y1 . > . y-max)) (clip-line-y #f y-max x1 y1 x2 y2)] - [else (values x1 y1 x2 y2)])) - -(define (clip-line v1 v2 x-min x-max y-min y-max) +(define (clip-line/bounds v1 v2 x-min x-max y-min y-max) (let/ec return (match-define (vector x1 y1) v1) (match-define (vector x2 y2) v2) @@ -70,90 +102,127 @@ (values (vector x1 y1) (vector x2 y2))))) ;; =================================================================================================== -;; Polygon clipping +;; Connected lines clipping -(define-syntax-rule (make-clip-polygon in-bounds? clip-line) - (λ (val xs ys) - (cond [(empty? xs) (values empty empty)] - [else - (for/fold ([res-xs empty] [res-ys empty]) ([x1 (in-list (cons (last xs) xs))] - [x2 (in-list xs)] - [y1 (in-list (cons (last ys) ys))] - [y2 (in-list ys)]) - (define v1-in-bounds? (in-bounds? x1 y1 val)) - (define v2-in-bounds? (in-bounds? x2 y2 val)) - (cond [(and v1-in-bounds? v2-in-bounds?) (values (cons x2 res-xs) - (cons y2 res-ys))] - [(and (not v1-in-bounds?) (not v2-in-bounds?)) (values res-xs res-ys)] - [else (let-values ([(x1 y1 x2 y2) (clip-line v1-in-bounds? val x1 y1 x2 y2)]) - (cond [v2-in-bounds? (values (list* x2 x1 res-xs) - (list* y2 y1 res-ys))] - [else (values (cons x2 res-xs) (cons y2 res-ys))]))]))]))) +(define-syntax-rule (make-clip-lines/plane a b c gte?) + (λ (vs) + (cond + [(empty? vs) empty] + [else + (define v1 (first vs)) + (match-define (vector x1 y1) v1) + (define v1? (gte? (plane-point-dist a b c x1 y1) 0)) + + (define init-vss (if v1? (list (list v1)) (list empty))) + + (define-values (vss _x1 _y1 _v1?) + (for/fold ([vss init-vss] [x1 x1] [y1 y1] [v1? v1?]) ([v2 (in-list (rest vs))]) + (define x2 (vector-ref v2 0)) + (define y2 (vector-ref v2 1)) + (define v2? (gte? (plane-point-dist a b c x2 y2) 0)) + (cond [(and v1? v2?) (values (cons (cons v2 (first vss)) (rest vss)) x2 y2 v2?)] + [(not (or v1? v2?)) (values vss x2 y2 v2?)] + [else + (define-values (x y) (plane-line-intersect a b c x1 y1 x2 y2)) + (if v1? + (values (cons (cons (vector x y) (first vss)) (rest vss)) x2 y2 v2?) + (values (cons (list v2 (vector x y)) vss) x2 y2 v2?))]))) + + (filter (compose not empty?) vss)]))) -(define-syntax-rule (x-min-in-bounds? x y x-min) (x . >= . x-min)) -(define-syntax-rule (x-max-in-bounds? x y x-max) (x . <= . x-max)) -(define-syntax-rule (y-min-in-bounds? x y y-min) (y . >= . y-min)) -(define-syntax-rule (y-max-in-bounds? x y y-max) (y . <= . y-max)) +(define (clip-lines/plane vs plane) + (match-define (vector a b c) plane) + ((make-clip-lines/plane a b c >=) vs)) -(define clip-polygon-x-min (make-clip-polygon x-min-in-bounds? clip-line-x)) -(define clip-polygon-x-max (make-clip-polygon x-max-in-bounds? clip-line-x)) -(define clip-polygon-y-min (make-clip-polygon y-min-in-bounds? clip-line-y)) -(define clip-polygon-y-max (make-clip-polygon y-max-in-bounds? clip-line-y)) +;; --------------------------------------------------------------------------------------------------- +;; Clipping inside axial bounding box -(define (clip-polygon vs x-min x-max y-min y-max) +;; Early accept: all endpoints in bounds (or empty vs) +(define (early-accept? vs x-min x-max y-min y-max) + (andmap (λ (v) (and (<= x-min (vector-ref v 0) x-max) + (<= y-min (vector-ref v 1) y-max))) + vs)) + +;; Early reject: all endpoints on the outside of the same plane +(define (early-reject? vs x-min x-max y-min y-max) + (or (andmap (λ (v) ((vector-ref v 0) . < . x-min)) vs) + (andmap (λ (v) ((vector-ref v 0) . > . x-max)) vs) + (andmap (λ (v) ((vector-ref v 1) . < . y-min)) vs) + (andmap (λ (v) ((vector-ref v 1) . > . y-max)) vs))) + +(define-syntax-rule (make-clip-lines/axis a b gte?) + (λ (val vs) + ((make-clip-lines/plane a b (- val) gte?) vs))) + +(define clip-lines-x-min (make-clip-lines/axis 1 0 >=)) +(define clip-lines-x-max (make-clip-lines/axis 1 0 <=)) +(define clip-lines-y-min (make-clip-lines/axis 0 1 >=)) +(define clip-lines-y-max (make-clip-lines/axis 0 1 <=)) + +(define (clip-lines/bounds vs x-min x-max y-min y-max) (let/ec return - ;; early reject: no polygon - (when (empty? vs) (return empty)) - (match-define (list (vector xs ys) ...) vs) - ;; early accept: all endpoints in bounds - (when (and (andmap (λ (x) (<= x-min x x-max)) xs) - (andmap (λ (y) (<= y-min y y-max)) ys)) - (return vs)) - ;; early reject: all endpoints on the outside of the same plane - (when (or (andmap (λ (x) (x . < . x-min)) xs) (andmap (λ (x) (x . > . x-max)) xs) - (andmap (λ (y) (y . < . y-min)) ys) (andmap (λ (y) (y . > . y-max)) ys)) - (return empty)) - (let*-values ([(xs ys) (clip-polygon-x-min x-min xs ys)] - [(xs ys) (clip-polygon-x-max x-max xs ys)] - [(xs ys) (clip-polygon-y-min y-min xs ys)] - [(xs ys) (clip-polygon-y-max y-max xs ys)]) - (map vector xs ys)))) + (when (early-accept? vs x-min x-max y-min y-max) (return (list vs))) + (when (early-reject? vs x-min x-max y-min y-max) (return empty)) + (let* ([vss (clip-lines-x-min x-min vs)] + [_ (when (empty? vss) (return empty))] + [vss (append* (map (λ (vs) (clip-lines-x-max x-max vs)) vss))] + [_ (when (empty? vss) (return empty))] + [vss (append* (map (λ (vs) (clip-lines-y-min y-min vs)) vss))] + [_ (when (empty? vss) (return empty))] + [vss (append* (map (λ (vs) (clip-lines-y-max y-max vs)) vss))]) + vss))) ;; =================================================================================================== -;; Lines clipping +;; Polygon clipping -(define (join-lines lines) - (let loop ([lines lines] [current-line empty]) - (cond [(empty? lines) (list (reverse current-line))] - [(empty? current-line) (loop (rest lines) (reverse (first lines)))] - [else - (match-define v (first current-line)) - (match-define (list v1 v2) (first lines)) - (cond [(equal? v v1) (loop (rest lines) (cons v2 current-line))] - [else (cons (reverse current-line) (loop lines empty))])]))) +(define-syntax-rule (make-clip-polygon/plane a b c gte?) + (λ (vs ls) + (define v1 (last vs)) + (define x1 (vector-ref v1 0)) + (define y1 (vector-ref v1 1)) + (define v1? (gte? (plane-point-dist a b c x1 y1) 0)) + + (define-values (new-vs new-ls _x1 _y1 _v1?) + (for/fold ([vs empty] [ls empty] [x1 x1] [y1 y1] [v1? v1?]) ([v2 (in-list vs)] + [l (in-list ls)]) + (define x2 (vector-ref v2 0)) + (define y2 (vector-ref v2 1)) + (define v2? (gte? (plane-point-dist a b c x2 y2) 0)) + (cond [(and v1? v2?) (values (cons v2 vs) (cons l ls) x2 y2 v2?)] + [(not (or v1? v2?)) (values vs ls x2 y2 v2?)] + [else + (define-values (x y) (plane-line-intersect a b c x1 y1 x2 y2)) + (if v1? + (values (cons (vector x y) vs) (cons l ls) x2 y2 v2?) + (values (list* v2 (vector x y) vs) (list* l #t ls) x2 y2 v2?))]))) + + (values (reverse new-vs) (reverse new-ls)))) -#; -(join-lines - '((#(1 2) #(3 4)) - (#(3 4) #(5 6)) - (#(5 7) #(7 8)) - (#(7 8) #(9 10)))) +(define (clip-polygon/plane vs ls plane) + (match-define (vector a b c) plane) + ((make-clip-polygon/plane a b c >=) vs ls)) -(define (clip-lines vs x-min x-max y-min y-max) +;; --------------------------------------------------------------------------------------------------- +;; Clipping inside axial bounding box + +(define-syntax-rule (make-clip-polygon/axis a b gte?) + (λ (vs ls val) + ((make-clip-polygon/plane a b (- val) gte?) vs ls))) + +(define clip-polygon-x-min (make-clip-polygon/axis 1 0 >=)) +(define clip-polygon-x-max (make-clip-polygon/axis 1 0 <=)) +(define clip-polygon-y-min (make-clip-polygon/axis 0 1 >=)) +(define clip-polygon-y-max (make-clip-polygon/axis 0 1 <=)) + +(define (clip-polygon/bounds vs ls x-min x-max y-min y-max) (let/ec return - ;; early reject: no lines - (when (empty? vs) (return empty)) - (match-define (list (vector xs ys) ...) vs) - ;; early accept: all endpoints in bounds - (when (and (andmap (λ (x) (<= x-min x x-max)) xs) (andmap (λ (y) (<= y-min y y-max)) ys)) - (return (list vs))) - ;; early reject: all endpoints on the outside of the same plane - (when (or (andmap (λ (x) (x . < . x-min)) xs) (andmap (λ (x) (x . > . x-max)) xs) - (andmap (λ (y) (y . < . y-min)) ys) (andmap (λ (y) (y . > . y-max)) ys)) - (return empty)) - (join-lines - (reverse - (for/fold ([res empty]) ([v1 (in-list vs)] [v2 (in-list (rest vs))]) - (let-values ([(v1 v2) (clip-line v1 v2 x-min x-max y-min y-max)]) - (cond [(and v1 v2) (cons (list v1 v2) res)] - [else res]))))))) + (when (early-accept? vs x-min x-max y-min y-max) (return vs ls)) + (when (early-reject? vs x-min x-max y-min y-max) (return empty empty)) + (let*-values ([(vs ls) (clip-polygon-x-min vs ls x-min)] + [(_) (when (empty? vs) (return empty empty))] + [(vs ls) (clip-polygon-x-max vs ls x-max)] + [(_) (when (empty? vs) (return empty empty))] + [(vs ls) (clip-polygon-y-min vs ls y-min)] + [(_) (when (empty? vs) (return empty empty))] + [(vs ls) (clip-polygon-y-max vs ls y-max)]) + (values vs ls)))) diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/contour.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/contour.rkt index 2a0169906c..d35e04673e 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/contour.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/contour.rkt @@ -2,7 +2,7 @@ ;; Renderers for contour lines and contour intervals -(require racket/contract racket/class racket/match racket/list racket/flonum racket/vector racket/math +(require racket/contract racket/class racket/match racket/list racket/vector unstable/latent-contract/defthing plot/utils "../common/utils.rkt") @@ -23,13 +23,10 @@ (when (<= z-min z z-max) (send area put-alpha alpha) (send area put-pen color width style) - (let* ([flonum-ok? (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max)] - [sample (if flonum-ok? (2d-sample-exact->inexact sample) sample)] - [z (if flonum-ok? (exact->inexact z) z)]) - (for-2d-sample - (xa xb ya yb z1 z2 z3 z4) sample - (for/list ([line (in-list (heights->lines xa xb ya yb z z1 z2 z3 z4))]) - (send/apply area put-line (map (λ (v) (vector-take v 2)) line)))))) + (for-2d-sample + (xa xb ya yb z1 z2 z3 z4) sample + (for/list ([line (in-list (heights->lines xa xb ya yb z z1 z2 z3 z4))]) + (send/apply area put-line (map (λ (v) (vector-take v 2)) line))))) (cond [label (line-legend-entry label color width style)] [else empty])) @@ -72,10 +69,7 @@ (let* ([colors (maybe-apply colors zs)] [widths (maybe-apply widths zs)] [styles (maybe-apply styles zs)] - [alphas (maybe-apply alphas zs)] - [flonum-ok? (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max)] - [sample (if flonum-ok? (2d-sample-exact->inexact sample) sample)] - [zs (if flonum-ok? (map exact->inexact zs) zs)]) + [alphas (maybe-apply alphas zs)]) (for ([z (in-list zs)] [color (in-cycle colors)] [width (in-cycle widths)] @@ -138,10 +132,7 @@ (send area put-pen 0 1 'transparent) (let* ([colors (map ->brush-color (maybe-apply colors z-ivls))] [styles (map ->brush-style (maybe-apply styles z-ivls))] - [alphas (maybe-apply alphas z-ivls)] - [flonum-ok? (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max)] - [sample (if flonum-ok? (2d-sample-exact->inexact sample) sample)] - [zs (if flonum-ok? (map exact->inexact zs) zs)]) + [alphas (maybe-apply alphas z-ivls)]) (for ([za (in-list zs)] [zb (in-list (rest zs))] [color (in-cycle colors)] diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/decoration.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/decoration.rkt index 40c6726f15..56435b9da7 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/decoration.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/decoration.rkt @@ -7,9 +7,6 @@ unstable/contract plot/utils "../common/utils.rkt" - "line.rkt" - "interval.rkt" - "point.rkt" "clip.rkt") (provide (all-defined-out)) @@ -98,8 +95,8 @@ (define-values (r-mins r-maxs) (for/lists (r-mins r-maxs) ([θ (in-list θs)]) (define-values (v1 v2) - (clip-line (vector 0 0) (vector (* max-r (cos θ)) (* max-r (sin θ))) - x-min x-max y-min y-max)) + (clip-line/bounds (vector 0 0) (vector (* max-r (cos θ)) (* max-r (sin θ))) + x-min x-max y-min y-max)) (values (if v1 (vmag v1) #f) (if v2 (vmag v2) #f)))) (for/lists (θs r-mins r-maxs) ([θ (in-list θs)] [r-min (in-list r-mins)] [r-max (in-list r-maxs)] diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/interval.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/interval.rkt index 59e4a41cf4..590b41cd80 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/interval.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/interval.rkt @@ -2,7 +2,7 @@ ;; Renderers for intervals between functions. -(require racket/contract racket/class racket/match racket/math racket/list racket/sequence +(require racket/contract racket/class racket/match racket/math racket/list unstable/latent-contract/defthing unstable/contract plot/utils diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/plot-area.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/plot-area.rkt index e86eee5b41..77ff41e35e 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/plot-area.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/plot-area.rkt @@ -1,18 +1,16 @@ #lang racket/base -(require racket/draw racket/class racket/contract racket/match racket/math racket/list racket/string - racket/flonum +(require racket/class racket/match racket/math racket/list racket/flonum + (only-in math/flonum fl) "../common/plot-device.rkt" "../common/ticks.rkt" - "../common/contract.rkt" "../common/math.rkt" "../common/draw.rkt" "../common/axis-transform.rkt" - "../common/sample.rkt" - "../common/legend.rkt" "../common/parameters.rkt" "../common/utils.rkt" - "clip.rkt") + "clip.rkt" + "vector.rkt") (provide (all-defined-out)) @@ -87,45 +85,20 @@ (and (equal? (plot-x-transform) id-transform) (equal? (plot-y-transform) id-transform))) - (define flonum-ok? (flonum-ok-for-2d? x-min x-max y-min y-max)) - (when flonum-ok? - (set! x-min (exact->inexact x-min)) - (set! x-max (exact->inexact x-max)) - (set! y-min (exact->inexact y-min)) - (set! y-max (exact->inexact y-max)) - (set! x-size (exact->inexact x-size)) - (set! y-size (exact->inexact y-size)) - (set! x-mid (exact->inexact x-mid)) - (set! y-mid (exact->inexact y-mid))) - (define plot->view - (if flonum-ok? - (let-map - (x-min y-min x-size y-size) exact->inexact - (if identity-transforms? - (match-lambda - [(vector x y) - (vector (fl/ (fl- (exact->inexact x) x-min) x-size) - (fl/ (fl- (exact->inexact y) y-min) y-size))]) - (match-lambda - [(vector (? rational? x) (? rational? y)) - (vector (fl/ (fl- (exact->inexact (fx x)) x-min) x-size) - (fl/ (fl- (exact->inexact (fy y)) y-min) y-size))] - [(vector x y) - (vector +nan.0 +nan.0)]))) - (if identity-transforms? - (match-lambda - [(vector (? rational? x) (? rational? y)) - (vector (exact->inexact (/ (- (inexact->exact x) x-min) x-size)) - (exact->inexact (/ (- (inexact->exact y) y-min) y-size)))] - [(vector x y) - (vector +nan.0 +nan.0)]) - (match-lambda - [(vector (? rational? x) (? rational? y)) - (vector (exact->inexact (/ (- (inexact->exact (fx x)) x-min) x-size)) - (exact->inexact (/ (- (inexact->exact (fy y)) y-min) y-size)))] - [(vector x y) - (vector +nan.0 +nan.0)])))) + (if identity-transforms? + (match-lambda + [(vector (? rational? x) (? rational? y)) + (vector (fl (/ (- (inexact->exact x) x-min) x-size)) + (fl (/ (- (inexact->exact y) y-min) y-size)))] + [(vector x y) + (vector +nan.0 +nan.0)]) + (match-lambda + [(vector (? rational? x) (? rational? y)) + (vector (fl (/ (- (inexact->exact (fx x)) x-min) x-size)) + (fl (/ (- (inexact->exact (fy y)) y-min) y-size)))] + [(vector x y) + (vector +nan.0 +nan.0)]))) (define view->dc #f) (define (plot->dc v) (view->dc (plot->view v))) @@ -144,7 +117,7 @@ (define area-per-view-x (/ (- area-x-max area-x-min) view-x-size)) (define area-per-view-y (/ (- area-y-max area-y-min) view-y-size)) (let-map - (area-x-min area-y-max area-per-view-x area-per-view-y) exact->inexact + (area-x-min area-y-max area-per-view-x area-per-view-y) fl (λ (v) (match-define (vector x y) v) (vector (fl+ area-x-min (fl* x area-per-view-x)) @@ -503,61 +476,69 @@ ;; Shapes (define/public (put-lines vs) - (for ([vs (vrational-sublists vs)]) - (for ([vs (if clipping? - (in-list (clip-lines vs clip-x-min clip-x-max - clip-y-min clip-y-max)) - (in-value vs))]) - (when (not (empty? vs)) - (let* ([vs (if identity-transforms? vs (subdivide-lines plot->dc vs))] - [vs (map (λ (v) (plot->dc v)) vs)]) - (send pd draw-lines vs)))))) + (for ([vs (in-list (exact-vector2d-sublists vs))]) + (let ([vss (if clipping? + (clip-lines/bounds vs + clip-x-min clip-x-max + clip-y-min clip-y-max) + (list vs))]) + (for ([vs (in-list vss)]) + (unless (empty? vs) + (let* ([vs (if identity-transforms? vs (subdivide-lines plot->dc vs))] + [vs (map plot->dc vs)]) + (send pd draw-lines vs))))))) (define/public (put-line v1 v2) - (when (and (vrational? v1) (vrational? v2)) - (let-values ([(v1 v2) (if clipping? - (clip-line v1 v2 clip-x-min clip-x-max - clip-y-min clip-y-max) - (values v1 v2))]) - (when (and v1 v2) - (if identity-transforms? - (send pd draw-line (plot->dc v1) (plot->dc v2)) - (send pd draw-lines (map (λ (v) (plot->dc v)) - (subdivide-line plot->dc v1 v2)))))))) + (let ([v1 (exact-vector2d v1)] + [v2 (exact-vector2d v2)]) + (when (and v1 v2) + (let-values ([(v1 v2) (if clipping? + (clip-line/bounds v1 v2 + clip-x-min clip-x-max + clip-y-min clip-y-max) + (values v1 v2))]) + (when (and v1 v2) + (if identity-transforms? + (send pd draw-line (plot->dc v1) (plot->dc v2)) + (send pd draw-lines (map plot->dc (subdivide-line plot->dc v1 v2))))))))) (define/public (put-polygon vs) - (when (andmap vrational? vs) - (let* ([vs (if clipping? - (clip-polygon vs clip-x-min clip-x-max - clip-y-min clip-y-max) - vs)]) - (when (not (empty? vs)) - (if identity-transforms? - (send pd draw-polygon (map (λ (v) (plot->dc v)) vs)) - (send pd draw-polygon (map (λ (v) (plot->dc v)) - (subdivide-polygon plot->dc vs)))))))) + (define ls (make-list (length vs) #t)) + (let-values ([(vs ls) (exact-polygon2d vs ls)]) + (unless (empty? vs) + (let-values ([(vs ls) (if clipping? + (clip-polygon/bounds vs ls + clip-x-min clip-x-max + clip-y-min clip-y-max) + (values vs ls))]) + (unless (empty? vs) + (if identity-transforms? + (send pd draw-polygon (map plot->dc vs)) + (send pd draw-polygon (map plot->dc (subdivide-polygon plot->dc vs))))))))) (define/public (put-rect r) (when (rect-rational? r) (match-define (vector (ivl x1 x2) (ivl y1 y2)) r) (put-polygon (list (vector x1 y1) (vector x2 y1) (vector x2 y2) (vector x1 y2))))) - (define/public (put-text str v [anchor 'top-left] [angle 0] [dist 0] - #:outline? [outline? #f]) - (when (and (vrational? v) (in-bounds? v)) - (send pd draw-text str (plot->dc v) anchor angle dist #:outline? outline?))) + (define/public (put-text str v [anchor 'top-left] [angle 0] [dist 0] #:outline? [outline? #f]) + (let ([v (exact-vector2d v)]) + (when (and v (in-bounds? v)) + (send pd draw-text str (plot->dc v) anchor angle dist #:outline? outline?)))) (define/public (put-glyphs vs symbol size) - (send pd draw-glyphs (map (λ (v) (plot->dc v)) - (filter (λ (v) (and (vrational? v) (in-bounds? v))) - vs)) - symbol size)) - + (let ([vs (filter (λ (v) (and v (in-bounds? v))) (map exact-vector2d vs))]) + (unless (empty? vs) + (send pd draw-glyphs (map plot->dc vs) symbol size)))) + (define/public (put-arrow v1 v2) - (when (and (vrational? v1) (vrational? v2) (in-bounds? v1)) - (send pd draw-arrow (plot->dc v1) (plot->dc v2)))) + (let ([v1 (exact-vector2d v1)] + [v2 (exact-vector2d v2)]) + (when (and v1 v2 (in-bounds? v1)) + (send pd draw-arrow (plot->dc v1) (plot->dc v2))))) (define/public (put-tick v r angle) - (when (and (vrational? v) (in-bounds? v)) - (send pd draw-tick (plot->dc v) r angle))) + (let ([v (exact-vector2d v)]) + (when (and v (in-bounds? v)) + (send pd draw-tick (plot->dc v) r angle)))) )) diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/point.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/point.rkt index daf87e5261..2db80ac1ce 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/point.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/point.rkt @@ -2,7 +2,7 @@ ;; Renderers for points and other point-like things. -(require racket/contract racket/class racket/match racket/math racket/list racket/sequence +(require racket/contract racket/class racket/match racket/math racket/list unstable/latent-contract/defthing unstable/contract plot/utils diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/vector.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/vector.rkt new file mode 100644 index 0000000000..aed08755e3 --- /dev/null +++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/vector.rkt @@ -0,0 +1,46 @@ +#lang racket/base + +(require racket/list + (only-in racket/unsafe/ops unsafe-vector-ref)) + +(provide exact-vector2d + exact-vector2d-sublists + exact-polygon2d) + +(define (exact-vector2d v) + (define n (vector-length v)) + (cond [(= n 2) + (define v1 (unsafe-vector-ref v 0)) + (define v2 (unsafe-vector-ref v 1)) + (cond [(and (exact? v1) (exact? v2)) v] + [(and (rational? v1) (rational? v2)) + (vector (inexact->exact v1) (inexact->exact v2))] + [else #f])] + [else #f])) + +(define (sublists vs) + (define vss + (for/fold ([vss (list null)]) ([v (in-list vs)]) + (cond [v (cons (cons v (car vss)) (cdr vss))] + [(null? (car vss)) vss] + [else (cons null vss)]))) + (cond [(null? (car vss)) (cdr vss)] + [else vss])) + +(define (exact-vector2d-sublists vs) + (sublists (map exact-vector2d vs))) + +(define (exact-polygon2d vs ls) + (cond + [(null? vs) (values null null)] + [else + (define-values (new-vs new-ls _) + (for/fold ([vs null] [ls null] [v1 (last vs)]) ([v2 (in-list vs)] + [l (in-list ls)]) + (cond [(equal? v1 v2) (values vs ls v2)] + [else + (define exact-v2 (exact-vector2d v2)) + (if exact-v2 + (values (cons exact-v2 vs) (cons l ls) v2) + (values vs ls v2))]))) + (values (reverse new-vs) (reverse new-ls))])) diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/bsp.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/bsp.rkt index 9c29b76f67..ff136c1e0f 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/bsp.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/bsp.rkt @@ -3,9 +3,7 @@ (require racket/list racket/match racket/promise - racket/math math/flonum - math/statistics "split.rkt") (provide diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/clip.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/clip.rkt index 9406621c41..fe75b7e0b7 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/clip.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/clip.rkt @@ -1,14 +1,45 @@ #lang racket/base -;; Small library for clipping points, lines and polygons against axial planes. +(require (for-syntax racket/base) + racket/match + racket/list) -(require racket/match racket/list racket/unsafe/ops - "../common/utils.rkt") +(provide point-in-bounds? + point-inside-plane? + clip-line/bounds + clip-line/plane + clip-lines/bounds + clip-lines/plane + clip-polygon/bounds + clip-polygon/plane) -(provide point-in-bounds? clip-line clip-polygon clip-lines - clip-polygon-x-min clip-polygon-x-max - clip-polygon-y-min clip-polygon-y-max - clip-polygon-z-min clip-polygon-z-max) +;; =================================================================================================== +;; Basic plane intersection and distance tests + +;; Need to use this instead of `*' to make axial plane tests as fast as hand-written because the +;; optimizer doesn't always figure out (* 0 e) = 0, etc. +(define-syntax (times stx) + (syntax-case stx () + [(_ 0 e) (syntax/loc stx 0)] + [(_ 1 e) (syntax/loc stx e)] + [(_ 0.0 e) (syntax/loc stx 0.0)] + [(_ 1.0 e) (syntax/loc stx (real->double-flonum e))] + [(_ . es) (syntax/loc stx (* . es))])) + +(define-syntax-rule (plane-line-intersect a b c d x1-stx y1-stx z1-stx x2-stx y2-stx z2-stx) + (let ([x1 x1-stx] [y1 y1-stx] [z1 z1-stx] [x2 x2-stx] [y2 y2-stx] [z2 z2-stx]) + (let ([dot1 (+ (times a x1) (times b y1) (times c z1))] + [dot2 (+ (times a x2) (times b y2) (times c z2))]) + (let ([denom (- dot1 dot2)]) + (if (zero? denom) + (values x2 y2 z2) + (let ([t (/ (+ dot1 d) denom)]) + (values (+ x1 (* t (- x2 x1))) + (+ y1 (* t (- y2 y1))) + (+ z1 (* t (- z2 z1)))))))))) + +(define-syntax-rule (plane-point-dist a b c d x y z) + (+ (times a x) (times b y) (times c z) d)) ;; =================================================================================================== ;; Points @@ -17,61 +48,46 @@ (match-define (vector x y z) v) (and (<= x-min x x-max) (<= y-min y y-max) (<= z-min z z-max))) +(define (point-inside-plane? v plane) + (match-define (vector x y z) v) + (match-define (vector a b c d) plane) + (>= (plane-point-dist a b c d x y z) 0)) + ;; =================================================================================================== ;; Line clipping -(define (clip-line-x start-in-bounds? x x1 y1 z1 x2 y2 z2) - (let-map - (x x1 y1 z1 x2 y2 z2) inexact->exact - (define t (/ (- x x1) (- x2 x1))) - (cond [start-in-bounds? (values x1 y1 z1 x (+ y1 (* t (- y2 y1))) (+ z1 (* t (- z2 z1))))] - [else (values x (+ y1 (* t (- y2 y1))) (+ z1 (* t (- z2 z1))) x2 y2 z2)]))) +(define (clip-line/plane v1 v2 plane) + (match-define (vector x1 y1 z1) v1) + (match-define (vector x2 y2 z2) v2) + (match-define (vector a b c d) plane) + (define v1? (>= (plane-point-dist a b c d x1 y1 z1) 0)) + (define v2? (>= (plane-point-dist a b c d x2 y2 z2) 0)) + (cond [(and v1? v2?) (values v1 v2)] + [(not (or v1? v2?)) (values #f #f)] + [else + (define-values (x y z) (plane-line-intersect a b c d x1 y1 z1 x2 y2 z2)) + (if v1? (values v1 (vector x y z)) (values (vector x y z) v2))])) -(define (clip-line-y start-in-bounds? y x1 y1 z1 x2 y2 z2) - (let-map - (y x1 y1 z1 x2 y2 z2) inexact->exact - (define t (/ (- y y1) (- y2 y1))) - (cond [start-in-bounds? (values x1 y1 z1 (+ x1 (* t (- x2 x1))) y (+ z1 (* t (- z2 z1))))] - [else (values (+ x1 (* t (- x2 x1))) y (+ z1 (* t (- z2 z1))) x2 y2 z2)]))) +;; --------------------------------------------------------------------------------------------------- +;; Clipping inside axial bounding box -(define (clip-line-z start-in-bounds? z x1 y1 z1 x2 y2 z2) - (let-map - (z x1 y1 z1 x2 y2 z2) inexact->exact - (define t (/ (- z z1) (- z2 z1))) - (cond [start-in-bounds? (values x1 y1 z1 (+ x1 (* t (- x2 x1))) (+ y1 (* t (- y2 y1))) z)] - [else (values (+ x1 (* t (- x2 x1))) (+ y1 (* t (- y2 y1))) z x2 y2 z2)]))) +(define-syntax-rule (make-clip-line/axis a b c gte?) + (λ (val x1 y1 z1 x2 y2 z2) + (define v1? (gte? (plane-point-dist a b c (- val) x1 y1 z1) 0)) + (define v2? (gte? (plane-point-dist a b c (- val) x2 y2 z2) 0)) + (cond [(or (and v1? v2?) (not (or v1? v2?))) (values x1 y1 z1 x2 y2 z2)] + [else + (define-values (x y z) (plane-line-intersect a b c (- val) x1 y1 z1 x2 y2 z2)) + (if v1? (values x1 y1 z1 x y z) (values x y z x2 y2 z2))]))) -(define (clip-line-x-min x-min x1 y1 z1 x2 y2 z2) - (cond [(and (x1 . >= . x-min) (x2 . < . x-min)) (clip-line-x #t x-min x1 y1 z1 x2 y2 z2)] - [(and (x2 . >= . x-min) (x1 . < . x-min)) (clip-line-x #f x-min x1 y1 z1 x2 y2 z2)] - [else (values x1 y1 z1 x2 y2 z2)])) +(define clip-line-x-min (make-clip-line/axis 1 0 0 >=)) +(define clip-line-x-max (make-clip-line/axis 1 0 0 <=)) +(define clip-line-y-min (make-clip-line/axis 0 1 0 >=)) +(define clip-line-y-max (make-clip-line/axis 0 1 0 <=)) +(define clip-line-z-min (make-clip-line/axis 0 0 1 >=)) +(define clip-line-z-max (make-clip-line/axis 0 0 1 <=)) -(define (clip-line-x-max x-max x1 y1 z1 x2 y2 z2) - (cond [(and (x1 . <= . x-max) (x2 . > . x-max)) (clip-line-x #t x-max x1 y1 z1 x2 y2 z2)] - [(and (x2 . <= . x-max) (x1 . > . x-max)) (clip-line-x #f x-max x1 y1 z1 x2 y2 z2)] - [else (values x1 y1 z1 x2 y2 z2)])) - -(define (clip-line-y-min y-min x1 y1 z1 x2 y2 z2) - (cond [(and (y1 . >= . y-min) (y2 . < . y-min)) (clip-line-y #t y-min x1 y1 z1 x2 y2 z2)] - [(and (y2 . >= . y-min) (y1 . < . y-min)) (clip-line-y #f y-min x1 y1 z1 x2 y2 z2)] - [else (values x1 y1 z1 x2 y2 z2)])) - -(define (clip-line-y-max y-max x1 y1 z1 x2 y2 z2) - (cond [(and (y1 . <= . y-max) (y2 . > . y-max)) (clip-line-y #t y-max x1 y1 z1 x2 y2 z2)] - [(and (y2 . <= . y-max) (y1 . > . y-max)) (clip-line-y #f y-max x1 y1 z1 x2 y2 z2)] - [else (values x1 y1 z1 x2 y2 z2)])) - -(define (clip-line-z-min z-min x1 y1 z1 x2 y2 z2) - (cond [(and (z1 . >= . z-min) (z2 . < . z-min)) (clip-line-z #t z-min x1 y1 z1 x2 y2 z2)] - [(and (z2 . >= . z-min) (z1 . < . z-min)) (clip-line-z #f z-min x1 y1 z1 x2 y2 z2)] - [else (values x1 y1 z1 x2 y2 z2)])) - -(define (clip-line-z-max z-max x1 y1 z1 x2 y2 z2) - (cond [(and (z1 . <= . z-max) (z2 . > . z-max)) (clip-line-z #t z-max x1 y1 z1 x2 y2 z2)] - [(and (z2 . <= . z-max) (z1 . > . z-max)) (clip-line-z #f z-max x1 y1 z1 x2 y2 z2)] - [else (values x1 y1 z1 x2 y2 z2)])) - -(define (clip-line v1 v2 x-min x-max y-min y-max z-min z-max) +(define (clip-line/bounds v1 v2 x-min x-max y-min y-max z-min z-max) (let/ec return (match-define (vector x1 y1 z1) v1) (match-define (vector x2 y2 z2) v2) @@ -92,117 +108,74 @@ [(x1 y1 z1 x2 y2 z2) (clip-line-z-max z-max x1 y1 z1 x2 y2 z2)]) (values (vector x1 y1 z1) (vector x2 y2 z2))))) -;; =================================================================================================== -;; Polygon clipping - -(define-syntax-rule (make-clip-polygon idx test? clip-line) - (λ (val vs ls) - (define-values (new-vs new-ls) - (for/fold ([vs empty] [ls empty]) ([v1 (in-list (cons (last vs) vs))] - [v2 (in-list vs)] - [l (in-list ls)]) - (define v1-in-bounds? (test? (unsafe-vector-ref v1 idx) val)) - (define v2-in-bounds? (test? (unsafe-vector-ref v2 idx) val)) - (cond [(and v1-in-bounds? v2-in-bounds?) - (values (cons v2 vs) (cons l ls))] - [(and (not v1-in-bounds?) (not v2-in-bounds?)) - (values vs ls)] - [else - (match-define (vector x1 y1 z1) v1) - (match-define (vector x2 y2 z2) v2) - (let-values ([(x1 y1 z1 x2 y2 z2) (clip-line v1-in-bounds? val x1 y1 z1 x2 y2 z2)]) - (cond [v2-in-bounds? - (values (list* (vector x2 y2 z2) (vector x1 y1 z1) vs) - (list* l #t ls))] - [else - (values (cons (vector x2 y2 z2) vs) - (cons l ls))]))]))) - (values (reverse new-vs) (reverse new-ls)))) - -(define clip-polygon-x-min (make-clip-polygon 0 >= clip-line-x)) -(define clip-polygon-x-max (make-clip-polygon 0 <= clip-line-x)) -(define clip-polygon-y-min (make-clip-polygon 1 >= clip-line-y)) -(define clip-polygon-y-max (make-clip-polygon 1 <= clip-line-y)) -(define clip-polygon-z-min (make-clip-polygon 2 >= clip-line-z)) -(define clip-polygon-z-max (make-clip-polygon 2 <= clip-line-z)) - -(define (clip-polygon vs ls x-min x-max y-min y-max z-min z-max) - (let/ec return - ;; early reject: no polygon - (when (empty? vs) (return empty empty)) - (match-define (list (vector xs ys zs) ...) vs) - ;; early accept: all endpoints in bounds - (when (and (andmap (λ (x) (<= x-min x x-max)) xs) - (andmap (λ (y) (<= y-min y y-max)) ys) - (andmap (λ (z) (<= z-min z z-max)) zs)) - (return vs ls)) - ;; early reject: all endpoints on the outside of the same plane - (when (or (andmap (λ (x) (x . < . x-min)) xs) (andmap (λ (x) (x . > . x-max)) xs) - (andmap (λ (y) (y . < . y-min)) ys) (andmap (λ (y) (y . > . y-max)) ys) - (andmap (λ (z) (z . < . z-min)) zs) (andmap (λ (z) (z . > . z-max)) zs)) - (return empty empty)) - (let*-values ([(vs ls) (clip-polygon-x-min x-min vs ls)] - [(_) (when (empty? vs) (return empty empty))] - [(vs ls) (clip-polygon-x-max x-max vs ls)] - [(_) (when (empty? vs) (return empty empty))] - [(vs ls) (clip-polygon-y-min y-min vs ls)] - [(_) (when (empty? vs) (return empty empty))] - [(vs ls) (clip-polygon-y-max y-max vs ls)] - [(_) (when (empty? vs) (return empty empty))] - [(vs ls) (clip-polygon-z-min z-min vs ls)] - [(_) (when (empty? vs) (return empty empty))] - [(vs ls) (clip-polygon-z-max z-max vs ls)]) - (values vs ls)))) - ;; =================================================================================================== ;; Connected lines clipping -(define-syntax-rule (make-clip-lines idx test? clip-line) +(define-syntax-rule (make-clip-lines/plane a b c d gte?) + (λ (vs) + (cond + [(empty? vs) empty] + [else + (define v1 (first vs)) + (match-define (vector x1 y1 z1) v1) + (define v1? (gte? (plane-point-dist a b c d x1 y1 z1) 0)) + + (define init-vss (if v1? (list (list v1)) (list empty))) + + (define-values (vss _x1 _y1 _z1 _v1?) + (for/fold ([vss init-vss] [x1 x1] [y1 y1] [z1 z1] [v1? v1?]) ([v2 (in-list (rest vs))]) + (define x2 (vector-ref v2 0)) + (define y2 (vector-ref v2 1)) + (define z2 (vector-ref v2 2)) + (define v2? (gte? (plane-point-dist a b c d x2 y2 z2) 0)) + (cond [(and v1? v2?) (values (cons (cons v2 (first vss)) (rest vss)) x2 y2 z2 v2?)] + [(not (or v1? v2?)) (values vss x2 y2 z2 v2?)] + [else + (define-values (x y z) (plane-line-intersect a b c d x1 y1 z1 x2 y2 z2)) + (if v1? + (values (cons (cons (vector x y z) (first vss)) (rest vss)) x2 y2 z2 v2?) + (values (cons (list v2 (vector x y z)) vss) x2 y2 z2 v2?))]))) + + (filter (compose not empty?) vss)]))) + +(define (clip-lines/plane vs plane) + (match-define (vector a b c d) plane) + ((make-clip-lines/plane a b c d >=) vs)) + +;; --------------------------------------------------------------------------------------------------- +;; Clipping inside axial bounding box + +;; Early accept: all endpoints in bounds (or empty vs) +(define (early-accept? vs x-min x-max y-min y-max z-min z-max) + (andmap (λ (v) (and (<= x-min (vector-ref v 0) x-max) + (<= y-min (vector-ref v 1) y-max) + (<= z-min (vector-ref v 2) z-max))) + vs)) + +;; Early reject: all endpoints on the outside of the same plane +(define (early-reject? vs x-min x-max y-min y-max z-min z-max) + (or (andmap (λ (v) ((vector-ref v 0) . < . x-min)) vs) + (andmap (λ (v) ((vector-ref v 0) . > . x-max)) vs) + (andmap (λ (v) ((vector-ref v 1) . < . y-min)) vs) + (andmap (λ (v) ((vector-ref v 1) . > . y-max)) vs) + (andmap (λ (v) ((vector-ref v 2) . < . z-min)) vs) + (andmap (λ (v) ((vector-ref v 2) . > . z-max)) vs))) + +(define-syntax-rule (make-clip-lines/axis a b c gte?) (λ (val vs) - (define v1 (first vs)) - (define v1? (test? (unsafe-vector-ref v1 idx) val)) - - (define-values (vss _v1 _v1?) - (for/fold ([vss (if v1? (list (list v1)) (list empty))] - [v1 v1] - [v1? v1?] - ) ([v2 (in-list (rest vs))]) - (define v2? (test? (unsafe-vector-ref v2 idx) val)) - (cond [(and v1? v2?) (values (cons (cons v2 (first vss)) (rest vss)) v2 v2?)] - [(not (or v1? v2?)) (values vss v2 v2?)] - [else - (match-define (vector x1 y1 z1) v1) - (match-define (vector x2 y2 z2) v2) - (let-values ([(x1 y1 z1 x2 y2 z2) (clip-line v1? val x1 y1 z1 x2 y2 z2)]) - (cond [v2? - (values (cons (list (vector x2 y2 z2) (vector x1 y1 z1)) vss) v2 v2?)] - [else - (values (cons (cons (vector x2 y2 z2) (first vss)) (rest vss)) v2 v2?)]))]))) - - (filter (compose not empty?) vss))) + ((make-clip-lines/plane a b c (- val) gte?) vs))) -(define clip-lines-x-min (make-clip-lines 0 >= clip-line-x)) -(define clip-lines-x-max (make-clip-lines 0 <= clip-line-x)) -(define clip-lines-y-min (make-clip-lines 1 >= clip-line-y)) -(define clip-lines-y-max (make-clip-lines 1 <= clip-line-y)) -(define clip-lines-z-min (make-clip-lines 2 >= clip-line-z)) -(define clip-lines-z-max (make-clip-lines 2 <= clip-line-z)) +(define clip-lines-x-min (make-clip-lines/axis 1 0 0 >=)) +(define clip-lines-x-max (make-clip-lines/axis 1 0 0 <=)) +(define clip-lines-y-min (make-clip-lines/axis 0 1 0 >=)) +(define clip-lines-y-max (make-clip-lines/axis 0 1 0 <=)) +(define clip-lines-z-min (make-clip-lines/axis 0 0 1 >=)) +(define clip-lines-z-max (make-clip-lines/axis 0 0 1 <=)) -(define (clip-lines vs x-min x-max y-min y-max z-min z-max) +(define (clip-lines/bounds vs x-min x-max y-min y-max z-min z-max) (let/ec return - ;; early reject: no polygon - (when (empty? vs) (return empty)) - (match-define (list (vector xs ys zs) ...) vs) - ;; early accept: all endpoints in bounds - (when (and (andmap (λ (x) (<= x-min x x-max)) xs) - (andmap (λ (y) (<= y-min y y-max)) ys) - (andmap (λ (z) (<= z-min z z-max)) zs)) - (return (list vs))) - ;; early reject: all endpoints on the outside of the same plane - (when (or (andmap (λ (x) (x . < . x-min)) xs) (andmap (λ (x) (x . > . x-max)) xs) - (andmap (λ (y) (y . < . y-min)) ys) (andmap (λ (y) (y . > . y-max)) ys) - (andmap (λ (z) (z . < . z-min)) zs) (andmap (λ (z) (z . > . z-max)) zs)) - (return empty)) + (when (early-accept? vs x-min x-max y-min y-max z-min z-max) (return (list vs))) + (when (early-reject? vs x-min x-max y-min y-max z-min z-max) (return empty)) (let* ([vss (clip-lines-x-min x-min vs)] [_ (when (empty? vss) (return empty))] [vss (append* (map (λ (vs) (clip-lines-x-max x-max vs)) vss))] @@ -215,3 +188,66 @@ [_ (when (empty? vss) (return empty))] [vss (append* (map (λ (vs) (clip-lines-z-max z-max vs)) vss))]) vss))) + +;; =================================================================================================== +;; Polygon clipping + +(define-syntax-rule (make-clip-polygon/plane a b c d gte?) + (λ (vs ls) + (define v1 (last vs)) + (define x1 (vector-ref v1 0)) + (define y1 (vector-ref v1 1)) + (define z1 (vector-ref v1 2)) + (define v1? (gte? (plane-point-dist a b c d x1 y1 z1) 0)) + + (define-values (new-vs new-ls _x1 _y1 _z1 _v1?) + (for/fold ([vs empty] [ls empty] [x1 x1] [y1 y1] [z1 z1] [v1? v1?]) ([v2 (in-list vs)] + [l (in-list ls)]) + (define x2 (vector-ref v2 0)) + (define y2 (vector-ref v2 1)) + (define z2 (vector-ref v2 2)) + (define v2? (gte? (plane-point-dist a b c d x2 y2 z2) 0)) + (cond [(and v1? v2?) (values (cons v2 vs) (cons l ls) x2 y2 z2 v2?)] + [(not (or v1? v2?)) (values vs ls x2 y2 z2 v2?)] + [else + (define-values (x y z) (plane-line-intersect a b c d x1 y1 z1 x2 y2 z2)) + (if v1? + (values (cons (vector x y z) vs) (cons l ls) x2 y2 z2 v2?) + (values (list* v2 (vector x y z) vs) (list* l #t ls) x2 y2 z2 v2?))]))) + + (values (reverse new-vs) (reverse new-ls)))) + +(define (clip-polygon/plane vs ls plane) + (match-define (vector a b c d) plane) + ((make-clip-polygon/plane a b c d >=) vs ls)) + +;; --------------------------------------------------------------------------------------------------- +;; Clipping inside axial bounding box + +(define-syntax-rule (make-clip-polygon/axis a b c gte?) + (λ (vs ls val) + ((make-clip-polygon/plane a b c (- val) gte?) vs ls))) + +(define clip-polygon-x-min (make-clip-polygon/axis 1 0 0 >=)) +(define clip-polygon-x-max (make-clip-polygon/axis 1 0 0 <=)) +(define clip-polygon-y-min (make-clip-polygon/axis 0 1 0 >=)) +(define clip-polygon-y-max (make-clip-polygon/axis 0 1 0 <=)) +(define clip-polygon-z-min (make-clip-polygon/axis 0 0 1 >=)) +(define clip-polygon-z-max (make-clip-polygon/axis 0 0 1 <=)) + +(define (clip-polygon/bounds vs ls x-min x-max y-min y-max z-min z-max) + (let/ec return + (when (early-accept? vs x-min x-max y-min y-max z-min z-max) (return vs ls)) + (when (early-reject? vs x-min x-max y-min y-max z-min z-max) (return empty empty)) + (let*-values ([(vs ls) (clip-polygon-x-min vs ls x-min)] + [(_) (when (empty? vs) (return empty empty))] + [(vs ls) (clip-polygon-x-max vs ls x-max)] + [(_) (when (empty? vs) (return empty empty))] + [(vs ls) (clip-polygon-y-min vs ls y-min)] + [(_) (when (empty? vs) (return empty empty))] + [(vs ls) (clip-polygon-y-max vs ls y-max)] + [(_) (when (empty? vs) (return empty empty))] + [(vs ls) (clip-polygon-z-min vs ls z-min)] + [(_) (when (empty? vs) (return empty empty))] + [(vs ls) (clip-polygon-z-max vs ls z-max)]) + (values vs ls)))) diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/contour.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/contour.rkt index 9040a1f717..66944701ec 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/contour.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/contour.rkt @@ -1,6 +1,6 @@ #lang racket/base -(require racket/class racket/match racket/list racket/flonum racket/contract +(require racket/class racket/match racket/list racket/contract unstable/latent-contract/defthing plot/utils "../common/utils.rkt") @@ -21,14 +21,11 @@ (when (<= z-min z z-max) (send area put-alpha alpha) (send area put-pen color width style) - (let* ([flonum-ok? (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max)] - [sample (if flonum-ok? (2d-sample-exact->inexact sample) sample)] - [z (if flonum-ok? (exact->inexact z) z)]) - (for-2d-sample - (xa xb ya yb z1 z2 z3 z4) sample - (for ([line (in-list (heights->lines xa xb ya yb z z1 z2 z3 z4))]) - (match-define (list v1 v2) line) - (send area put-line v1 v2))))) + (for-2d-sample + (xa xb ya yb z1 z2 z3 z4) sample + (for ([line (in-list (heights->lines xa xb ya yb z z1 z2 z3 z4))]) + (match-define (list v1 v2) line) + (send area put-line v1 v2)))) (cond [label (line-legend-entry label color width style)] [else empty])) @@ -71,10 +68,7 @@ (let* ([colors (maybe-apply colors zs)] [widths (maybe-apply widths zs)] [styles (maybe-apply styles zs)] - [alphas (maybe-apply alphas zs)] - [flonum-ok? (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max)] - [sample (if flonum-ok? (2d-sample-exact->inexact sample) sample)] - [zs (if flonum-ok? (map exact->inexact zs) zs)]) + [alphas (maybe-apply alphas zs)]) (for ([z (in-list zs)] [color (in-cycle colors)] [width (in-cycle widths)] @@ -141,10 +135,7 @@ [alphas (maybe-apply alphas z-ivls)] [line-colors (maybe-apply line-colors z-ivls)] [line-widths (maybe-apply line-widths z-ivls)] - [line-styles (maybe-apply line-styles z-ivls)] - [flonum-ok? (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max)] - [sample (if flonum-ok? (2d-sample-exact->inexact sample) sample)] - [zs (if flonum-ok? (map exact->inexact zs) zs)]) + [line-styles (maybe-apply line-styles z-ivls)]) (for ([za (in-list zs)] [zb (in-list (rest zs))] [color (in-cycle colors)] diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/isosurface.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/isosurface.rkt index c08993af1e..7739d32c1a 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/isosurface.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/isosurface.rkt @@ -2,7 +2,6 @@ (require racket/class racket/match racket/list racket/flonum racket/contract racket/math unstable/latent-contract/defthing - unstable/flonum plot/utils) (provide (all-defined-out)) @@ -24,13 +23,10 @@ (send area put-alpha alpha) (send area put-brush color style) (send area put-pen line-color line-width line-style) - (let* ([flonum-ok? (flonum-ok-for-4d? x-min x-max y-min y-max z-min z-max d-min d-max)] - [sample (if flonum-ok? (3d-sample-exact->inexact sample) sample)] - [d (if flonum-ok? (exact->inexact d) d)]) - (for-3d-sample - (xa xb ya yb za zb d1 d2 d3 d4 d5 d6 d7 d8) sample - (for ([vs (in-list (heights->cube-polys xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))]) - (send area put-polygon vs)))) + (for-3d-sample + (xa xb ya yb za zb d1 d2 d3 d4 d5 d6 d7 d8) sample + (for ([vs (in-list (heights->cube-polys xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))]) + (send area put-polygon vs))) (cond [label (rectangle-legend-entry label color style line-color line-width line-style)] @@ -85,10 +81,7 @@ [alphas (maybe-apply alphas ds)] [line-colors (maybe-apply line-colors ds)] [line-widths (maybe-apply line-widths ds)] - [line-styles (maybe-apply line-styles ds)] - [flonum-ok? (flonum-ok-for-4d? x-min x-max y-min y-max z-min z-max d-min d-max)] - [sample (if flonum-ok? (3d-sample-exact->inexact sample) sample)] - [ds (if flonum-ok? (map exact->inexact ds) ds)]) + [line-styles (maybe-apply line-styles ds)]) (for ([d (in-list ds)] [color (in-cycle colors)] [style (in-cycle styles)] @@ -153,26 +146,24 @@ (send area put-alpha alpha) (send area put-brush color style) (send area put-pen line-color line-width line-style) - (let* ([flonum-ok? (flonum-ok-for-4d? x-min x-max y-min y-max z-min z-max d-min d-max)] - [sample (if flonum-ok? (3d-sample-exact->inexact sample) sample)] - [d (if flonum-ok? 0.0 0)]) - (for-3d-sample - (xa xb ya yb za zb d1 d2 d3 d4 d5 d6 d7 d8) sample - (cond [(and (xb . > . 0) (ya . < . 0) (yb . > . 0)) - (let* ([yb -0.00001] - [d3 (f xb yb za)] - [d4 (f xa yb za)] - [d7 (f xb yb zb)] - [d8 (f xa yb zb)]) - (draw-cube xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8)) - (let* ([ya 0.00001] - [d1 (f xa ya za)] - [d2 (f xb ya za)] - [d5 (f xa ya zb)] - [d6 (f xb ya zb)]) - (draw-cube xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))] - [else - (draw-cube xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8)]))) + (define d 0) + (for-3d-sample + (xa xb ya yb za zb d1 d2 d3 d4 d5 d6 d7 d8) sample + (cond [(and (xb . > . 0) (ya . < . 0) (yb . > . 0)) + (let* ([yb -0.00001] + [d3 (f xb yb za)] + [d4 (f xa yb za)] + [d7 (f xb yb zb)] + [d8 (f xa yb zb)]) + (draw-cube xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8)) + (let* ([ya 0.00001] + [d1 (f xa ya za)] + [d2 (f xb ya za)] + [d5 (f xa ya zb)] + [d6 (f xb ya zb)]) + (draw-cube xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))] + [else + (draw-cube xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8)])) (cond [label (rectangle-legend-entry label color style line-color line-width line-style)] diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/matrix.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/matrix.rkt deleted file mode 100644 index 5e861fe768..0000000000 --- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/matrix.rkt +++ /dev/null @@ -1,61 +0,0 @@ -#lang racket/base - -;; A small rotation matrix library, used to transform plot coordinates into view coordinates. - -(require racket/match - (only-in math/flonum fl) - racket/flonum) - -(provide m3-apply m3-transpose m3* m3-rotate-z m3-rotate-x m3-scale - fl3-dot) - -(define-syntax-rule (dot x1 y1 z1 x2 y2 z2) (fl+ (fl+ (fl* x1 x2) (fl* y1 y2)) (fl* z1 z2))) - -(define (m3-apply m v) - (match-define (vector v1 v2 v3) m) - (define x (flvector-ref v 0)) - (define y (flvector-ref v 1)) - (define z (flvector-ref v 2)) - (flvector (dot x y z (flvector-ref v1 0) (flvector-ref v1 1) (flvector-ref v1 2)) - (dot x y z (flvector-ref v2 0) (flvector-ref v2 1) (flvector-ref v2 2)) - (dot x y z (flvector-ref v3 0) (flvector-ref v3 1) (flvector-ref v3 2)))) - -(define (m3-transpose m) - (match-define (vector v1 v2 v3) m) - (vector (flvector (flvector-ref v1 0) (flvector-ref v2 0) (flvector-ref v3 0)) - (flvector (flvector-ref v1 1) (flvector-ref v2 1) (flvector-ref v3 1)) - (flvector (flvector-ref v1 2) (flvector-ref v2 2) (flvector-ref v3 2)))) - -(define (m3* m1 m2) - (match-define (vector v1 v2 v3) m1) - (define m (m3-transpose m2)) - (vector (m3-apply m v1) (m3-apply m v2) (m3-apply m v3))) - -(define (m3-rotate-z theta) - (let ([theta (fl theta)]) - (define cos-theta (flcos theta)) - (define sin-theta (flsin theta)) - (vector (flvector cos-theta (- sin-theta) 0.0) - (flvector sin-theta cos-theta 0.0) - (flvector 0.0 0.0 1.0)))) - -(define (m3-rotate-x rho) - (let ([rho (fl rho)]) - (define cos-rho (flcos rho)) - (define sin-rho (flsin rho)) - (vector (flvector 1.0 0.0 0.0) - (flvector 0.0 cos-rho (- sin-rho)) - (flvector 0.0 sin-rho cos-rho)))) - -(define (m3-scale x-scale y-scale z-scale) - (vector (flvector (fl x-scale) 0.0 0.0) - (flvector 0.0 (fl y-scale) 0.0) - (flvector 0.0 0.0 (fl z-scale)))) - -(define (fl3-dot v1 v2) - (fl+ (fl* (flvector-ref v1 0) - (flvector-ref v2 0)) - (fl+ (fl* (flvector-ref v1 1) - (flvector-ref v2 1)) - (fl* (flvector-ref v1 2) - (flvector-ref v2 2))))) diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/plot-area.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/plot-area.rkt index 3caa75f0fe..2799c7c159 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/plot-area.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/plot-area.rkt @@ -1,20 +1,15 @@ #lang racket/base -(require racket/class racket/match racket/list racket/math racket/contract racket/vector - racket/flonum +(require racket/class racket/match racket/list racket/math racket/flonum (only-in math fl flvector->vector vector->flvector) - racket/promise - unstable/flonum "../common/math.rkt" "../common/plot-device.rkt" "../common/ticks.rkt" "../common/draw.rkt" - "../common/contract.rkt" "../common/axis-transform.rkt" "../common/parameters.rkt" - "../common/sample.rkt" "../common/utils.rkt" - "matrix.rkt" + "vector.rkt" "clip.rkt" "bsp-trees.rkt" "bsp.rkt") @@ -104,7 +99,8 @@ (define/public (clear-clip-rect) (set! clipping? #f)) (define (in-bounds? v) - (or (not clipping?) (point-in-bounds? v clip-x-min clip-x-max + (or (not clipping?) (point-in-bounds? v + clip-x-min clip-x-max clip-y-min clip-y-max clip-z-min clip-z-max))) @@ -132,7 +128,7 @@ ;; There are four coordinate systems: ;; 1. Plot coordinates (original, user-facing coordinate system) ;; 2. Normalized coordinates (from plot coordinates: for each axis: transform, center, and scale - ;; to [-1/2,1/2]) - these are always vectors of flonum + ;; to [-0.5,0.5]) - these are always flvectors ;; 3. View coordinates (from normalized coordinates: rotate) ;; 4. Device context coordinates (from view coordinates: project to 2D) @@ -145,40 +141,22 @@ (equal? (plot-y-transform) id-transform) (equal? (plot-z-transform) id-transform))) - (define flonum-ok? (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max)) - (define plot->norm - (if flonum-ok? - (let-map - (x-mid y-mid z-mid x-size y-size z-size) fl - (if identity-transforms? - (match-lambda - [(vector x y z) - (flvector (fl/ (fl- (fl x) x-mid) x-size) - (fl/ (fl- (fl y) y-mid) y-size) - (fl/ (fl- (fl z) z-mid) z-size))]) - (match-lambda - [(vector (? rational? x) (? rational? y) (? rational? z)) - (flvector (fl/ (fl- (fl (fx x)) x-mid) x-size) - (fl/ (fl- (fl (fy y)) y-mid) y-size) - (fl/ (fl- (fl (fz z)) z-mid) z-size))] - [(vector x y z) - (flvector +nan.0 +nan.0 +nan.0)]))) - (if identity-transforms? - (match-lambda - [(vector (? rational? x) (? rational? y) (? rational? z)) - (flvector (fl (/ (- (inexact->exact x) x-mid) x-size)) - (fl (/ (- (inexact->exact y) y-mid) y-size)) - (fl (/ (- (inexact->exact z) z-mid) z-size)))] - [(vector x y z) - (flvector +nan.0 +nan.0 +nan.0)]) - (match-lambda - [(vector (? rational? x) (? rational? y) (? rational? z)) - (flvector (fl (/ (- (inexact->exact (fx x)) x-mid) x-size)) - (fl (/ (- (inexact->exact (fy y)) y-mid) y-size)) - (fl (/ (- (inexact->exact (fz z)) z-mid) z-size)))] - [(vector x y z) - (flvector +nan.0 +nan.0 +nan.0)])))) + (if identity-transforms? + (match-lambda + [(vector (? rational? x) (? rational? y) (? rational? z)) + (flvector (fl (/ (- x x-mid) x-size)) + (fl (/ (- y y-mid) y-size)) + (fl (/ (- z z-mid) z-size)))] + [(vector x y z) + (flvector +nan.0 +nan.0 +nan.0)]) + (match-lambda + [(vector (? rational? x) (? rational? y) (? rational? z)) + (flvector (fl (/ (- (inexact->exact (fx x)) x-mid) x-size)) + (fl (/ (- (inexact->exact (fy y)) y-mid) y-size)) + (fl (/ (- (inexact->exact (fz z)) z-mid) z-size)))] + [(vector x y z) + (flvector +nan.0 +nan.0 +nan.0)]))) (define rotate-theta-matrix (m3-rotate-z theta)) (define rotate-rho-matrix (m3-rotate-x rho)) @@ -789,11 +767,11 @@ ;; Diffuse lighting: typical Lambertian surface model (using absolute value because we ;; can't expect surface normals to point the right direction) (define diff - (cond [diffuse-light? (flabs (fl3-dot normal light-dir))] + (cond [diffuse-light? (flabs (flv3-dot normal light-dir))] [else 1.0])) ;; Specular highlighting: Blinn-Phong model (define spec - (cond [specular-light? (fl* 32.0 (expt (flabs (fl3-dot normal half-dir)) 20.0))] + (cond [specular-light? (fl* 32.0 (expt (flabs (flv3-dot normal half-dir)) 20.0))] [else 0.0])) ;; Blend ambient light with diffuse light, return specular as it is ;; As ambient-light -> 1.0, contribution of diffuse -> 0.0 @@ -809,7 +787,7 @@ vs ls normal) s) (define view-normal (norm->view normal)) - (define cos-view (fl3-dot view-dir view-normal)) + (define cos-view (flv3-dot view-dir view-normal)) (cond [(and (cos-view . < . 0.0) (eq? face 'front)) (void)] [(and (cos-view . > . 0.0) (eq? face 'back)) (void)] @@ -819,7 +797,7 @@ (let ([pen-color (map (λ (v) (+ (* v diff) spec)) pen-color)] [brush-color (map (λ (v) (+ (* v diff) spec)) brush-color)] [vs (map norm->dc vs)]) - ;(send pd set-pen "black" 0.5 'solid) + ;(send pd set-pen "black" 0.5 'solid) ; for BSP debugging (send pd set-pen "black" 0 'transparent) (send pd set-brush brush-color brush-style) (send pd draw-polygon vs) @@ -990,79 +968,79 @@ ;; Drawing shapes (define/public (put-line v1 v2) - (let/ec return - (unless (and (vrational? v1) (vrational? v2)) (return (void))) - (let-values ([(v1 v2) (if clipping? - (clip-line v1 v2 clip-x-min clip-x-max - clip-y-min clip-y-max - clip-z-min clip-z-max) - (values v1 v2))]) - (unless (and v1 v2) (return (void))) - (cond [identity-transforms? - (add-shape! plot3d-area-layer - (line (line-data alpha pen-color pen-width pen-style) - (plot->norm v1) - (plot->norm v2)))] - [else - (define vs (subdivide-line plot->dc v1 v2)) - (add-shape! plot3d-area-layer - (lines (line-data alpha pen-color pen-width pen-style) - (map plot->norm vs)))])))) - - (define/public (put-lines vs) - (for ([vs (vrational-sublists vs)]) - (unless (empty? vs) - (let ([vss (if clipping? - (clip-lines vs clip-x-min clip-x-max - clip-y-min clip-y-max - clip-z-min clip-z-max) - (list vs))]) - (cond [identity-transforms? - (for ([vs (in-list vss)]) + (let ([v1 (exact-vector3d v1)] + [v2 (exact-vector3d v2)]) + (when (and v1 v2) + (let-values ([(v1 v2) (if clipping? + (clip-line/bounds v1 v2 + clip-x-min clip-x-max + clip-y-min clip-y-max + clip-z-min clip-z-max) + (values v1 v2))]) + (when (and v1 v2) + (cond [identity-transforms? + (add-shape! plot3d-area-layer + (line (line-data alpha pen-color pen-width pen-style) + (plot->norm v1) + (plot->norm v2)))] + [else + (define vs (subdivide-line plot->dc v1 v2)) (add-shape! plot3d-area-layer (lines (line-data alpha pen-color pen-width pen-style) - (map plot->norm vs))))] - [else - (for ([vs (in-list vss)]) - (let ([vs (subdivide-lines plot->dc vs)]) - (add-shape! plot3d-area-layer - (lines (line-data alpha pen-color pen-width pen-style) - (map plot->norm vs)))))]))))) + (map plot->norm vs)))])))))) - (define (add-polygon! vs face ls) - (let/ec return - (when (or (empty? vs) (not (andmap vrational? vs))) - (return empty)) - - (define norm-vs (map (λ (v) (flvector->vector (plot->norm v))) vs)) - (define normal (vector->flvector (vnormal norm-vs))) - (define center (vector->flvector (vcenter norm-vs))) - (let*-values - ([(vs ls) (if clipping? - (clip-polygon vs ls - clip-x-min clip-x-max - clip-y-min clip-y-max - clip-z-min clip-z-max) - (values vs ls))] - [(vs ls) (if identity-transforms? - (values vs ls) - (subdivide-polygon plot->dc vs ls))]) - (when (empty? vs) (return empty)) - (add-shape! plot3d-area-layer - (poly (poly-data alpha center - pen-color pen-width pen-style - brush-color brush-style face) - (map plot->norm vs) ls normal))))) + (define/public (put-lines vs) + (for ([vs (in-list (exact-vector3d-sublists vs))]) + (let ([vss (if clipping? + (clip-lines/bounds vs + clip-x-min clip-x-max + clip-y-min clip-y-max + clip-z-min clip-z-max) + (list vs))]) + (cond [identity-transforms? + (for ([vs (in-list vss)]) + (add-shape! plot3d-area-layer + (lines (line-data alpha pen-color pen-width pen-style) + (map plot->norm vs))))] + [else + (for ([vs (in-list vss)]) + (let ([vs (subdivide-lines plot->dc vs)]) + (add-shape! plot3d-area-layer + (lines (line-data alpha pen-color pen-width pen-style) + (map plot->norm vs)))))])))) (define/public (put-polygon vs [face 'both] [ls (make-list (length vs) #t)]) - (add-polygon! vs face ls)) + (let-values ([(vs ls) (exact-polygon3d vs ls)]) + (unless (empty? vs) + (let*-values ([(vs ls) (if clipping? + (clip-polygon/bounds vs ls + clip-x-min clip-x-max + clip-y-min clip-y-max + clip-z-min clip-z-max) + (values vs ls))] + [(vs ls) (if identity-transforms? + (values vs ls) + (subdivide-polygon plot->dc vs ls))]) + (unless (empty? vs) + (define norm-vs (map plot->norm vs)) + (define normal (flv3-normal norm-vs)) + (define center (flv3-center norm-vs)) + (add-shape! plot3d-area-layer + (poly (poly-data alpha center + pen-color pen-width pen-style + brush-color brush-style face) + norm-vs ls normal))))))) (define/public (put-rect r) (when (rect-rational? r) (let ([r (rect-meet r bounds-rect)]) (match-define (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max)) r) - (define v-min (plot->norm (vector x-min y-min z-min))) - (define v-max (plot->norm (vector x-max y-max z-max))) + (define v-min (plot->norm (vector (inexact->exact x-min) + (inexact->exact y-min) + (inexact->exact z-min)))) + (define v-max (plot->norm (vector (inexact->exact x-max) + (inexact->exact y-max) + (inexact->exact z-max)))) (let () (define x-min (flvector-ref v-min 0)) (define y-min (flvector-ref v-min 1)) @@ -1112,13 +1090,14 @@ (define/public (put-text str v [anchor 'center] [angle 0] [dist 0] #:outline? [outline? #f] #:layer [layer plot3d-area-layer]) - (when (and (vrational? v) (in-bounds? v)) - (add-shape! layer (points (text-data alpha anchor angle dist str - font-size font-family text-foreground outline?) - (list (plot->norm v)))))) + (let ([v (exact-vector3d v)]) + (when (and v (in-bounds? v)) + (add-shape! layer (points (text-data alpha anchor angle dist str + font-size font-family text-foreground outline?) + (list (plot->norm v))))))) (define/public (put-glyphs vs symbol size #:layer [layer plot3d-area-layer]) - (let ([vs (filter (λ (v) (and (vrational? v) (in-bounds? v))) vs)]) + (let ([vs (filter (λ (v) (and v (in-bounds? v))) (map exact-vector3d vs))]) (unless (empty? vs) (add-shape! layer (points (glyph-data alpha symbol size pen-color pen-width pen-style @@ -1126,14 +1105,16 @@ (map plot->norm vs)))))) (define/public (put-arrow v1 v2) - (when (and (vrational? v1) (vrational? v2) (in-bounds? v1)) - (cond [(in-bounds? v2) - (define c (v* (v+ v1 v2) 1/2)) - (define outline-color (->brush-color (plot-background))) - (add-shape! plot3d-area-layer - (points (arrow-data alpha (plot->norm v1) (plot->norm v2) - outline-color pen-color pen-width pen-style) - (list (plot->norm c))))] - [else - (put-line v1 v2)]))) + (let ([v1 (exact-vector3d v1)] + [v2 (exact-vector3d v2)]) + (when (and v1 v2 (in-bounds? v1)) + (cond [(in-bounds? v2) + (define c (v* (v+ v1 v2) 1/2)) + (define outline-color (->brush-color (plot-background))) + (add-shape! plot3d-area-layer + (points (arrow-data alpha (plot->norm v1) (plot->norm v2) + outline-color pen-color pen-width pen-style) + (list (plot->norm c))))] + [else + (put-line v1 v2)])))) )) ; end class diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/split.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/split.rkt index fc42973205..a785cc1a85 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/split.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/split.rkt @@ -1,7 +1,6 @@ #lang typed/racket/base -(require racket/match - racket/list +(require racket/list racket/flonum) (provide point3d-plane-dist diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/surface.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/surface.rkt index 7b0bf5469d..ca00ac0190 100644 --- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/surface.rkt +++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/surface.rkt @@ -1,6 +1,6 @@ #lang racket/base -(require racket/class racket/match racket/list racket/flonum racket/contract +(require racket/class racket/match racket/list racket/contract unstable/latent-contract/defthing plot/utils) @@ -22,12 +22,10 @@ (send area put-alpha alpha) (send area put-brush color style) (send area put-pen line-color line-width line-style) - (let* ([flonum-ok? (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max)] - [sample (if flonum-ok? (2d-sample-exact->inexact sample) sample)]) - (for-2d-sample - (xa xb ya yb z1 z2 z3 z4) sample - (define vs (list (vector xa ya z1) (vector xb ya z2) (vector xb yb z3) (vector xa yb z4))) - (send area put-polygon vs))) + (for-2d-sample + (xa xb ya yb z1 z2 z3 z4) sample + (define vs (list (vector xa ya z1) (vector xb ya z2) (vector xb yb z3) (vector xa yb z4))) + (send area put-polygon vs)) (cond [label (rectangle-legend-entry label color style line-color line-width line-style)] [else empty])) diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/vector.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/vector.rkt new file mode 100644 index 0000000000..026a1b9f5b --- /dev/null +++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/vector.rkt @@ -0,0 +1,144 @@ +#lang racket/base + +(require racket/match + racket/list + (only-in math/flonum fl flvector+ flvector-) + (only-in racket/unsafe/ops unsafe-vector-ref) + racket/flonum) + +(provide m3-apply m3-transpose m3* m3-rotate-z m3-rotate-x + flv3-dot + flv3-normal + flv3-center + exact-vector3d + exact-vector3d-sublists + exact-polygon3d) + +(define-syntax-rule (dot x1 y1 z1 x2 y2 z2) (fl+ (fl+ (fl* x1 x2) (fl* y1 y2)) (fl* z1 z2))) + +(define (m3-apply m v) + (match-define (vector v1 v2 v3) m) + (define x (flvector-ref v 0)) + (define y (flvector-ref v 1)) + (define z (flvector-ref v 2)) + (flvector (dot x y z (flvector-ref v1 0) (flvector-ref v1 1) (flvector-ref v1 2)) + (dot x y z (flvector-ref v2 0) (flvector-ref v2 1) (flvector-ref v2 2)) + (dot x y z (flvector-ref v3 0) (flvector-ref v3 1) (flvector-ref v3 2)))) + +(define (m3-transpose m) + (match-define (vector v1 v2 v3) m) + (vector (flvector (flvector-ref v1 0) (flvector-ref v2 0) (flvector-ref v3 0)) + (flvector (flvector-ref v1 1) (flvector-ref v2 1) (flvector-ref v3 1)) + (flvector (flvector-ref v1 2) (flvector-ref v2 2) (flvector-ref v3 2)))) + +(define (m3* m1 m2) + (match-define (vector v1 v2 v3) m1) + (define m (m3-transpose m2)) + (vector (m3-apply m v1) (m3-apply m v2) (m3-apply m v3))) + +(define (m3-rotate-z theta) + (let ([theta (fl theta)]) + (define cos-theta (flcos theta)) + (define sin-theta (flsin theta)) + (vector (flvector cos-theta (- sin-theta) 0.0) + (flvector sin-theta cos-theta 0.0) + (flvector 0.0 0.0 1.0)))) + +(define (m3-rotate-x rho) + (let ([rho (fl rho)]) + (define cos-rho (flcos rho)) + (define sin-rho (flsin rho)) + (vector (flvector 1.0 0.0 0.0) + (flvector 0.0 cos-rho (- sin-rho)) + (flvector 0.0 sin-rho cos-rho)))) + +(define (flv3-dot v1 v2) + (fl+ (fl* (flvector-ref v1 0) + (flvector-ref v2 0)) + (fl+ (fl* (flvector-ref v1 1) + (flvector-ref v2 1)) + (fl* (flvector-ref v1 2) + (flvector-ref v2 2))))) + +(define default-normal (flvector 0.0 -1.0 0.0)) + +(define (flv3-normal vs) + (define n (length vs)) + (cond + [(n . < . 3) default-normal] + [else + (match-define (list v1 v2) (take-right vs 2)) + (define x1 (flvector-ref v1 0)) + (define y1 (flvector-ref v1 1)) + (define z1 (flvector-ref v1 2)) + (define x2 (flvector-ref v2 0)) + (define y2 (flvector-ref v2 1)) + (define z2 (flvector-ref v2 2)) + (define-values (x y z _x1 _y1 _z1 _x2 _y2 _z2) + (for/fold ([x 0.0] [y 0.0] [z 0.0] [x1 x1] [y1 y1] [z1 z1] [x2 x2] [y2 y2] [z2 z2] + ) ([v3 (in-list vs)]) + (define x3 (flvector-ref v3 0)) + (define y3 (flvector-ref v3 1)) + (define z3 (flvector-ref v3 2)) + (define x32 (fl- x3 x2)) + (define y32 (fl- y3 y2)) + (define z32 (fl- z3 z2)) + (define x12 (fl- x1 x2)) + (define y12 (fl- y1 y2)) + (define z12 (fl- z1 z2)) + (values (+ x (fl- (fl* y32 z12) (fl* z32 y12))) + (+ y (fl- (fl* z32 x12) (fl* x32 z12))) + (+ z (fl- (fl* x32 y12) (fl* y32 x12))) + x2 y2 z2 + x3 y3 z3))) + (define m (flsqrt (fl+ (fl* x x) (fl+ (fl* y y) (fl* z z))))) + (if (fl> m 0.0) + (flvector (fl/ x m) (fl/ y m) (fl/ z m)) + default-normal)])) + +(define (flv3-center vs) + (define xs (map (λ (v) (flvector-ref v 0)) vs)) + (define ys (map (λ (v) (flvector-ref v 1)) vs)) + (define zs (map (λ (v) (flvector-ref v 2)) vs)) + (flvector (* 0.5 (+ (apply min xs) (apply max xs))) + (* 0.5 (+ (apply min ys) (apply max ys))) + (* 0.5 (+ (apply min zs) (apply max zs))))) + +(define (exact-vector3d v) + (define n (vector-length v)) + (cond [(= n 3) + (define v1 (unsafe-vector-ref v 0)) + (define v2 (unsafe-vector-ref v 1)) + (define v3 (unsafe-vector-ref v 2)) + (cond [(and (exact? v1) (exact? v2) (exact? v3)) v] + [(and (rational? v1) (rational? v2) (rational? v3)) + (vector (inexact->exact v1) (inexact->exact v2) (inexact->exact v3))] + [else #f])] + [else #f])) + +(define (sublists vs) + (define vss + (for/fold ([vss (list null)]) ([v (in-list vs)]) + (cond [v (cons (cons v (car vss)) (cdr vss))] + [(null? (car vss)) vss] + [else (cons null vss)]))) + (cond [(null? (car vss)) (cdr vss)] + [else vss])) + +(define (exact-vector3d-sublists vs) + (sublists (map exact-vector3d vs))) + +(define (exact-polygon3d vs ls) + (cond + [(null? vs) (values null null)] + [else + (define-values (new-vs new-ls _) + (for/fold ([vs null] [ls null] [v1 (last vs)]) ([v2 (in-list vs)] + [l (in-list ls)]) + (cond [(equal? v1 v2) (values vs ls v2)] + [else + (define exact-v2 (exact-vector3d v2)) + (if exact-v2 + (values (cons exact-v2 vs) (cons l ls) v2) + (values vs ls v2))]))) + (values (reverse new-vs) (reverse new-ls))])) diff --git a/pkgs/plot-pkgs/plot-test/plot/tests/extreme-bounds-tests.rkt b/pkgs/plot-pkgs/plot-test/plot/tests/extreme-bounds-tests.rkt index 0d8bfe2e5d..2f2c1432ef 100644 --- a/pkgs/plot-pkgs/plot-test/plot/tests/extreme-bounds-tests.rkt +++ b/pkgs/plot-pkgs/plot-test/plot/tests/extreme-bounds-tests.rkt @@ -6,18 +6,20 @@ (module config info (define timeout 150))) +(printf "Point at 0,0:~n") (plot (points '(#(0 0))) #:x-min +min.0 #:x-max (flstep +min.0 1000) #:y-min 0 #:y-max 1) +(printf "Point at 0,0,0:~n") (plot3d (points3d '(#(0 0 0))) #:x-min +min.0 #:x-max (flstep +min.0 1000) #:y-min 0 #:y-max 1 #:z-min 0 #:z-max 1) -;; this should go all the way from the bottom to the top +(printf "Steps should appear all the way from the bottom to the top:~n") (plot (function exact->inexact (flstep 0.0 -1) (flstep 0.0 2))) -;; this should go all the way from the bottom to the top +(printf "Steps should appear all the way from the bottom to the top:~n") (plot3d (surface3d (λ (x y) (exact->inexact x)) (flstep 0.0 -1) (flstep 0.0 2) (flstep 0.0 -1) (flstep 0.0 2))) diff --git a/pkgs/plot-pkgs/plot-test/plot/tests/plot2d-tests.rkt b/pkgs/plot-pkgs/plot-test/plot/tests/plot2d-tests.rkt index fa831585f5..5c535bb5c6 100644 --- a/pkgs/plot-pkgs/plot-test/plot/tests/plot2d-tests.rkt +++ b/pkgs/plot-pkgs/plot-test/plot/tests/plot2d-tests.rkt @@ -89,6 +89,7 @@ (plot (list (points (map vector ys xs) #:sym 'full6star #:color 1) (points (map vector xs ys) #:sym "B" #:color 3)))) +(printf "This plot should be blank:~n") (time (plot (vector-field (λ (x y) (vector +nan.0 +nan.0))) #:x-min -2 #:x-max 2 #:y-min -2 #:y-max 2)) diff --git a/pkgs/plot-pkgs/plot-test/plot/tests/plot3d-bsp-tests.rkt b/pkgs/plot-pkgs/plot-test/plot/tests/plot3d-bsp-tests.rkt index be8ed2a60c..4fbe061bc9 100644 --- a/pkgs/plot-pkgs/plot-test/plot/tests/plot3d-bsp-tests.rkt +++ b/pkgs/plot-pkgs/plot-test/plot/tests/plot3d-bsp-tests.rkt @@ -12,18 +12,10 @@ (plot-y-transform cbrt-transform) (plot-z-transform cbrt-transform)) -#; -(plot3d (contour-intervals3d * -1 1 -1 1 #:samples 3 #:line-widths '(2))) - -(plot3d (list (surface3d * -1 1 -1 1 #:samples 2) +(plot3d (list (contour-intervals3d * -1 1 -1 1 #:samples 2) (isosurface3d (λ (x y z) z) 0 -1 1 -1 1 -1 1 - #:samples 2 #:color 3) - #; - (isosurface3d (λ (x y z) (+ x y -0.5)) - 0 -1 1 -1 1 -1 1 - #:samples 2 #:color 2) - )) + #:samples 2 #:color 3))) (plot3d (list (parametric3d (λ (t) (list (* 0.5 (cos (* 2 pi (- t)))) @@ -42,8 +34,7 @@ #:samples 53 #:color "red" #:width 15 - #:alpha 0.6) - )) + #:alpha 0.6))) (plot3d (list (isosurface3d (λ (x y z) (+ x y z)) 0 -1 1 -1 1 -1 1 #:samples 2 @@ -114,6 +105,10 @@ ;#:out-file "test.pdf" ))) +(define (f2 x y) + (let ([x (fl x)] [y (fl y)]) + (- (sqrt (+ (abs y) (abs x)))))) + (plot3d (list (surface3d * -1 1 -1 1 #:samples 6 #:alpha 0.75 #:color 1) (surface3d (λ (x y) (+ 0.1 (* x y))) -1 1 -1 1 #:samples 6 #:alpha 0.75 #:color 2) (surface3d (λ (x y) (+ 0.2 (* x y))) -1 1 -1 1 #:samples 6 #:alpha 0.75 #:color 3 diff --git a/pkgs/plot-pkgs/plot-test/plot/tests/typed/parameter-regression-tests.rkt b/pkgs/plot-pkgs/plot-test/plot/tests/typed/parameter-regression-tests.rkt index d355e850b2..fe864caebd 100644 --- a/pkgs/plot-pkgs/plot-test/plot/tests/typed/parameter-regression-tests.rkt +++ b/pkgs/plot-pkgs/plot-test/plot/tests/typed/parameter-regression-tests.rkt @@ -14,7 +14,8 @@ (list (ivl 1 2) (ivl 1 2))))) ;; --------------------------------------------------------------------------------------------------- -;; The following two plots should be the same + +(printf "The following two plots should be the same:~n") (contour-colors (λ (n) '((255 0 0) (0 0 255)))) (contour-widths (λ (n) '(.25 .5 .75))) @@ -33,7 +34,8 @@ (plot (contour-intervals + 0 1 0 1)) ;; --------------------------------------------------------------------------------------------------- -;; The following two plots should be the same + +(printf "The following two plots should be the same:~n") (stacked-histogram-alphas (λ (n) '(0.25 0.5 0.75))) (stacked-histogram-colors (λ (n) '(1 2 3))) @@ -50,7 +52,8 @@ (list 'c '()) (list 'd '(1/2))))) ;; --------------------------------------------------------------------------------------------------- -;; The following two plots should be the same + +(printf "The following two plots should be the same:~n") (contour-interval-line-colors (λ (n) '(1 2 3))) (contour-interval-line-widths (λ (n) '(0.25 0.5 0.75))) @@ -63,7 +66,8 @@ (plot3d (contour-intervals3d + 0 1 0 1)) ;; --------------------------------------------------------------------------------------------------- -;; The following two plots should be the same + +(printf "The following two plots should be the same:~n") (isosurface-colors (λ (n) '(1 2 3))) (isosurface-styles (λ (n) '(0 1 2)))