Exposed marching squares and cubes algorithms through plot/utils, speed improvements, fixes
This commit is contained in:
Neil Toronto
parent
ee71e01c13
commit
56f70fb4f2
|
|
@ -1,22 +1,12 @@
|
|||
#lang racket/base
|
||||
|
||||
(require racket/flonum racket/fixnum racket/list racket/match
|
||||
(require racket/contract racket/flonum racket/fixnum racket/list racket/match racket/unsafe/ops
|
||||
(for-syntax racket/base racket/syntax racket/match racket/list)
|
||||
"math.rkt")
|
||||
"math.rkt"
|
||||
"marching-utils.rkt"
|
||||
"contract-doc.rkt")
|
||||
|
||||
(provide scale-normalized-poly
|
||||
heights->cube-polys)
|
||||
|
||||
(define (scale-normalized-poly poly xa xb ya yb za zb)
|
||||
(for/list ([uvw (in-list poly)])
|
||||
(match-define (vector u v w) uvw)
|
||||
(vector (blend xb xa u) (blend yb ya v) (blend zb za w))))
|
||||
|
||||
(define-syntax-rule (solve-t d da db)
|
||||
(fl/ (fl- d da) (fl- db da)))
|
||||
|
||||
(define (flavg4 d1 d2 d3 d4)
|
||||
(fl* 0.25 (fl+ (fl+ (fl+ d1 d2) d3) d4)))
|
||||
(provide heights->cube-polys heights->cube-polys:doc)
|
||||
|
||||
(define-for-syntax (->datum x)
|
||||
(if (syntax? x) (syntax->datum x) x))
|
||||
|
|
@ -81,7 +71,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 (flavg4 d1 d2 d3 d4))
|
||||
(define da (unsafe-flavg4 d1 d2 d3 d4))
|
||||
(cond
|
||||
[(test? da d)
|
||||
(list (list (edge-1-2 d d1 d2) (edge-2-3 d d2 d3)
|
||||
|
|
@ -92,8 +82,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 fl>=))
|
||||
(define known-cube-0101-1111 (make-known-cube-1010-0000 fl<))
|
||||
(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-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))
|
||||
|
|
@ -102,7 +92,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 (flavg4 d2 d6 d7 d3))
|
||||
(define da (unsafe-flavg4 d2 d6 d7 d3))
|
||||
(cond
|
||||
[(test? da d)
|
||||
(list (list (edge-1-5 d d1 d5) (edge-2-6 d d2 d6)
|
||||
|
|
@ -115,12 +105,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 fl>=))
|
||||
(define known-cube-0011-1101 (make-known-cube-1100-0010 fl<))
|
||||
(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 ((make-known-cube-1100-0011 test?) d d1 d2 d3 d4 d5 d6 d7 d8)
|
||||
(define da (flavg4 d1 d5 d8 d4))
|
||||
(define db (flavg4 d2 d6 d7 d3))
|
||||
(define da (unsafe-flavg4 d1 d5 d8 d4))
|
||||
(define db (unsafe-flavg4 d2 d6 d7 d3))
|
||||
(cond
|
||||
[(and (test? da d) (test? db d))
|
||||
(list (list (edge-1-5 d d1 d5) (edge-5-8 d d5 d8)
|
||||
|
|
@ -161,8 +151,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 fl>=))
|
||||
(define known-cube-0011-1100 (make-known-cube-1100-0011 fl<))
|
||||
(define known-cube-1100-0011 (make-known-cube-1100-0011 unsafe-fl>=))
|
||||
(define known-cube-0011-1100 (make-known-cube-1100-0011 unsafe-fl<))
|
||||
|
||||
#|
|
||||
Cube vertex numbers:
|
||||
|
|
@ -178,12 +168,12 @@ Cube vertex numbers:
|
|||
|
||||
|
||||
(define ((make-known-cube-1010-0101 test?) d d1 d2 d3 d4 d5 d6 d7 d8)
|
||||
(define da (flavg4 d1 d2 d3 d4))
|
||||
(define db (flavg4 d1 d5 d8 d4))
|
||||
(define dc (flavg4 d3 d4 d8 d7))
|
||||
(define dd (flavg4 d1 d2 d6 d5))
|
||||
(define de (flavg4 d2 d3 d7 d6))
|
||||
(define df (flavg4 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))
|
||||
(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))
|
||||
|
|
@ -214,25 +204,25 @@ 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 fl>=))
|
||||
(define known-cube-0101-1010 (make-known-cube-1010-0101 fl<))
|
||||
(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 ((make-known-cube-1110-0001 fl>=) d d1 d2 d3 d4 d5 d6 d7 d8)
|
||||
(define da (flavg4 d1 d5 d8 d4))
|
||||
(define db (flavg4 d7 d8 d4 d3))
|
||||
(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))
|
||||
(cond
|
||||
[(and (da . fl>= . d) (db . fl>= . d))
|
||||
[(and (da . unsafe-fl>= . d) (db . unsafe-fl>= . 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 . fl>= . d)
|
||||
[(da . unsafe-fl>= . d)
|
||||
(define ec (v* (v+ (edge-1-5 d d1 d5) (edge-3-7 d d3 d7)) 0.5))
|
||||
(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 . fl>= . d)
|
||||
[(db . unsafe-fl>= . d)
|
||||
(define ec (v* (v+ (edge-1-5 d d1 d5) (edge-3-7 d d3 d7)) 0.5))
|
||||
(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))
|
||||
|
|
@ -244,8 +234,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 fl>=))
|
||||
(define known-cube-0001-1110 (make-known-cube-1110-0001 fl<))
|
||||
(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-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)
|
||||
|
|
@ -260,9 +250,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 (flavg4 d1 d2 d3 d4))
|
||||
(define db (flavg4 d1 d5 d8 d4))
|
||||
(define dc (flavg4 d3 d4 d8 d7))
|
||||
(define da (unsafe-flavg4 d1 d2 d3 d4))
|
||||
(define db (unsafe-flavg4 d1 d5 d8 d4))
|
||||
(define dc (unsafe-flavg4 d3 d4 d8 d7))
|
||||
(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))
|
||||
|
|
@ -280,8 +270,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 fl>=))
|
||||
(define known-cube-0101-1110 (make-known-cube-1010-0001 fl<))
|
||||
(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-for-syntax known-cubes
|
||||
'((0 0 0 0 0 0 0 0)
|
||||
|
|
@ -317,7 +307,7 @@ Cube vertex numbers:
|
|||
|
||||
(define (mirror-vec-d v)
|
||||
(match-define (vector x y d) v)
|
||||
(vector x y (fl- 1.0 d)))
|
||||
(vector x y (unsafe-fl- 1.0 d)))
|
||||
|
||||
(define ((mirror-cube-d f) d d1 d2 d3 d4 d5 d6 d7 d8)
|
||||
(map (λ (poly) (map mirror-vec-d poly))
|
||||
|
|
@ -325,7 +315,7 @@ Cube vertex numbers:
|
|||
|
||||
(define (mirror-vec-y v)
|
||||
(match-define (vector x y d) v)
|
||||
(vector x (fl- 1.0 y) d))
|
||||
(vector x (unsafe-fl- 1.0 y) d))
|
||||
|
||||
(define ((mirror-cube-y f) d d1 d2 d3 d4 d5 d6 d7 d8)
|
||||
(map (λ (poly) (map mirror-vec-y poly))
|
||||
|
|
@ -333,7 +323,7 @@ Cube vertex numbers:
|
|||
|
||||
(define (mirror-vec-x v)
|
||||
(match-define (vector x y d) v)
|
||||
(vector (fl- 1.0 x) y d))
|
||||
(vector (unsafe-fl- 1.0 x) y d))
|
||||
|
||||
(define ((mirror-cube-x f) d d1 d2 d3 d4 d5 d6 d7 d8)
|
||||
(map (λ (poly) (map mirror-vec-x poly))
|
||||
|
|
@ -343,7 +333,7 @@ Cube vertex numbers:
|
|||
|
||||
(define (unrotate-vec-d v)
|
||||
(match-define (vector x y d) v)
|
||||
(vector (fl- 1.0 y) x d))
|
||||
(vector (unsafe-fl- 1.0 y) x d))
|
||||
|
||||
(define ((rotate-cube-d f) d d1 d2 d3 d4 d5 d6 d7 d8)
|
||||
(map (λ (poly) (map unrotate-vec-d poly))
|
||||
|
|
@ -351,7 +341,7 @@ Cube vertex numbers:
|
|||
|
||||
(define (unrotate-vec-y v)
|
||||
(match-define (vector x y d) v)
|
||||
(vector d y (fl- 1.0 x)))
|
||||
(vector d y (unsafe-fl- 1.0 x)))
|
||||
|
||||
(define ((rotate-cube-y f) d d1 d2 d3 d4 d5 d6 d7 d8)
|
||||
(map (λ (poly) (map unrotate-vec-y poly))
|
||||
|
|
@ -359,7 +349,7 @@ Cube vertex numbers:
|
|||
|
||||
(define (unrotate-vec-x v)
|
||||
(match-define (vector x y d) v)
|
||||
(vector x (fl- 1.0 d) y))
|
||||
(vector x (unsafe-fl- 1.0 d) y))
|
||||
|
||||
(define ((rotate-cube-x f) d d1 d2 d3 d4 d5 d6 d7 d8)
|
||||
(map (λ (poly) (map unrotate-vec-x poly))
|
||||
|
|
@ -369,7 +359,7 @@ Cube vertex numbers:
|
|||
|
||||
(define (rotate-vec-d v)
|
||||
(match-define (vector x y d) v)
|
||||
(vector y (fl- 1.0 x) d))
|
||||
(vector y (unsafe-fl- 1.0 x) d))
|
||||
|
||||
(define ((unrotate-cube-d f) d d1 d2 d3 d4 d5 d6 d7 d8)
|
||||
(map (λ (poly) (map rotate-vec-d poly))
|
||||
|
|
@ -377,7 +367,7 @@ Cube vertex numbers:
|
|||
|
||||
(define (rotate-vec-y v)
|
||||
(match-define (vector x y d) v)
|
||||
(vector (fl- 1.0 d) y x))
|
||||
(vector (unsafe-fl- 1.0 d) y x))
|
||||
|
||||
(define ((unrotate-cube-y f) d d1 d2 d3 d4 d5 d6 d7 d8)
|
||||
(map (λ (poly) (map rotate-vec-y poly))
|
||||
|
|
@ -385,7 +375,7 @@ Cube vertex numbers:
|
|||
|
||||
(define (rotate-vec-x v)
|
||||
(match-define (vector x y d) v)
|
||||
(vector x d (fl- 1.0 y)))
|
||||
(vector x d (unsafe-fl- 1.0 y)))
|
||||
|
||||
(define ((unrotate-cube-x f) d d1 d2 d3 d4 d5 d6 d7 d8)
|
||||
(map (λ (poly) (map rotate-vec-x poly))
|
||||
|
|
@ -518,44 +508,34 @@ Cube vertex numbers:
|
|||
(make-cube-dispatch-table))
|
||||
|
||||
(define-syntax-rule (add-digit j idx)
|
||||
(fx+ (fx* idx 2) j))
|
||||
(unsafe-fx+ (unsafe-fx* idx 2) j))
|
||||
|
||||
(define (heights->cube-polys d d1 d2 d3 d4 d5 d6 d7 d8)
|
||||
(define j1 (if (d1 . fl< . d) 0 1))
|
||||
(define j2 (if (d2 . fl< . d) 0 1))
|
||||
(define j3 (if (d3 . fl< . d) 0 1))
|
||||
(define j4 (if (d4 . fl< . d) 0 1))
|
||||
(define j5 (if (d5 . fl< . d) 0 1))
|
||||
(define j6 (if (d6 . fl< . d) 0 1))
|
||||
(define j7 (if (d7 . fl< . d) 0 1))
|
||||
(define j8 (if (d8 . fl< . d) 0 1))
|
||||
(define facet-num
|
||||
(add-digit
|
||||
j1 (add-digit
|
||||
j2 (add-digit
|
||||
j3 (add-digit
|
||||
j4 (add-digit
|
||||
j5 (add-digit
|
||||
j6 (add-digit j7 j8))))))))
|
||||
(define f (vector-ref cube-dispatch-table facet-num))
|
||||
(f d d1 d2 d3 d4 d5 d6 d7 d8))
|
||||
|
||||
#|
|
||||
(require "../plot3d.rkt")
|
||||
|
||||
(define ((shift-x dx) v)
|
||||
(match-define (vector x y d) v)
|
||||
(vector (+ x dx) y d))
|
||||
|
||||
(define (polys->lines lsts [dx 0.0] #:color [color 1] #:alpha [alpha 1])
|
||||
(map (λ (lst) (lines3d (map (shift-x dx)
|
||||
(append lst (list (first lst))))
|
||||
#:color color #:alpha alpha))
|
||||
lsts))
|
||||
|
||||
(plot3d (polys->lines (heights->cube-polys 0.5
|
||||
1.0 0.0 1.0 0.0
|
||||
0.0 1.0 0.0 1.0)
|
||||
#:alpha 1/2)
|
||||
#:x-min 0 #:x-max 1 #:y-min 0 #:y-max 1 #:d-min 0 #:d-max 1)
|
||||
|#
|
||||
(defproc (heights->cube-polys [xa real?] [xb real?] [ya real?] [yb real?] [za real?] [zb real?]
|
||||
[d real?]
|
||||
[d1 real?] [d2 real?] [d3 real?] [d4 real?]
|
||||
[d5 real?] [d6 real?] [d7 real?] [d8 real?]
|
||||
) (listof (vector/c real? real? real?))
|
||||
(check-all-real! 'heights->cube-polys xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8)
|
||||
(let-exact->inexact
|
||||
(xa xb ya yb za zb 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 facet-num
|
||||
(add-digit
|
||||
j1 (add-digit
|
||||
j2 (add-digit
|
||||
j3 (add-digit
|
||||
j4 (add-digit
|
||||
j5 (add-digit
|
||||
j6 (add-digit j7 j8))))))))
|
||||
(define f (vector-ref cube-dispatch-table facet-num))
|
||||
(for/list ([poly (in-list (f d d1 d2 d3 d4 d5 d6 d7 d8))])
|
||||
(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))))))
|
||||
|
|
|
|||
|
|
@ -3,22 +3,12 @@
|
|||
(require racket/flonum racket/fixnum racket/list racket/match racket/unsafe/ops racket/contract
|
||||
(for-syntax racket/base)
|
||||
"math.rkt"
|
||||
"marching-utils.rkt"
|
||||
"contract-doc.rkt")
|
||||
|
||||
(provide heights->lines heights->polys
|
||||
heights->lines:doc heights->polys:doc)
|
||||
|
||||
;; Returns the interpolated distance of z from za toward zb
|
||||
;; Examples: if z = za, this returns 0.0
|
||||
;; 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)
|
||||
(unsafe-fl/ (unsafe-fl- z za) (unsafe-fl- zb za)))
|
||||
|
||||
(define-syntax-rule (solve-z za zb t)
|
||||
(unsafe-fl+ (unsafe-fl* t zb) (unsafe-fl* (unsafe-fl- 1.0 t) za)))
|
||||
|
||||
#|
|
||||
Z values are at these normalized coordinates:
|
||||
|
||||
|
|
@ -50,9 +40,6 @@ above.
|
|||
(match-let ([(vector x y z) v])
|
||||
(vector x (fl- 1.0 y) z)))
|
||||
|
||||
(define-syntax-rule (flavg4 z1 z2 z3 z4)
|
||||
(unsafe-fl* 0.25 (unsafe-fl+ (unsafe-fl+ (unsafe-fl+ z1 z2) z3) z4)))
|
||||
|
||||
;; =============================================================================
|
||||
;; Contour lines
|
||||
|
||||
|
|
@ -116,7 +103,7 @@ 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 (flavg4 z1 z2 z3 z4)])
|
||||
(let ([z5 (unsafe-flavg4 z1 z2 z3 z4)])
|
||||
(if (test? z5 z)
|
||||
(list (vector (solve-t z z1 z2) 0.0
|
||||
1.0 (solve-t z z2 z3))
|
||||
|
|
@ -130,19 +117,6 @@ above.
|
|||
(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 (check-all-real! stx)
|
||||
(syntax-case stx ()
|
||||
[(_ name id ...)
|
||||
(with-syntax ([(i ...) (build-list (length (syntax->list #'(id ...))) values)])
|
||||
(define ids (syntax->list #'(id ...)))
|
||||
#`(begin (unless (real? id)
|
||||
(raise-type-error name "real number" i #,@ids))
|
||||
...))]))
|
||||
|
||||
(define-syntax-rule (let-exact->inexact (id ...) body0 body ...)
|
||||
(let ([id (exact->inexact id)] ...)
|
||||
body0 body ...))
|
||||
|
||||
(defproc (heights->lines [xa real?] [xb real?] [ya real?] [yb real?]
|
||||
[z real?] [z1 real?] [z2 real?] [z3 real?] [z4 real?]
|
||||
) (list/c (vector/c real? real? real?) (vector/c real? real? real?))
|
||||
|
|
@ -187,8 +161,8 @@ above.
|
|||
(lines0000 z z1 z2 z3 z4))))))
|
||||
(for/list ([line (in-list lines)])
|
||||
(match-define (vector u1 v1 u2 v2) line)
|
||||
(list (vector (solve-z xa xb u1) (solve-z ya yb v1) z)
|
||||
(vector (solve-z xa xb u2) (solve-z ya yb v2) z)))))
|
||||
(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
|
||||
|
|
@ -417,7 +391,7 @@ above.
|
|||
(vector 0.0 (solve-t za z1 z4) za))))
|
||||
|
||||
(define (polys1010 za zb z1 z2 z3 z4)
|
||||
(define z5 (flavg4 z1 z2 z3 z4))
|
||||
(define z5 (unsafe-flavg4 z1 z2 z3 z4))
|
||||
(cond [(z5 . unsafe-fl< . za) (polys10100 za zb z1 z2 z3 z4)]
|
||||
; (z5 . >= . zb) is impossible
|
||||
[else (polys10101 za zb z1 z2 z3 z4)]))
|
||||
|
|
@ -441,7 +415,7 @@ above.
|
|||
(vector 0.0 (solve-t zb z1 z4) zb))))
|
||||
|
||||
(define (polys1212 za zb z1 z2 z3 z4)
|
||||
(define z5 (flavg4 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)]
|
||||
; (z5 . < . za) is impossible
|
||||
[else (polys1212-1 za zb z1 z2 z3 z4)]))
|
||||
|
|
@ -470,7 +444,7 @@ above.
|
|||
(vector (solve-t zb z4 z3) 1.0 zb))))
|
||||
|
||||
(define (polys0212 za zb z1 z2 z3 z4)
|
||||
(define z5 (flavg4 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)]
|
||||
; handling (z5 . < . za) separately results in a non-convex polygon
|
||||
[else (polys0212-2 za zb z1 z2 z3 z4)]))
|
||||
|
|
@ -498,7 +472,7 @@ above.
|
|||
(vector (solve-t za z4 z3) 1.0 za))))
|
||||
|
||||
(define (polys2010 za zb z1 z2 z3 z4)
|
||||
(define z5 (flavg4 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)]
|
||||
; handling (z5 . >= . zb) separately results in a non-convex polygon
|
||||
[else (polys2010-0 za zb z1 z2 z3 z4)]))
|
||||
|
|
@ -541,7 +515,7 @@ above.
|
|||
(vector (solve-t zb z4 z3) 1.0 zb))))
|
||||
|
||||
(define (polys0202 za zb z1 z2 z3 z4)
|
||||
(define z5 (flavg4 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)]
|
||||
[else (polys0202-2 za zb z1 z2 z3 z4)]))
|
||||
|
|
@ -593,7 +567,7 @@ above.
|
|||
[za real?] [zb real?]
|
||||
[z1 real?] [z2 real?] [z3 real?] [z4 real?]
|
||||
) (listof (vector/c real? real? real?))
|
||||
(check-all-real! xa xb ya yb za zb z1 z2 z3 z4)
|
||||
(check-all-real! 'heights->polys xa xb ya yb za zb z1 z2 z3 z4)
|
||||
(let-exact->inexact
|
||||
(xa xb ya yb za zb z1 z2 z3 z4)
|
||||
(define t1 (if (z1 . unsafe-fl< . za) 0 (if (z1 . unsafe-fl< . zb) 1 2)))
|
||||
|
|
@ -611,4 +585,4 @@ above.
|
|||
(vector xb yb z3) (vector xa yb z4))]
|
||||
[else (for/list ([uv (in-list poly)])
|
||||
(match-define (vector u v z) uv)
|
||||
(vector (solve-z xa xb u) (solve-z ya yb v) z))]))))
|
||||
(vector (unsolve-t xa xb u) (unsolve-t ya yb v) z))]))))
|
||||
|
|
|
|||
32
collects/plot/common/marching-utils.rkt
Normal file
32
collects/plot/common/marching-utils.rkt
Normal file
|
|
@ -0,0 +1,32 @@
|
|||
#lang racket/base
|
||||
|
||||
(require (for-syntax racket/base) racket/unsafe/ops)
|
||||
|
||||
(provide (all-defined-out))
|
||||
|
||||
;; Returns the interpolated distance of z from za toward zb
|
||||
;; Examples: if z = za, this returns 0.0
|
||||
;; 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)
|
||||
(unsafe-fl/ (unsafe-fl- z za) (unsafe-fl- zb za)))
|
||||
|
||||
(define-syntax-rule (unsolve-t za zb t)
|
||||
(unsafe-fl+ (unsafe-fl* t zb) (unsafe-fl* (unsafe-fl- 1.0 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)))
|
||||
|
||||
(define-syntax (check-all-real! stx)
|
||||
(syntax-case stx ()
|
||||
[(_ name id ...)
|
||||
(with-syntax ([(i ...) (build-list (length (syntax->list #'(id ...))) values)])
|
||||
(define ids (syntax->list #'(id ...)))
|
||||
#`(begin (unless (real? id)
|
||||
(raise-type-error name "real number" i #,@ids))
|
||||
...))]))
|
||||
|
||||
(define-syntax-rule (let-exact->inexact (id ...) body0 body ...)
|
||||
(let ([id (exact->inexact id)] ...)
|
||||
body0 body ...))
|
||||
|
|
@ -21,7 +21,7 @@
|
|||
(cond [(not (flonum? x)) (raise-type-error 'flblend "flonum" 0 x y α)]
|
||||
[(not (flonum? y)) (raise-type-error 'flblend "flonum" 1 x y α)]
|
||||
[(not (flonum? α)) (raise-type-error 'flblend "flonum" 2 x y α)]
|
||||
[else (unsafe-fl+ (unsafe-fl* α x) (unsafe-fl* (unsafe-fl- 1 α) y))]))
|
||||
[else (unsafe-fl+ (unsafe-fl* α x) (unsafe-fl* (unsafe-fl- 1.0 α) y))]))
|
||||
|
||||
(defproc (flatan2 [y flonum?] [x flonum?]) flonum?
|
||||
(cond [(not (flonum? y)) (raise-type-error 'flatan2 "flonum" 0 x y)]
|
||||
|
|
@ -35,13 +35,13 @@
|
|||
(unsafe-fl+ sum y))]))
|
||||
|
||||
(defproc (flmodulo [x flonum?] [y flonum?]) flonum?
|
||||
(cond [(not (flonum? x)) (raise-type-error 'real-modulo "flonum" 0 x y)]
|
||||
[(not (flonum? y)) (raise-type-error 'real-modulo "flonum" 1 x y)]
|
||||
(cond [(not (flonum? x)) (raise-type-error 'flmodulo "flonum" 0 x y)]
|
||||
[(not (flonum? y)) (raise-type-error 'flmodulo "flonum" 1 x y)]
|
||||
[else (unsafe-fl- x (unsafe-fl* y (unsafe-flfloor (unsafe-fl/ x y))))]))
|
||||
|
||||
(define fldistance
|
||||
(case-lambda
|
||||
[() 0]
|
||||
[() 0.0]
|
||||
[(x) (if (flonum? x) (abs x) (raise-type-error 'fldistance "flonum" x))]
|
||||
[(x y) (cond [(not (flonum? x)) (raise-type-error 'fldistance "flonum" 0 x y)]
|
||||
[(not (flonum? y)) (raise-type-error 'fldistance "flonum" 1 x y)]
|
||||
|
|
@ -406,18 +406,24 @@
|
|||
(defproc (vcenter [vs (listof (vectorof real?))]) (vectorof real?)
|
||||
(match vs
|
||||
[(list (vector xs ys) ...)
|
||||
(define mins (vector (apply min* xs) (apply min* ys)))
|
||||
(define maxs (vector (apply max* xs) (apply max* ys)))
|
||||
(unrolled-vmap2 'center-coord (λ (x1 x2) (* 1/2 (+ x1 x2))) mins maxs)]
|
||||
(define x-min (apply min* xs))
|
||||
(define x-max (apply max* xs))
|
||||
(define y-min (apply min* ys))
|
||||
(define y-max (apply max* ys))
|
||||
(vector (* 1/2 (+ x-min x-max)) (* 1/2 (+ y-min y-max)))]
|
||||
[(list (vector xs ys zs) ...)
|
||||
(define mins (vector (apply min* xs) (apply min* ys) (apply min* zs)))
|
||||
(define maxs (vector (apply max* xs) (apply max* ys) (apply max* zs)))
|
||||
(unrolled-vmap2 'center-coord (λ (x1 x2) (* 1/2 (+ x1 x2))) mins maxs)]
|
||||
(define x-min (apply min* xs))
|
||||
(define x-max (apply max* xs))
|
||||
(define y-min (apply min* ys))
|
||||
(define y-max (apply max* ys))
|
||||
(define z-min (apply min* zs))
|
||||
(define z-max (apply max* zs))
|
||||
(vector (* 1/2 (+ x-min x-max)) (* 1/2 (+ y-min y-max)) (* 1/2 (+ z-min z-max)))]
|
||||
[_
|
||||
(define xss (apply vector-map list vs))
|
||||
(define mins (vector-map (λ (xs) (apply min xs)) xss))
|
||||
(define maxs (vector-map (λ (xs) (apply max xs)) xss))
|
||||
(unrolled-vmap2 'center-coord (λ (x1 x2) (* 1/2 (+ x1 x2))) mins maxs)]))
|
||||
(unrolled-vmap2 'vcenter (λ (x1 x2) (* 1/2 (+ x1 x2))) mins maxs)]))
|
||||
|
||||
(define (vregular-sublists vs)
|
||||
(define res
|
||||
|
|
|
|||
|
|
@ -15,7 +15,9 @@
|
|||
"common/format.rkt"
|
||||
"common/sample.rkt"
|
||||
"common/draw.rkt"
|
||||
"common/date-time.rkt")
|
||||
"common/date-time.rkt"
|
||||
"common/marching-squares.rkt"
|
||||
"common/marching-cubes.rkt")
|
||||
|
||||
(provide (only-doc-out
|
||||
(combine-out (all-from-out "common/parameters.rkt")
|
||||
|
|
@ -28,7 +30,9 @@
|
|||
(all-from-out "common/format.rkt")
|
||||
(all-from-out "common/sample.rkt")
|
||||
(all-from-out "common/draw.rkt")
|
||||
(all-from-out "common/date-time.rkt"))))
|
||||
(all-from-out "common/date-time.rkt")
|
||||
(all-from-out "common/marching-squares.rkt")
|
||||
(all-from-out "common/marching-cubes.rkt"))))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; 2D exports
|
||||
|
|
|
|||
|
|
@ -2,7 +2,6 @@
|
|||
|
||||
(require racket/class racket/match racket/list racket/flonum racket/contract racket/math
|
||||
plot/utils
|
||||
"../common/marching-cubes.rkt"
|
||||
"../common/contract-doc.rkt")
|
||||
|
||||
(provide (all-defined-out))
|
||||
|
|
@ -10,10 +9,6 @@
|
|||
;; ===================================================================================================
|
||||
;; Surfaces of constant value (isosurfaces)
|
||||
|
||||
(define (scale-normalized-polys polys xa xb ya yb za zb)
|
||||
(map (λ (poly) (scale-normalized-poly poly xa xb ya yb za zb))
|
||||
polys))
|
||||
|
||||
(define ((isosurface3d-render-proc f d samples color line-color line-width line-style alpha label)
|
||||
area)
|
||||
(match-define (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max))
|
||||
|
|
@ -48,16 +43,11 @@
|
|||
[d6 (in-vector ds10 1)]
|
||||
[d7 (in-vector ds11 1)]
|
||||
[d8 (in-vector ds11)])
|
||||
(define polys
|
||||
(heights->cube-polys
|
||||
(exact->inexact d)
|
||||
(exact->inexact d1) (exact->inexact d2) (exact->inexact d3) (exact->inexact d4)
|
||||
(exact->inexact d5) (exact->inexact d6) (exact->inexact d7) (exact->inexact d8)))
|
||||
(define polys (heights->cube-polys xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))
|
||||
|
||||
(when (not (empty? polys))
|
||||
(send area put-polygons
|
||||
(scale-normalized-polys polys xa xb ya yb za zb)
|
||||
(vcenter (list (vector xa ya za) (vector xb yb zb))))))
|
||||
(send area put-polygons polys
|
||||
(vector (* 1/2 (+ xa xb)) (* 1/2 (+ ya yb)) (* 1/2 (+ za zb))))))
|
||||
|
||||
(cond [label (rectangle-legend-entry
|
||||
label color 'solid line-color line-width line-style)]
|
||||
|
|
@ -136,16 +126,10 @@
|
|||
[d6 (in-vector ds10 1)]
|
||||
[d7 (in-vector ds11 1)]
|
||||
[d8 (in-vector ds11)])
|
||||
(define polys
|
||||
(heights->cube-polys
|
||||
(exact->inexact d)
|
||||
(exact->inexact d1) (exact->inexact d2) (exact->inexact d3) (exact->inexact d4)
|
||||
(exact->inexact d5) (exact->inexact d6) (exact->inexact d7) (exact->inexact d8)))
|
||||
|
||||
(define polys (heights->cube-polys xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))
|
||||
(when (not (empty? polys))
|
||||
(send area put-polygons
|
||||
(scale-normalized-polys polys xa xb ya yb za zb)
|
||||
(vcenter (list (vector xa ya za) (vector xb yb zb)))))))
|
||||
(send area put-polygons polys
|
||||
(vector (* 1/2 (+ xa xb)) (* 1/2 (+ ya yb)) (* 1/2 (+ za zb)))))))
|
||||
|
||||
(cond
|
||||
[label (rectangle-legend-entries
|
||||
|
|
@ -213,15 +197,10 @@
|
|||
[d7 (in-vector ds11 1)]
|
||||
[d8 (in-vector ds11)])
|
||||
(define (draw-cube xa xb ya yb za zb d1 d2 d3 d4 d5 d6 d7 d8)
|
||||
(define polys
|
||||
(heights->cube-polys
|
||||
0.0
|
||||
(exact->inexact d1) (exact->inexact d2) (exact->inexact d3) (exact->inexact d4)
|
||||
(exact->inexact d5) (exact->inexact d6) (exact->inexact d7) (exact->inexact d8)))
|
||||
(define polys (heights->cube-polys xa xb ya yb za zb 0.0 d1 d2 d3 d4 d5 d6 d7 d8))
|
||||
(when (not (empty? polys))
|
||||
(send area put-polygons
|
||||
(scale-normalized-polys polys xa xb ya yb za zb)
|
||||
(vcenter (list (vector xa ya za) (vector xb yb zb))))))
|
||||
(send area put-polygons polys
|
||||
(vector (* 1/2 (+ xa xb)) (* 1/2 (+ ya yb)) (* 1/2 (+ za zb))))))
|
||||
(cond [(and (xb . > . 0) (ya . < . 0) (yb . > . 0))
|
||||
(let* ([yb -0.00001]
|
||||
[d3 (f xb yb za)]
|
||||
|
|
@ -251,8 +230,8 @@
|
|||
(define-values (θ ρ)
|
||||
(cond [(and (fl= x 0.0) (fl= y 0.0)) (values 0.0 0.0)]
|
||||
[else (values (flmodulo (flatan2 y x) 2pi)
|
||||
(flatan (fl/ z (distance x y))))]))
|
||||
(fl- (exact->inexact (f θ ρ)) (distance x y z))))
|
||||
(flatan (fl/ z (fldistance x y))))]))
|
||||
(fl- (exact->inexact (f θ ρ)) (fldistance x y z))))
|
||||
|
||||
(defproc (polar3d [f (real? real? . -> . real?)]
|
||||
[#:x-min x-min (or/c regular-real? #f) #f]
|
||||
|
|
@ -269,7 +248,7 @@
|
|||
[#:alpha alpha (real-in 0 1) (surface-alpha)]
|
||||
[#:label label (or/c string? #f) #f]
|
||||
) renderer3d?
|
||||
(define vs (for*/list ([θ (in-list (linear-seq 0 2pi (* 4 samples)))]
|
||||
(define vs (for*/list ([θ (in-list (linear-seq 0.0 2pi (* 4 samples)))]
|
||||
[ρ (in-list (linear-seq (* -1/2 pi) (* 1/2 pi) (* 2 samples)))])
|
||||
(3d-polar->3d-cartesian θ ρ (f θ ρ))))
|
||||
(define rvs (filter vregular? vs))
|
||||
|
|
|
|||
|
|
@ -790,7 +790,7 @@
|
|||
|
||||
;; Drawing shapes
|
||||
|
||||
(define/public (put-line v1 v2 [c (vcenter (list v1 v2))])
|
||||
(define/public (put-line v1 v2 [c (v* (v+ v1 v2) 1/2)])
|
||||
(let/ec return
|
||||
(unless (and (vregular? v1) (vregular? v2)) (return (void)))
|
||||
(let-values ([(v1 v2) (if clipping?
|
||||
|
|
|
|||
|
|
@ -7,10 +7,13 @@
|
|||
@itemlist[
|
||||
@item{Planned new renderers
|
||||
@itemlist[
|
||||
@item{Functions with integer domains}
|
||||
@item{2D kernel density estimator}
|
||||
@item{3D kernel density estimator}
|
||||
@item{3D decorations: labeled points, axes, grids}
|
||||
@item{3D vector field}
|
||||
@item{2D and 3D stacked histograms}
|
||||
@item{2D grouped histograms}
|
||||
]
|
||||
}
|
||||
@item{Possible minor enhancements
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
#lang racket
|
||||
|
||||
(require plot plot/utils)
|
||||
(require plot plot/utils racket/flonum)
|
||||
|
||||
(time
|
||||
(plot3d (isosurface3d (λ (x y z) (sqrt (+ (sqr x) (sqr y) (sqr z)))) 1
|
||||
|
|
@ -41,9 +41,11 @@
|
|||
#:y-min -2 #:y-max 2
|
||||
#:z-min -2 #:z-max 2))
|
||||
|
||||
(define 2pi (* 2 pi))
|
||||
|
||||
(time
|
||||
(define (f1 θ ρ) (+ 1 (/ θ 2 pi) (* 1/8 (sin (* 8 ρ)))))
|
||||
(define (f2 θ ρ) (+ 0 (/ θ 2 pi) (* 1/8 (sin (* 8 ρ)))))
|
||||
(define (f1 θ ρ) (+ 1 (/ θ 2pi) (* 1/8 (sin (* 8 ρ)))))
|
||||
(define (f2 θ ρ) (+ (/ θ 2pi) (* 1/8 (sin (* 8 ρ)))))
|
||||
|
||||
(plot3d (list (polar3d f1 #:samples 41 #:color "navajowhite" #:line-style 'transparent #:alpha 2/3)
|
||||
(polar3d f2 #:samples 41 #:color "navajowhite" #:line-style 'transparent #:alpha 2/3))
|
||||
|
|
|
|||
|
|
@ -2,6 +2,7 @@
|
|||
|
||||
(require "common/contract.rkt"
|
||||
"common/marching-squares.rkt"
|
||||
"common/marching-cubes.rkt"
|
||||
"contracted/parameters.rkt"
|
||||
"contracted/math.rkt"
|
||||
"contracted/axis-transform.rkt"
|
||||
|
|
@ -16,6 +17,7 @@
|
|||
|
||||
(provide (all-from-out "common/contract.rkt")
|
||||
(all-from-out "common/marching-squares.rkt")
|
||||
(all-from-out "common/marching-cubes.rkt")
|
||||
(all-from-out "contracted/parameters.rkt")
|
||||
(all-from-out "contracted/math.rkt")
|
||||
(all-from-out "contracted/axis-transform.rkt")
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user