suzanne.soy/racket
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Plots with intervals too small or too large for flonums (2D and 3D)

Browse Source 
Bounds fixpoint computation now uses only exact rationals
Speed improvements (3d-plot-area% now uses flonums internally as much as possible)
This commit is contained in:
Neil Toronto 2011-11-14 16:15:12 -07:00
parent 20e09c9f7f
commit 75f76986c3
23 changed files with 860 additions and 516 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
collects/plot/common/draw.rkt
View File

@ -279,8 +279,6 @@
;; As long as this process is monotone and bounded, the distance off the dc is zero in the limit. In
;; practice, only a few iterations drives this distance to less than 1 drawing unit.
(define (appx= x y) ((abs (- x y)) . < . 1/2))
(define (margin-fixpoint x-min x-max y-min y-max
init-left init-right init-top init-bottom
get-vs)
@ -288,21 +286,20 @@
(for/fold ([left init-left] [right init-right] [top init-top] [bottom init-bottom]
) ([i (in-range 3)])
(match-define (list (vector xs ys) ...) (get-vs left right top bottom))
(define param-x-min (apply min x-min xs))
(define param-x-max (apply max (sub1 x-max) xs))
(define param-y-min (apply min y-min ys))
(define param-y-max (apply max (sub1 y-max) ys))
(define new-left (+ left (- x-min param-x-min)))
(define new-right (- right (- (sub1 x-max) param-x-max)))
(define new-top (+ top (- y-min param-y-min)))
(define new-bottom (- bottom (- (sub1 y-max) param-y-max)))
(define new-left (round (+ left (- x-min param-x-min))))
(define new-right (round (- right (- (sub1 x-max) param-x-max))))
(define new-top (round (+ top (- y-min param-y-min))))
(define new-bottom (round (- bottom (- (sub1 y-max) param-y-max))))
;; Early out: if the margins haven't changed much, another iteration won't change them more
;; (hopefully)
(when (and (appx= left new-left) (appx= right new-right)
(appx= top new-top) (appx= bottom new-bottom))
(when (and (= left new-left) (= right new-right)
(= top new-top) (= bottom new-bottom))
(return new-left new-right new-top new-bottom))
(values new-left new-right new-top new-bottom))))
@ -380,30 +377,30 @@
(match-define (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max)) r)
(list
;; Top (z-max) face
(list (vector 0 0 1)
(list (vector 0.0 0.0 1.0)
(vector x-min y-min z-max) (vector x-max y-min z-max)
(vector x-max y-max z-max) (vector x-min y-max z-max))
;; Front (y-min) face
(if ((cos theta) . > . 0)
(list (vector 0 -1 0)
(if ((cos theta) . > . 0.0)
(list (vector 0.0 -1.0 0.0)
(vector x-min y-min z-min) (vector x-max y-min z-min)
(vector x-max y-min z-max) (vector x-min y-min z-max))
empty)
;; Back (y-max) face
(if ((cos theta) . < . 0)
(list (vector 0 1 0)
(if ((cos theta) . < . 0.0)
(list (vector 0.0 1.0 0.0)
(vector x-min y-max z-min) (vector x-max y-max z-min)
(vector x-max y-max z-max) (vector x-min y-max z-max))
empty)
;; Left (x-min) face
(if ((sin theta) . > . 0)
(list (vector -1 0 0)
(if ((sin theta) . > . 0.0)
(list (vector -1.0 0.0 0.0)
(vector x-min y-min z-min) (vector x-min y-max z-min)
(vector x-min y-max z-max) (vector x-min y-min z-max))
empty)
;; Right (x-max) face
(if ((sin theta) . < . 0)
(list (vector 1 0 0)
(if ((sin theta) . < . 0.0)
(list (vector 1.0 0.0 0.0)
(vector x-max y-min z-min) (vector x-max y-max z-min)
(vector x-max y-max z-max) (vector x-max y-min z-max))
empty)))

113
collects/plot/common/marching-cubes.rkt
View File

@ -3,6 +3,7 @@
(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"
"utils.rkt"
"marching-utils.rkt"
"contract-doc.rkt")
@ -16,20 +17,20 @@
;; edge vertexes
(define-syntax-rule (edge-1-2 d d1 d2) (vector (solve-t d d1 d2) 0.0 0.0))
(define-syntax-rule (edge-2-3 d d2 d3) (vector 1.0 (solve-t d d2 d3) 0.0))
(define-syntax-rule (edge-3-4 d d3 d4) (vector (solve-t d d4 d3) 1.0 0.0))
(define-syntax-rule (edge-1-4 d d1 d4) (vector 0.0 (solve-t d d1 d4) 0.0))
(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-5-6 d d5 d6) (vector (solve-t d d5 d6) 0.0 1.0))
(define-syntax-rule (edge-6-7 d d6 d7) (vector 1.0 (solve-t d d6 d7) 1.0))
(define-syntax-rule (edge-7-8 d d7 d8) (vector (solve-t d d7 d8) 1.0 1.0))
(define-syntax-rule (edge-5-8 d d5 d8) (vector 0.0 (solve-t d d5 d8) 1.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-1-5 d d1 d5) (vector 0.0 0.0 (solve-t d d1 d5)))
(define-syntax-rule (edge-2-6 d d2 d6) (vector 1.0 0.0 (solve-t d d2 d6)))
(define-syntax-rule (edge-3-7 d d3 d7) (vector 1.0 1.0 (solve-t d d3 d7)))
(define-syntax-rule (edge-4-8 d d4 d8) (vector 0.0 1.0 (solve-t d d4 d8)))
(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)))
#|
Cube vertex numbers:
@ -74,10 +75,10 @@ Cube vertex numbers:
(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)
(edge-3-7 d d3 d7) (edge-1-5 d d1 d5))
(list (edge-3-4 d d3 d4) (edge-1-4 d d1 d4)
(edge-1-5 d d1 d5) (edge-3-7 d d3 d7)))]
(list (list (edge-1-2 d d1 d2) (edge-2-3 d d2 d3)
(edge-3-7 d d3 d7) (edge-1-5 d d1 d5))
(list (edge-3-4 d d3 d4) (edge-1-4 d d1 d4)
(edge-1-5 d d1 d5) (edge-3-7 d d3 d7)))]
[else
(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)))]))
@ -104,7 +105,7 @@ Cube vertex numbers:
(list (list (edge-1-5 d d1 d5) (edge-2-6 d d2 d6)
(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<))
@ -510,32 +511,62 @@ 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 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))
(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))))))
(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))
(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)))))]))

445
collects/plot/common/marching-squares.rkt
View File

@ -3,6 +3,7 @@
(require racket/flonum racket/fixnum racket/list racket/match racket/unsafe/ops racket/contract
(for-syntax racket/base)
"math.rkt"
"utils.rkt"
"marching-utils.rkt"
"contract-doc.rkt")
@ -58,23 +59,23 @@ above.
;(: lines1000 (FType Lines))
(define-syntax-rule (lines1000 z z1 z2 z3 z4)
(list (vector (solve-t z z1 z2) 0.0
0.0 (solve-t z z1 z4))))
(list (vector (unsafe-solve-t z z1 z2) 0.0
0.0 (unsafe-solve-t z z1 z4))))
;(: lines0100 (FType Lines))
(define-syntax-rule (lines0100 z z1 z2 z3 z4)
(list (vector (solve-t z z1 z2) 0.0
1.0 (solve-t z z2 z3))))
(list (vector (unsafe-solve-t z z1 z2) 0.0
1.0 (unsafe-solve-t z z2 z3))))
;(: lines0010 (FType Lines))
(define-syntax-rule (lines0010 z z1 z2 z3 z4)
(list (vector 1.0 (solve-t z z2 z3)
(solve-t z z4 z3) 1.0)))
(list (vector 1.0 (unsafe-solve-t z z2 z3)
(unsafe-solve-t z z4 z3) 1.0)))
;(: lines0001 (FType Lines))
(define-syntax-rule (lines0001 z z1 z2 z3 z4)
(list (vector 0.0 (solve-t z z1 z4)
(solve-t z z4 z3) 1.0)))
(list (vector 0.0 (unsafe-solve-t z z1 z4)
(unsafe-solve-t z z4 z3) 1.0)))
(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))
@ -86,13 +87,13 @@ above.
;(: lines1100 (FType Lines))
(define-syntax-rule (lines1100 z z1 z2 z3 z4)
(list (vector 0.0 (solve-t z z1 z4)
1.0 (solve-t z z2 z3))))
(list (vector 0.0 (unsafe-solve-t z z1 z4)
1.0 (unsafe-solve-t z z2 z3))))
;(: lines0110 (FType Lines))
(define-syntax-rule (lines0110 z z1 z2 z3 z4)
(list (vector (solve-t z z1 z2) 0.0
(solve-t z z4 z3) 1.0)))
(list (vector (unsafe-solve-t z z1 z2) 0.0
(unsafe-solve-t z z4 z3) 1.0)))
(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))
@ -105,64 +106,83 @@ above.
; disambiguate using average of corners as guess for center value
(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))
(vector 0.0 (solve-t z z1 z4)
(solve-t z z4 z3) 1.0))
(list (vector (solve-t z z1 z2) 0.0
0.0 (solve-t z z1 z4))
(vector 1.0 (solve-t z z2 z3)
(solve-t z z4 z3) 1.0)))))
(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)))))
(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 (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))
(if p1
(if p2
(if p3
(if p4
(lines1111 z z1 z2 z3 z4)
(lines1110 z z1 z2 z3 z4))
(if p4
(lines1101 z z1 z2 z3 z4)
(lines1100 z z1 z2 z3 z4)))
(if p3
(if p4
(lines1011 z z1 z2 z3 z4)
(lines1010 z z1 z2 z3 z4))
(if p4
(lines1001 z z1 z2 z3 z4)
(lines1000 z z1 z2 z3 z4))))
(if p2
(if p3
(if p4
(lines0111 z z1 z2 z3 z4)
(lines0110 z z1 z2 z3 z4))
(if p4
(lines0101 z z1 z2 z3 z4)
(lines0100 z z1 z2 z3 z4)))
(if p3
(if p4
(lines0011 z z1 z2 z3 z4)
(lines0010 z z1 z2 z3 z4))
(if p4
(lines0001 z z1 z2 z3 z4)
(lines0000 z z1 z2 z3 z4))))))
(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?))
(check-all-real! 'heights->lines xa xb ya yb z z1 z2 z3 z4)
(let-exact->inexact
(xa xb ya yb 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 lines
(if p1
(if p2
(if p3
(if p4
(lines1111 z z1 z2 z3 z4)
(lines1110 z z1 z2 z3 z4))
(if p4
(lines1101 z z1 z2 z3 z4)
(lines1100 z z1 z2 z3 z4)))
(if p3
(if p4
(lines1011 z z1 z2 z3 z4)
(lines1010 z z1 z2 z3 z4))
(if p4
(lines1001 z z1 z2 z3 z4)
(lines1000 z z1 z2 z3 z4))))
(if p2
(if p3
(if p4
(lines0111 z z1 z2 z3 z4)
(lines0110 z z1 z2 z3 z4))
(if p4
(lines0101 z z1 z2 z3 z4)
(lines0100 z z1 z2 z3 z4)))
(if p3
(if p4
(lines0011 z z1 z2 z3 z4)
(lines0010 z z1 z2 z3 z4))
(if p4
(lines0001 z z1 z2 z3 z4)
(lines0000 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)))))
(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))]
[(= 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))))]))
;; =============================================================================
;; Isoband marching squares: polygonizes contour between two isoline values
@ -183,16 +203,16 @@ above.
;; all corners same
(define (polys0000 za zb z1 z2 z3 z4) empty)
(define (polys1111 za zb z1 z2 z3 z4) (list 'full))
(define (polys2222 za zb z1 z2 z3 z4) empty)
(define (polys1111 za zb z1 z2 z3 z4) (list 'full))
;; -----------------------------------------------------------------------------
;; single triangle
(define (polys1000 za zb z1 z2 z3 z4)
(list (list (vector 0.0 0.0 z1)
(vector (solve-t za z1 z2) 0.0 za)
(vector 0.0 (solve-t za z1 z4) za))))
(vector (unsafe-solve-t za z1 z2) 0.0 za)
(vector 0.0 (unsafe-solve-t za z1 z4) za))))
(define polys0100 (rotate-facet polys1000))
(define polys0010 (rotate-facet polys0100))
@ -200,8 +220,8 @@ above.
(define (polys1222 za zb z1 z2 z3 z4)
(list (list (vector 0.0 0.0 z1)
(vector (solve-t zb z1 z2) 0.0 zb)
(vector 0.0 (solve-t zb z1 z4) zb))))
(vector (unsafe-solve-t zb z1 z2) 0.0 zb)
(vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
(define polys2122 (rotate-facet polys1222))
(define polys2212 (rotate-facet polys2122))
@ -211,20 +231,20 @@ above.
;; single trapezoid
(define (polys2000 za zb z1 z2 z3 z4)
(list (list (vector (solve-t zb z1 z2) 0.0 zb)
(vector (solve-t za z1 z2) 0.0 za)
(vector 0.0 (solve-t za z1 z4) za)
(vector 0.0 (solve-t zb z1 z4) zb))))
(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))))
(define polys0200 (rotate-facet polys2000))
(define polys0020 (rotate-facet polys0200))
(define polys0002 (rotate-facet polys0020))
(define (polys0222 za zb z1 z2 z3 z4)
(list (list (vector (solve-t za z1 z2) 0.0 za)
(vector (solve-t zb z1 z2) 0.0 zb)
(vector 0.0 (solve-t zb z1 z4) zb)
(vector 0.0 (solve-t za z1 z4) za))))
(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))))
(define polys2022 (rotate-facet polys0222))
(define polys2202 (rotate-facet polys2022))
@ -236,8 +256,8 @@ above.
(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 (solve-t za z2 z3) za)
(vector 0.0 (solve-t za z1 z4) za))))
(vector 1.0 (unsafe-solve-t za z2 z3) za)
(vector 0.0 (unsafe-solve-t za z1 z4) za))))
(define polys0110 (rotate-facet polys1100))
(define polys0011 (rotate-facet polys0110))
@ -246,18 +266,18 @@ above.
(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 (solve-t zb z2 z3) zb)
(vector 0.0 (solve-t zb z1 z4) zb))))
(vector 1.0 (unsafe-solve-t zb z2 z3) zb)
(vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
(define polys2112 (rotate-facet polys1122))
(define polys2211 (rotate-facet polys2112))
(define polys1221 (rotate-facet polys2211))
(define (polys0022 za zb z1 z2 z3 z4)
(list (list (vector 0.0 (solve-t za z1 z4) za)
(vector 1.0 (solve-t za z2 z3) za)
(vector 1.0 (solve-t zb z2 z3) zb)
(vector 0.0 (solve-t zb z1 z4) zb))))
(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))))
(define polys2002 (rotate-facet polys0022))
(define polys2200 (rotate-facet polys2002))
@ -267,22 +287,22 @@ above.
;; single pentagon
(define (polys0111 za zb z1 z2 z3 z4)
(list (list (vector (solve-t za z1 z2) 0.0 za)
(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 (solve-t za z1 z4) za))))
(vector 0.0 (unsafe-solve-t za z1 z4) za))))
(define polys1011 (rotate-facet polys0111))
(define polys1101 (rotate-facet polys1011))
(define polys1110 (rotate-facet polys1101))
(define (polys2111 za zb z1 z2 z3 z4)
(list (list (vector (solve-t zb z1 z2) 0.0 zb)
(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 (solve-t zb z1 z4) zb))))
(vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
(define polys1211 (rotate-facet polys2111))
(define polys1121 (rotate-facet polys1211))
@ -290,10 +310,10 @@ above.
(define (polys1002 za zb z1 z2 z3 z4)
(list (list (vector 0.0 0.0 z1)
(vector (solve-t za z1 z2) 0.0 za)
(vector (solve-t za z4 z3) 1.0 za)
(vector (solve-t zb z4 z3) 1.0 zb)
(vector 0.0 (solve-t zb z1 z4) zb))))
(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))))
(define polys2100 (rotate-facet polys1002))
(define polys0210 (rotate-facet polys2100))
@ -301,10 +321,10 @@ above.
(define (polys1220 za zb z1 z2 z3 z4)
(list (list (vector 0.0 0.0 z1)
(vector (solve-t zb z1 z2) 0.0 zb)
(vector (solve-t zb z4 z3) 1.0 zb)
(vector (solve-t za z4 z3) 1.0 za)
(vector 0.0 (solve-t za z1 z4) za))))
(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))))
(define polys0122 (rotate-facet polys1220))
(define polys2012 (rotate-facet polys0122))
@ -312,10 +332,10 @@ above.
(define (polys1200 za zb z1 z2 z3 z4)
(list (list (vector 0.0 0.0 z1)
(vector (solve-t zb z1 z2) 0.0 zb)
(vector 1.0 (solve-t zb z2 z3) zb)
(vector 1.0 (solve-t za z2 z3) za)
(vector 0.0 (solve-t za z1 z4) 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 0.0 (unsafe-solve-t za z1 z4) za))))
(define polys0120 (rotate-facet polys1200))
(define polys0012 (rotate-facet polys0120))
@ -323,10 +343,10 @@ above.
(define (polys1022 za zb z1 z2 z3 z4)
(list (list (vector 0.0 0.0 z1)
(vector (solve-t za z1 z2) 0.0 za)
(vector 1.0 (solve-t za z2 z3) za)
(vector 1.0 (solve-t zb z2 z3) zb)
(vector 0.0 (solve-t zb z1 z4) zb))))
(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))))
(define polys2102 (rotate-facet polys1022))
(define polys2210 (rotate-facet polys2102))
@ -336,36 +356,36 @@ above.
;; single hexagon
(define (polys0112 za zb z1 z2 z3 z4)
(list (list (vector (solve-t za z1 z2) 0.0 za)
(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 (solve-t zb z4 z3) 1.0 zb)
(vector 0.0 (solve-t zb z1 z4) zb)
(vector 0.0 (solve-t za z1 z4) 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))))
(define polys2011 (rotate-facet polys0112))
(define polys1201 (rotate-facet polys2011))
(define polys1120 (rotate-facet polys1201))
(define (polys2110 za zb z1 z2 z3 z4)
(list (list (vector (solve-t zb z1 z2) 0.0 zb)
(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 (solve-t za z4 z3) 1.0 za)
(vector 0.0 (solve-t za z1 z4) za)
(vector 0.0 (solve-t zb z1 z4) zb))))
(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))))
(define polys0211 (rotate-facet polys2110))
(define polys1021 (rotate-facet polys0211))
(define polys1102 (rotate-facet polys1021))
(define (polys0121 za zb z1 z2 z3 z4)
(list (list (vector (solve-t za z1 z2) 0.0 za)
(list (list (vector (unsafe-solve-t za z1 z2) 0.0 za)
(vector 1.0 0.0 z2)
(vector 1.0 (solve-t zb z2 z3) zb)
(vector (solve-t zb z4 z3) 1.0 zb)
(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 (solve-t za z1 z4) za))))
(vector 0.0 (unsafe-solve-t za z1 z4) za))))
(define polys1012 (rotate-facet polys0121))
(define polys2101 (rotate-facet polys1012))
@ -376,19 +396,19 @@ above.
(define (polys10100 za zb z1 z2 z3 z4)
(list (list (vector 0.0 0.0 z1)
(vector (solve-t za z1 z2) 0.0 za)
(vector 0.0 (solve-t za z1 z4) za))
(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 (solve-t za z4 z3) 1.0 za)
(vector 1.0 (solve-t za z2 z3) za))))
(vector (unsafe-solve-t za z4 z3) 1.0 za)
(vector 1.0 (unsafe-solve-t za z2 z3) za))))
(define (polys10101 za zb z1 z2 z3 z4)
(list (list (vector 0.0 0.0 z1)
(vector (solve-t za z1 z2) 0.0 za)
(vector 1.0 (solve-t za z2 z3) za)
(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 (solve-t za z4 z3) 1.0 za)
(vector 0.0 (solve-t za z1 z4) za))))
(vector (unsafe-solve-t za z4 z3) 1.0 za)
(vector 0.0 (unsafe-solve-t za z1 z4) za))))
(define (polys1010 za zb z1 z2 z3 z4)
(define z5 (unsafe-flavg4 z1 z2 z3 z4))
@ -400,19 +420,19 @@ above.
(define (polys1212-2 za zb z1 z2 z3 z4)
(list (list (vector 0.0 0.0 z1)
(vector (solve-t zb z1 z2) 0.0 zb)
(vector 0.0 (solve-t zb z1 z4) zb))
(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 (solve-t zb z4 z3) 1.0 zb)
(vector 1.0 (solve-t zb z2 z3) zb))))
(vector (unsafe-solve-t zb z4 z3) 1.0 zb)
(vector 1.0 (unsafe-solve-t zb z2 z3) zb))))
(define (polys1212-1 za zb z1 z2 z3 z4)
(list (list (vector 0.0 0.0 z1)
(vector (solve-t zb z1 z2) 0.0 zb)
(vector 1.0 (solve-t zb z2 z3) zb)
(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 (solve-t zb z4 z3) 1.0 zb)
(vector 0.0 (solve-t zb z1 z4) zb))))
(vector (unsafe-solve-t zb z4 z3) 1.0 zb)
(vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
(define (polys1212 za zb z1 z2 z3 z4)
(define z5 (unsafe-flavg4 z1 z2 z3 z4))
@ -426,22 +446,22 @@ above.
;; 7-sided saddle
(define (polys0212-1 za zb z1 z2 z3 z4)
(list (list (vector (solve-t za z1 z2) 0.0 za)
(vector (solve-t zb z1 z2) 0.0 zb)
(vector 1.0 (solve-t zb z2 z3) zb)
(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 (solve-t zb z4 z3) 1.0 zb)
(vector 0.0 (solve-t zb z1 z4) zb)
(vector 0.0 (solve-t za z1 z4) 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))))
(define (polys0212-2 za zb z1 z2 z3 z4)
(list (list (vector (solve-t za z1 z2) 0.0 za)
(vector (solve-t zb z1 z2) 0.0 zb)
(vector 0.0 (solve-t zb z1 z4) zb)
(vector 0.0 (solve-t za z1 z4) za))
(list (vector 1.0 (solve-t zb z2 z3) zb)
(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 (solve-t zb z4 z3) 1.0 zb))))
(vector (unsafe-solve-t zb z4 z3) 1.0 zb))))
(define (polys0212 za zb z1 z2 z3 z4)
(define z5 (unsafe-flavg4 z1 z2 z3 z4))
@ -454,22 +474,22 @@ above.
(define polys2120 (rotate-facet polys1202))
(define (polys2010-1 za zb z1 z2 z3 z4)
(list (list (vector (solve-t zb z1 z2) 0.0 zb)
(vector (solve-t za z1 z2) 0.0 za)
(vector 1.0 (solve-t za z2 z3) za)
(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 (solve-t za z4 z3) 1.0 za)
(vector 0.0 (solve-t za z1 z4) za)
(vector 0.0 (solve-t zb z1 z4) zb))))
(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))))
(define (polys2010-0 za zb z1 z2 z3 z4)
(list (list (vector (solve-t zb z1 z2) 0.0 zb)
(vector (solve-t za z1 z2) 0.0 za)
(vector 0.0 (solve-t za z1 z4) za)
(vector 0.0 (solve-t zb z1 z4) zb))
(list (vector 1.0 (solve-t za z2 z3) za)
(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 (solve-t za z4 z3) 1.0 za))))
(vector (unsafe-solve-t za z4 z3) 1.0 za))))
(define (polys2010 za zb z1 z2 z3 z4)
(define z5 (unsafe-flavg4 z1 z2 z3 z4))
@ -485,34 +505,34 @@ above.
;; 8-sided saddle
(define (polys0202-0 za zb z1 z2 z3 z4)
(list (list (vector (solve-t za z1 z2) 0.0 za)
(vector (solve-t zb z1 z2) 0.0 zb)
(vector 1.0 (solve-t zb z2 z3) zb)
(vector 1.0 (solve-t za z2 z3) za))
(list (vector 0.0 (solve-t za z1 z4) za)
(vector (solve-t za z4 z3) 1.0 za)
(vector (solve-t zb z4 z3) 1.0 zb)
(vector 0.0 (solve-t zb z1 z4) zb))))
(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))))
(define (polys0202-1 za zb z1 z2 z3 z4)
(list (list (vector (solve-t za z1 z2) 0.0 za)
(vector (solve-t zb z1 z2) 0.0 zb)
(vector 1.0 (solve-t zb z2 z3) zb)
(vector 1.0 (solve-t za z2 z3) za)
(vector (solve-t za z4 z3) 1.0 za)
(vector (solve-t zb z4 z3) 1.0 zb)
(vector 0.0 (solve-t zb z1 z4) zb)
(vector 0.0 (solve-t za z1 z4) za))))
(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))))
(define (polys0202-2 za zb z1 z2 z3 z4)
(list (list (vector (solve-t za z1 z2) 0.0 za)
(vector (solve-t zb z1 z2) 0.0 zb)
(vector 0.0 (solve-t zb z1 z4) zb)
(vector 0.0 (solve-t za z1 z4) za))
(list (vector 1.0 (solve-t zb z2 z3) zb)
(vector 1.0 (solve-t za z2 z3) za)
(vector (solve-t za z4 z3) 1.0 za)
(vector (solve-t zb z4 z3) 1.0 zb))))
(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))))
(define (polys0202 za zb z1 z2 z3 z4)
(define z5 (unsafe-flavg4 z1 z2 z3 z4))
@ -563,26 +583,53 @@ 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 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))
(unsafe-fx* t3 3))
t4))
(define f (vector-ref polys-dispatch-table facet-num))
(f za zb z1 z2 z3 z4))
(defproc (heights->polys [xa real?] [xb real?] [ya real?] [yb real?]
[za real?] [zb real?]
[z1 real?] [z2 real?] [z3 real?] [z4 real?]
) (listof (vector/c real? real? real?))
(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)))
(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 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))
(unsafe-fx* t3 3))
t4))
(define f (vector-ref polys-dispatch-table facet-num))
(for/list ([poly (in-list (f za zb z1 z2 z3 z4))])
(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))]))))
(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 (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))))])))]))

29
collects/plot/common/marching-utils.rkt
View File

@ -1,6 +1,7 @@
#lang racket/base
(require (for-syntax racket/base) racket/unsafe/ops)
(require (for-syntax racket/base) racket/unsafe/ops
"utils.rkt")
(provide (all-defined-out))
@ -9,24 +10,24 @@
;; 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)
(define-syntax-rule (unsafe-solve-t z za zb)
(unsafe-fl/ (unsafe-fl- z za) (unsafe-fl- zb za)))
(define-syntax-rule (unsolve-t za zb t)
(define-syntax-rule (unsafe-unsolve-t za zb t)
(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))))
(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 (all pred? obj ...) (and (pred? obj) ...))
(define-syntax-rule (let-exact->inexact (id ...) body0 body ...)
(let ([id (exact->inexact id)] ...)
body0 body ...))
(define-syntax (find-failure-index stx)
(syntax-case stx ()
[(_ pred? id ...)
(with-syntax ([(i ...) (build-list (length (syntax->list #'(id ...))) values)])
(syntax/loc stx
(or (and (not (pred? id)) i) ...)))]))

100
collects/plot/common/math.rkt
View File

@ -54,12 +54,79 @@
[xs (cond [(not (andmap flonum? xs)) (raise-type-error 'fldistance "flonums" xs)]
[else (unsafe-flsqrt (flsum (λ (x) (unsafe-fl* x x)) xs))])]))
(define (flonum->bit-field x)
(integer-bytes->integer (real->floating-point-bytes x 8) #f))
(define (bit-field->flonum i)
(floating-point-bytes->real (integer->integer-bytes i 8 #f)))
(defproc (flonum->ordinal [x flonum?]) integer?
(cond [(x . < . 0) (- (flonum->bit-field (- x)))]
[else (flonum->bit-field (abs x))]))
(defproc (ordinal->flonum [i (integer-in #x-7FFFFFFFFFFFFFFF #x7FFFFFFFFFFFFFFF)]) flonum?
(cond [(i . < . 0) (- (bit-field->flonum (- i)))]
[else (bit-field->flonum i)]))
(define +inf-ordinal (flonum->ordinal +inf.0))
(define -inf-ordinal (flonum->ordinal -inf.0))
(defproc (flstep [x flonum?] [n exact-integer?]) flonum?
(cond [(eqv? x +nan.0) +nan.0]
[(and (eqv? x +inf.0) (n . >= . 0)) +inf.0]
[(and (eqv? x -inf.0) (n . <= . 0)) -inf.0]
[else
(define i (+ n (flonum->ordinal x)))
(cond [(i . < . -inf-ordinal) -inf.0]
[(i . > . +inf-ordinal) +inf.0]
[else (ordinal->flonum i)])]))
(defproc (flnext [x flonum?]) flonum? #:document-body
(flstep x 1))
(defproc (flprev [x flonum?]) flonum? #:document-body
(flstep x -1))
(defthing +min.0 flonum? #:document-value (flnext 0.0))
(defthing -min.0 flonum? #:document-value (flprev 0.0))
(defthing +max.0 flonum? #:document-value (flprev +inf.0))
(defthing -max.0 flonum? #:document-value (flnext -inf.0))
;; ===================================================================================================
;; Reals
(defproc (regular-real? [x any/c]) boolean?
(and (real? x) (not (special-real? x))))
(defproc (maybe-inexact->exact [x (or/c regular-real? #f)]) (or/c regular-real? #f)
(cond [x (unless (regular-real? x)
(raise-type-error 'maybe-inexact->exact "regular real or #f" x))
(inexact->exact x)]
[else #f]))
(defproc (flonum-ok-for-range? [x-min regular-real?] [x-max regular-real?]
[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 (regular-real? inexact-x-min) (return #f))
(define inexact-x-max (exact->inexact x-max))
(unless (regular-real? inexact-x-max) (return #f))
(define inexact-x-max-prev (flprev inexact-x-max))
(unless (regular-real? inexact-x-max-prev) (return #f))
(define inexact-x-min-next (flnext inexact-x-min))
(unless (regular-real? 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)))))
(define equal?*
(case-lambda
[() #t]
@ -154,25 +221,23 @@
(cond [(not (and (exact-integer? b) (b . >= . 2))) (raise-type-error 'floor-log/base
"exact integer >= 2" 0 b x)]
[(not (and (real? x) (x . > . 0))) (raise-type-error 'floor-log/base "real > 0" 1 b x)]
[else (define y (inexact->exact (floor (/ (log x) (log b)))))
(cond [(exact? x)
(let loop ([y y] [x (/ x (expt b y))])
[else (cond [(exact? x)
(let loop ([y 0] [x x])
(cond [(x . >= . b) (loop (add1 y) (/ x b))]
[(x . < . 1) (loop (sub1 y) (* x b))]
[else y]))]
[else y])]))
[else (inexact->exact (floor (/ (log x) (log b))))])]))
(defproc (ceiling-log/base [b (and/c exact-integer? (>=/c 2))] [x (>/c 0)]) exact-integer?
(cond [(not (and (exact-integer? b) (b . >= . 2))) (raise-type-error 'floor-log/base
"exact integer >= 2" 0 b x)]
[(not (and (real? x) (x . > . 0))) (raise-type-error 'floor-log/base "real > 0" 1 b x)]
[else (define y (inexact->exact (ceiling (/ (log x) (log b)))))
(cond [(exact? x)
(let loop ([y y] [x (/ x (expt b y))])
[else (cond [(exact? x)
(let loop ([y 0] [x x])
(cond [(x . > . 1) (loop (add1 y) (/ x b))]
[(x . <= . (/ 1 b)) (loop (sub1 y) (* x b))]
[else y]))]
[else y])]))
[else (inexact->exact (ceiling (/ (log x) (log b))))])]))
(defproc (polar->cartesian [θ real?] [r real?]) (vector/c real? real?)
(cond [(not (real? θ)) (raise-type-error 'polar->cartesian "real number" 0 θ r)]
@ -367,14 +432,12 @@
(defproc (vregular? [v (vectorof real?)]) boolean?
(match v
[(vector (? real? x) (? real? y))
(cond [(flonum? x) (unsafe-flregular? x)]
[(flonum? y) (unsafe-flregular? y)]
[else #t])]
(not (or (and (flonum? x) (unsafe-flspecial? x))
(and (flonum? y) (unsafe-flspecial? y))))]
[(vector (? real? x) (? real? y) (? real? z))
(cond [(flonum? x) (unsafe-flregular? x)]
[(flonum? y) (unsafe-flregular? y)]
[(flonum? z) (unsafe-flregular? z)]
[else #t])]
(not (or (and (flonum? x) (unsafe-flspecial? x))
(and (flonum? y) (unsafe-flspecial? y))
(and (flonum? z) (unsafe-flspecial? z))))]
[_ (cond [(vector-andmap real? v) (let/ec break
(for ([x (in-vector v)])
(when (and (flonum? x) (unsafe-flspecial? x))
@ -435,8 +498,6 @@
(cond [(and (not (null? res)) (null? (car res))) (cdr res)]
[else res]))
(define default-normal (vector 0 -1 0))
(define (remove-degenerate-edges vs)
(cond
[(empty? vs) empty]
@ -450,6 +511,8 @@
(cond [(v= last (first vs)) (rest vs)]
[else vs]))]))
(define default-normal (vector 0.0 -1.0 0.0))
(define (vnormal vs)
(let ([vs (remove-degenerate-edges vs)])
(cond
@ -515,7 +578,8 @@
(defproc (ivl-inexact->exact [i ivl?]) ivl?
(match-define (ivl a b) i)
(ivl (inexact->exact a) (inexact->exact b)))
(ivl (and a (if (nan? a) a (inexact->exact a)))
(and b (if (nan? b) b (inexact->exact b)))))
(defproc (ivl-contains? [i ivl?] [x real?]) boolean?
(match-define (ivl a b) i)

13
collects/plot/common/plot-device.rkt
View File

@ -7,6 +7,7 @@
;; It is up to callers to transform view or plot coordinates into dc coordinates.
(require racket/draw racket/class racket/match racket/math racket/bool racket/list racket/contract
racket/vector
"contract.rkt"
"draw.rkt"
"math.rkt"
@ -308,9 +309,11 @@
(draw-text/anchor dc str x y anchor #t 0 angle)))
(define/public (get-text-corners str v [anchor 'top-left] [angle 0])
(when (vregular? v)
(match-define (vector x y) v)
(get-text-corners/anchor dc str x y anchor #t 0 angle)))
(cond [(vregular? v)
(match-define (vector x y) v)
(map (λ (v) (vector-map inexact->exact v))
(get-text-corners/anchor dc str x y anchor #t 0 angle))]
[else empty]))
(define/public (draw-arrow v1 v2)
(when (and (vregular? v1) (vregular? v2))
@ -370,8 +373,8 @@
(define/public (get-tick-endpoints v r angle)
(match-define (vector x y) v)
(define dx (* (cos angle) r))
(define dy (* (sin angle) r))
(define dx (* (inexact->exact (cos angle)) r))
(define dy (* (inexact->exact (sin angle)) r))
(list (vector (- x dx) (- y dy)) (vector (+ x dx) (+ y dy))))
(define/public (make-draw-tick r angle)

63
collects/plot/common/plot-element.rkt
View File

@ -84,30 +84,34 @@
;; More precisely, starting with the given plot bounds, it attempts to compute a fixpoint of
;; (apply-bounds* elems), overridden at every iteration by the plot bounds (if given). Because a
;; fixpoint doesn't always exist, or only exists in the limit, it stops after max-iters.
(define (bounds-fixpoint elems plot-bounds-rect [max-iters 4])
(define (bounds-fixpoint elems given-bounds-rect [max-iters 4])
(define res
(let/ec break
;; Shortcut eval: if the plot bounds are all known, the code below just returns them anyway
(when (rect-known? plot-bounds-rect) (break plot-bounds-rect))
;; Objective: find the fixpoint of F starting at plot-bounds-rect
(define (F bounds-rect) (rect-meet plot-bounds-rect (apply-bounds* elems bounds-rect)))
;; Iterate joint bounds to (hopefully) a fixpoint
(define-values (bounds-rect area delta-area)
(for/fold ([bounds-rect plot-bounds-rect]
[area (rect-area plot-bounds-rect)] [delta-area #f]
) ([n (in-range max-iters)])
;(printf "bounds-rect = ~v~n" bounds-rect)
;; Get new bounds from the elements' bounds functions
(define new-bounds-rect (F bounds-rect))
(define new-area (rect-area new-bounds-rect))
(define new-delta-area (and area new-area (- new-area area)))
(cond
;; Shortcut eval: if the bounds haven't changed, we have a fixpoint
[(equal? bounds-rect new-bounds-rect) (break bounds-rect)]
;; If the area grew more this iteration than last, it may not converge, so stop now
[(and delta-area new-delta-area (new-delta-area . > . delta-area)) (break bounds-rect)]
;; All good - one more iteration
[else (values new-bounds-rect new-area new-delta-area)])))
bounds-rect))
;; Shortcut eval: if the plot bounds are all given, the code below just returns them anyway
(when (rect-known? given-bounds-rect) (break given-bounds-rect))
(let ([given-bounds-rect (rect-inexact->exact given-bounds-rect)])
;; Objective: find the fixpoint of F starting at given-bounds-rect
(define (F bounds-rect) (rect-meet given-bounds-rect (apply-bounds* elems bounds-rect)))
;; Iterate joint bounds to (hopefully) a fixpoint
(define-values (bounds-rect area delta-area)
(for/fold ([bounds-rect given-bounds-rect]
[area (rect-area given-bounds-rect)] [delta-area #f]
) ([n (in-range max-iters)])
;(printf "bounds-rect = ~v~n" bounds-rect)
;; Get new bounds from the elements' bounds functions
(define new-bounds-rect (F bounds-rect))
(define new-area (rect-area new-bounds-rect))
(define new-delta-area (and area new-area (- new-area area)))
(cond
;; Shortcut eval: if the bounds haven't changed, we have a fixpoint
[(equal? bounds-rect new-bounds-rect) (break bounds-rect)]
;; If the area grew more this iteration than last, it may not converge, so stop now
[(and delta-area new-delta-area (new-delta-area . > . delta-area)) (break bounds-rect)]
;; All good - one more iteration
[else (values new-bounds-rect new-area new-delta-area)])))
bounds-rect)))
;(printf "fixpoint bounds-rect = ~v~n" res)
res)
;; Applies the bounds functions of multiple plot elements, in parallel, and returns the smallest
;; bounds containing all the new bounds. This function is monotone and increasing regardless of
@ -120,7 +124,12 @@
;; what bounds will you try to use?
(define (apply-bounds elem bounds-rect)
(match-define (plot-element elem-bounds-rect elem-bounds-fun _) elem)
(let ([elem-bounds-rect (cond [elem-bounds-rect (rect-meet bounds-rect elem-bounds-rect)]
[else bounds-rect])])
(cond [elem-bounds-fun (elem-bounds-fun elem-bounds-rect)]
[else elem-bounds-rect])))
;(printf "elem-bounds-rect = ~v~n" elem-bounds-rect)
(let ([elem-bounds-rect (if elem-bounds-rect
(rect-meet bounds-rect (rect-inexact->exact elem-bounds-rect))
bounds-rect)])
(if elem-bounds-fun
(let ([new-elem-bounds-rect (elem-bounds-fun elem-bounds-rect)])
;(printf "new-elem-bounds-rect = ~v~n" new-elem-bounds-rect)
(rect-inexact->exact new-elem-bounds-rect))
elem-bounds-rect)))

58
collects/plot/common/sample.rkt
View File

@ -2,7 +2,7 @@
;; Functions that sample from functions, and functions that create memoized samplers.
(require racket/match racket/flonum racket/math racket/contract racket/list
(require racket/match racket/flonum racket/math racket/contract racket/list racket/vector
"contract-doc.rkt"
"math.rkt"
"axis-transform.rkt")
@ -75,7 +75,6 @@
(struct 2d-sample (xs ys zss z-min z-max) #:transparent)
(struct 3d-sample (xs ys zs dsss d-min d-max) #:transparent)
(defcontract sample/c (list/c (listof real?) (listof real?)))
(defcontract sampler/c (real? real? exact-nonnegative-integer? . -> . sample?))
(defcontract 2d-sampler/c (real? real? exact-nonnegative-integer?
@ -101,7 +100,9 @@
(define-values (y-min y-max)
(cond [(empty? rys) (values #f #f)]
[else (values (apply min* rys) (apply max* rys))]))
(sample xs ys y-min y-max))))))
(sample xs ys
(maybe-inexact->exact y-min)
(maybe-inexact->exact y-max)))))))
(defproc (make-2d-function->sampler [transform-x-thnk (-> axis-transform/c)]
[transform-y-thnk (-> axis-transform/c)]
@ -124,7 +125,9 @@
(unless (and z-min (z . >= . z-min)) (set! z-min z))
(unless (and z-max (z . <= . z-max)) (set! z-max z)))
z))))
(2d-sample xs ys zss z-min z-max))))))
(2d-sample xs ys zss
(maybe-inexact->exact z-min)
(maybe-inexact->exact z-max)))))))
(defproc (make-3d-function->sampler [transform-x-thnk (-> axis-transform/c)]
[transform-y-thnk (-> axis-transform/c)]
@ -153,7 +156,9 @@
(unless (and d-min (d . >= . d-min)) (set! d-min d))
(unless (and d-max (d . <= . d-max)) (set! d-max d)))
d)))))
(3d-sample xs ys zs dsss d-min d-max))))))
(3d-sample xs ys zs dsss
(maybe-inexact->exact d-min)
(maybe-inexact->exact d-max)))))))
(define-syntax-rule (for-2d-sample (xa xb ya yb z1 z2 z3 z4) sample expr ...)
(let ()
@ -199,3 +204,46 @@
(values xb d2 d3 d6 d7))
(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 regular-real?] [x-max regular-real?]
[y-min regular-real?] [y-max regular-real?]) 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 regular-real?] [x-max regular-real?]
[y-min regular-real?] [y-max regular-real?]
[z-min regular-real?] [z-max regular-real?]) 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 regular-real?] [x-max regular-real?]
[y-min regular-real?] [y-max regular-real?]
[z-min regular-real?] [z-max regular-real?]
[d-min regular-real?] [d-max regular-real?]) 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)))

4
collects/plot/common/utils.rkt
View File

@ -81,3 +81,7 @@
(values (cons (reverse lst) res)
rest-xs)))
(reverse res)))
(define-syntax-rule (let-map (id ...) fun body0 body ...)
(let ([id (fun id)] ...)
body0 body ...))

5
collects/plot/contracted/math.rkt
View File

@ -7,8 +7,11 @@
;; Flonums
nan? infinite? special-real?
flblend flatan2 flsum flmodulo fldistance
(activate-contract-out flonum->ordinal ordinal->flonum flstep flnext flprev
flonum-ok-for-range?)
-max.0 -min.0 +min.0 +max.0
;; Reals
regular-real?
regular-real? maybe-inexact->exact
min* max* degrees->radians radians->degrees blend atan2 sum real-modulo distance
floor-log/base ceiling-log/base
polar->cartesian 3d-polar->3d-cartesian

9
collects/plot/contracted/sample.rkt
View File

@ -23,7 +23,14 @@
sampler/c 2d-sampler/c 3d-sampler/c
make-function->sampler
make-2d-function->sampler
make-3d-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))])))
map*

31
collects/plot/plot2d/contour.rkt
View File

@ -21,10 +21,13 @@
(when (<= z-min z z-max)
(send area put-alpha alpha)
(send area put-pen color width style)
(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)))))
(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))))))
(cond [label (line-legend-entry label color width style)]
[else empty]))
@ -57,10 +60,13 @@
(match-define (2d-sample xs ys zss z-min z-max) sample)
(match-define (list (tick zs _ labels) ...) (contour-ticks z-min z-max levels #f))
(let ([colors (maybe-apply colors zs)]
[widths (maybe-apply widths zs)]
[styles (maybe-apply styles zs)]
[alphas (maybe-apply alphas zs)])
(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)])
(for ([z (in-list zs)]
[color (in-cycle colors)]
[width (in-cycle widths)]
@ -116,9 +122,12 @@
(values (ivl za zb) (format "[~a,~a]" la lb))))
(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)])
(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)])
(for ([za (in-list zs)]
[zb (in-list (rest zs))]
[color (in-cycle colors)]

5
collects/plot/plot2d/plot.rkt
View File

@ -35,14 +35,12 @@
(define (get-bounds-rect renderer-list x-min x-max y-min y-max)
(define given-bounds-rect (vector (ivl x-min x-max) (ivl y-min y-max)))
(define plot-bounds-rect (bounds-fixpoint renderer-list given-bounds-rect))
(when (or (not (rect-regular? plot-bounds-rect))
(rect-zero-area? plot-bounds-rect))
(match-define (vector (ivl x-min x-max) (ivl y-min y-max)) plot-bounds-rect)
(error 'plot "could not determine sensible plot bounds; got x ∈ [~a,~a], y ∈ [~a,~a]"
x-min x-max y-min y-max))
(rect-inexact->exact plot-bounds-rect))
plot-bounds-rect)
(define (get-ticks renderer-list bounds-rect)
(define-values (all-x-ticks all-x-far-ticks all-y-ticks all-y-far-ticks)
@ -51,7 +49,6 @@
(define ticks-fun (plot-element-ticks-fun r))
(cond [ticks-fun (ticks-fun bounds-rect)]
[else (values empty empty empty empty)])))
(values (remove-duplicates (append* all-x-ticks))
(remove-duplicates (append* all-x-far-ticks))
(remove-duplicates (append* all-y-ticks))

48
collects/plot/plot3d/contour.rkt
View File

@ -19,13 +19,16 @@
(when (<= z-min z z-max)
(send area put-alpha alpha)
(send area put-pen color width style)
(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)
(define center (vector (* 1/2 (+ xa xb)) (* 1/2 (+ ya yb))
(* 1/2 (+ (min z1 z2 z3 z4) (max z1 z2 z3 z4)))))
(send area put-line v1 v2 center))))
(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)
(define center (vector (* 1/2 (+ xa xb)) (* 1/2 (+ ya yb))
(* 1/2 (+ (min z1 z2 z3 z4) (max z1 z2 z3 z4)))))
(send area put-line v1 v2 center)))))
(cond [label (line-legend-entry label color width style)]
[else empty]))
@ -62,10 +65,13 @@
y-min y-max (animated-samples samples)))
(match-define (list (tick zs _ labels) ...) (contour-ticks z-min z-max levels #f))
(let ([colors (maybe-apply colors zs)]
[widths (maybe-apply widths zs)]
[styles (maybe-apply styles zs)]
[alphas (maybe-apply alphas zs)])
(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)])
(for ([z (in-list zs)]
[color (in-cycle colors)]
[width (in-cycle widths)]
@ -115,9 +121,8 @@
area)
(match-define (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max))
(send area get-bounds-rect))
(define sample (f x-min x-max (animated-samples samples) y-min y-max (animated-samples samples)))
(match-define (2d-sample xs ys zss fz-min fz-max) sample)
(define sample (f x-min x-max (animated-samples samples)
y-min y-max (animated-samples samples)))
(match-define (list (tick zs _ labels) ...) (contour-ticks z-min z-max levels #t))
(define-values (z-ivls ivl-labels)
@ -127,12 +132,15 @@
[lb (in-list (rest labels))])
(values (ivl za zb) (format "[~a,~a]" la lb))))
(let ([colors (maybe-apply colors z-ivls)]
[styles (maybe-apply styles z-ivls)]
[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)])
(let* ([colors (maybe-apply colors z-ivls)]
[styles (maybe-apply styles z-ivls)]
[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)])
(for ([za (in-list zs)]
[zb (in-list (rest zs))]
[color (in-cycle colors)]

73
collects/plot/plot3d/isosurface.rkt
View File

@ -22,12 +22,15 @@
(send area put-alpha alpha)
(send area put-brush color style)
(send area put-pen line-color line-width line-style)
(for-3d-sample
(xa xb ya yb za zb d1 d2 d3 d4 d5 d6 d7 d8) sample
(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 polys
(vector (* 1/2 (+ xa xb)) (* 1/2 (+ ya yb)) (* 1/2 (+ za zb))))))
(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
(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 polys
(vector (* 1/2 (+ xa xb)) (* 1/2 (+ ya yb)) (* 1/2 (+ za zb)))))))
(cond [label (rectangle-legend-entry
label color style line-color line-width line-style)]
@ -76,13 +79,16 @@
(match-define (list (tick ds _ labels) ...) (isosurface-ticks d-min d-max levels))
#;(define ds (linear-seq d-min d-max levels #:start? (and rd-min #t) #:end? (and rd-max #t)))
(let ([colors (maybe-apply colors ds)]
[styles (maybe-apply styles ds)]
[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)])
(for ([d (in-list ds)]
(let* ([colors (maybe-apply colors ds)]
[styles (maybe-apply styles ds)]
[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)])
(for ([d (in-list ds)]
[color (in-cycle colors)]
[style (in-cycle styles)]
[alpha (in-cycle alphas)]
@ -137,8 +143,8 @@
z-min z-max (animated-samples samples)))
(match-define (3d-sample xs ys zs dsss d-min d-max) sample)
(define (draw-cube xa xb ya yb za zb d1 d2 d3 d4 d5 d6 d7 d8)
(define polys (heights->cube-polys xa xb ya yb za zb 0.0 d1 d2 d3 d4 d5 d6 d7 d8))
(define (draw-cube xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 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 polys
(vector (* 1/2 (+ xa xb)) (* 1/2 (+ ya yb)) (* 1/2 (+ za zb))))))
@ -146,23 +152,26 @@
(send area put-alpha alpha)
(send area put-brush color style)
(send area put-pen line-color line-width line-style)
(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 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 d1 d2 d3 d4 d5 d6 d7 d8))]
[else
(draw-cube xa xb ya yb za zb d1 d2 d3 d4 d5 d6 d7 d8)]))
(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)])))
(cond [label (rectangle-legend-entry
label color style line-color line-width line-style)]

30
collects/plot/plot3d/matrix.rkt
View File

@ -7,7 +7,7 @@
(provide m3-apply m3-transpose m3* m3-rotate-z m3-rotate-x m3-scale)
(define-syntax-rule (dot x1 y1 z1 x2 y2 z2) (+ (* x1 x2) (* y1 y2) (* z1 z2)))
(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 (vector v11 v12 v13) (vector v21 v22 v23) (vector v31 v32 v33)) m)
@ -28,20 +28,22 @@
(vector (m3-apply m v1) (m3-apply m v2) (m3-apply m v3)))
(define (m3-rotate-z theta)
(define cos-theta (cos theta))
(define sin-theta (sin theta))
(vector (vector cos-theta (- sin-theta) 0.0)
(vector sin-theta cos-theta 0.0)
(vector 0.0 0.0 1.0)))
(let ([theta (exact->inexact theta)])
(define cos-theta (flcos theta))
(define sin-theta (flsin theta))
(vector (vector cos-theta (- sin-theta) 0.0)
(vector sin-theta cos-theta 0.0)
(vector 0.0 0.0 1.0))))
(define (m3-rotate-x rho)
(define cos-rho (cos rho))
(define sin-rho (sin rho))
(vector (vector 1.0 0.0 0.0)
(vector 0.0 cos-rho (- sin-rho))
(vector 0.0 sin-rho cos-rho)))
(let ([rho (exact->inexact rho)])
(define cos-rho (flcos rho))
(define sin-rho (flsin rho))
(vector (vector 1.0 0.0 0.0)
(vector 0.0 cos-rho (- sin-rho))
(vector 0.0 sin-rho cos-rho))))
(define (m3-scale x-scale y-scale z-scale)
(vector (vector x-scale 0.0 0.0)
(vector 0.0 y-scale 0.0)
(vector 0.0 0.0 z-scale)))
(vector (vector (exact->inexact x-scale) 0.0 0.0)
(vector 0.0 (exact->inexact y-scale) 0.0)
(vector 0.0 0.0 (exact->inexact z-scale))))

224
collects/plot/plot3d/plot-area.rkt
View File

@ -1,6 +1,6 @@
#lang racket/base
(require racket/class racket/match racket/list racket/math racket/contract racket/vector
(require racket/class racket/match racket/list racket/math racket/contract racket/vector racket/flonum
"../common/math.rkt"
"../common/plot-device.rkt"
"../common/ticks.rkt"
@ -9,6 +9,7 @@
"../common/axis-transform.rkt"
"../common/parameters.rkt"
"../common/sample.rkt"
"../common/utils.rkt"
"matrix.rkt"
"shape.rkt"
"clip.rkt")
@ -98,53 +99,60 @@
(define theta (+ (degrees->radians angle) 0.00001))
(define rho (degrees->radians altitude))
(define identity-transforms?
(and (equal? (plot-x-transform) id-transform)
(equal? (plot-y-transform) id-transform)
(equal? (plot-z-transform) id-transform)))
;; 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
;; 3. View coordinates (from normalized coordinates: rotate)
;; 4. Device context coordinates (from view coordinates: project to 2D)
(match-define (invertible-function fx _) (apply-axis-transform (plot-x-transform) x-min x-max))
(match-define (invertible-function fy _) (apply-axis-transform (plot-y-transform) y-min y-max))
(match-define (invertible-function fz _) (apply-axis-transform (plot-z-transform) z-min z-max))
(define axis-transform
(cond
[identity-transforms? (λ (v) v)]
[else (λ (v)
(match-define (vector x y z) v)
(vector (fx x) (fy y) (fz z)))]))
(define identity-transforms?
(and (equal? (plot-x-transform) id-transform)
(equal? (plot-y-transform) id-transform)
(equal? (plot-z-transform) id-transform)))
(define (center v)
(match-define (vector x y z) v)
(vector (- x x-mid) (- y y-mid) (- z z-mid)))
(define flonum-ok? (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max))
;; There are four coordinate systems:
;; 1. Plot coordinates (original, user-facing coordinate system)
;; 2. Normalized coordinates (axis-transformed, centered, and scaled to [-1,1] on each axis)
;; 3. View coordinates (normalized coordinates, rotated)
;; 4. DC coordinates (view coordinates, projected to 2D)
(define plot->norm
(cond [identity-transforms?
(cond [flonum-ok?
(let-map
(x-mid y-mid z-mid x-size y-size z-size) exact->inexact
(λ (v)
(match-define (vector x y z) v)
(vector (fl/ (fl- (exact->inexact x) x-mid) x-size)
(fl/ (fl- (exact->inexact y) y-mid) y-size)
(fl/ (fl- (exact->inexact z) z-mid) z-size))))]
[else
(λ (v)
(match-define (vector x y z) v)
(vector (exact->inexact (/ (- x x-mid) x-size))
(exact->inexact (/ (- y y-mid) y-size))
(exact->inexact (/ (- z z-mid) z-size))))])]
[else
(λ (v)
(match-define (vector x y z) v)
(vector (exact->inexact (/ (- (fx x) x-mid) x-size))
(exact->inexact (/ (- (fy y) y-mid) y-size))
(exact->inexact (/ (- (fz z) z-mid) z-size))))]))
;; View coordinates show up mostly in tick decorations. Normalized coordinates are only used for
;; surface normals. Most user vertexes get transformed from plot coordinates directly to DC
;; coordinates.
(define scale-matrix (m3-scale (/ x-size) (/ y-size) (/ z-size)))
(define rotate-theta-matrix (m3-rotate-z theta))
(define rotate-rho-matrix (m3-rotate-x rho))
(define rotation-matrix (m3* rotate-rho-matrix rotate-theta-matrix))
(define transform-matrix (m3* rotation-matrix scale-matrix))
(define (plot->norm v) (m3-apply scale-matrix (center (axis-transform v))))
(define (norm->view v) (m3-apply rotation-matrix v))
(define (plot->view v) (m3-apply transform-matrix (center (axis-transform v))))
(define transform-matrix/no-rho (m3* (m3-rotate-z theta) scale-matrix))
(define (plot->view/no-rho v) (m3-apply transform-matrix/no-rho (center (axis-transform v))))
(define (plot->view v) (norm->view (plot->norm v)))
(define (plot->view/no-rho v) (m3-apply rotate-theta-matrix (plot->norm v)))
(define (norm->view/no-rho v) (m3-apply rotate-theta-matrix v))
(define (rotate/rho v) (m3-apply rotate-rho-matrix v))
(define view->dc #f)
(define (plot->dc/no-axis-trans v) (view->dc (m3-apply transform-matrix (center v))))
(define (plot->dc v) (view->dc (plot->view v)))
(define (norm->dc v) (view->dc (norm->view v)))
(define-values (view-x-size view-y-size view-z-size)
(match-let ([(vector view-x-ivl view-y-ivl view-z-ivl)
@ -164,30 +172,34 @@
(define area-y-mid (* 1/2 (+ area-y-min area-y-max)))
(define area-per-view-x (/ (- area-x-max area-x-min) view-x-size))
(define area-per-view-z (/ (- area-y-max area-y-min) view-z-size))
(λ (v)
(match-define (vector x _ z) v)
(vector (+ area-x-mid (* x area-per-view-x))
(- area-y-mid (* z area-per-view-z)))))
(let-map
(area-x-mid area-y-mid area-per-view-x area-per-view-z) exact->inexact
(λ (v)
(match-define (vector x _ z) v)
(vector (fl+ area-x-mid (fl* x area-per-view-x))
(fl- area-y-mid (fl* z area-per-view-z))))))
;; Initial view->dc
(define init-top-margin (if (and (plot-decorations?) (plot-title)) (* 3/2 char-height) 0))
(set! view->dc (make-view->dc 0 0 init-top-margin 0))
(define (plot-dir->dc-angle v)
(match-define (vector dx dy)
(v- (plot->dc/no-axis-trans (v+ v (vector x-mid y-mid z-mid)))
(plot->dc/no-axis-trans (vector x-mid y-mid z-mid))))
(define (x-axis-angle)
(match-define (vector dx dy) (v- (norm->dc (vector 0.5 0.0 0.0))
(norm->dc (vector -0.5 0.0 0.0))))
(- (atan2 (- dy) dx)))
(define (x-axis-angle) (plot-dir->dc-angle #(1 0 0)))
(define (y-axis-angle) (plot-dir->dc-angle #(0 1 0)))
(define (y-axis-angle)
(match-define (vector dx dy) (v- (norm->dc (vector 0.0 0.5 0.0))
(norm->dc (vector 0.0 -0.5 0.0))))
(- (atan2 (- dy) dx)))
(define/public (plot-dir->dc-dir v)
(vnormalize (v- (plot->dc/no-axis-trans (v+ v (vector x-mid y-mid z-mid)))
(plot->dc/no-axis-trans (vector x-mid y-mid z-mid)))))
(define (x-axis-dir)
(vnormalize (v- (norm->dc (vector 0.5 0.0 0.0))
(norm->dc (vector -0.5 0.0 0.0)))))
(define (x-axis-dir) (plot-dir->dc-dir #(1 0 0)))
(define (y-axis-dir) (plot-dir->dc-dir #(0 1 0)))
(define (y-axis-dir)
(vnormalize (v- (norm->dc (vector 0.0 0.5 0.0))
(norm->dc (vector 0.0 -0.5 0.0)))))
;; ===============================================================================================
;; Tick and label constants
@ -208,6 +220,16 @@
(define z-far-axis-x (if x-axis-y-min? x-max x-min))
(define z-far-axis-y (if y-axis-x-min? y-min y-max))
(define x-axis-norm-y (if x-axis-y-min? -0.5 0.5))
(define y-axis-norm-x (if y-axis-x-min? -0.5 0.5))
(define z-axis-norm-x (if x-axis-y-min? -0.5 0.5))
(define z-axis-norm-y (if y-axis-x-min? 0.5 -0.5))
(define x-far-axis-norm-y (if x-axis-y-min? 0.5 -0.5))
(define y-far-axis-norm-x (if y-axis-x-min? 0.5 -0.5))
(define z-far-axis-norm-x (if x-axis-y-min? 0.5 -0.5))
(define z-far-axis-norm-y (if y-axis-x-min? -0.5 0.5))
(define near-dist^2 (sqr (* 3 (plot-line-width))))
(define (vnear? v1 v2)
((vmag^2 (v- (plot->dc v1) (plot->dc v2))) . <= . near-dist^2))
@ -432,13 +454,13 @@
0))
(define (get-x-label-params)
(define v0 (plot->dc/no-axis-trans (vector x-mid x-axis-y z-min)))
(define v0 (norm->dc (vector 0.0 x-axis-norm-y -0.5)))
(define dist (+ max-x-tick-offset (max-x-tick-label-diag) half-char-height))
(list (plot-x-label) (v+ v0 (v* (y-axis-dir) (if x-axis-y-min? (- dist) dist)))
'top (- (if x-axis-y-min? 0 pi) (x-axis-angle))))
(define (get-y-label-params)
(define v0 (plot->dc/no-axis-trans (vector y-axis-x y-mid z-min)))
(define v0 (norm->dc (vector y-axis-norm-x 0.0 -0.5)))
(define dist (+ max-y-tick-offset (max-y-tick-label-diag) half-char-height))
(list (plot-y-label) (v+ v0 (v* (x-axis-dir) (if y-axis-x-min? (- dist) dist)))
'top (- (if y-axis-x-min? pi 0) (y-axis-angle))))
@ -449,13 +471,13 @@
'bottom-left 0))
(define (get-x-far-label-params)
(define v0 (plot->dc/no-axis-trans (vector x-mid x-far-axis-y z-min)))
(define v0 (norm->dc (vector 0.0 x-far-axis-norm-y -0.5)))
(define dist (+ max-x-far-tick-offset (max-x-far-tick-label-diag) half-char-height))
(list (plot-x-far-label) (v+ v0 (v* (y-axis-dir) (if x-axis-y-min? dist (- dist))))
'bottom (- (if x-axis-y-min? 0 pi) (x-axis-angle))))
(define (get-y-far-label-params)
(define v0 (plot->dc/no-axis-trans (vector y-far-axis-x y-mid z-min)))
(define v0 (norm->dc (vector y-far-axis-norm-x 0.0 -0.5)))
(define dist (+ max-y-far-tick-offset (max-y-far-tick-label-diag) half-char-height))
(list (plot-y-far-label) (v+ v0 (v* (x-axis-dir) (if y-axis-x-min? dist (- dist))))
'bottom (- (if y-axis-x-min? pi 0) (y-axis-angle))))
@ -519,6 +541,7 @@
;(printf "label params = ~v~n" (get-all-label-params))
;(printf "tick params = ~v~n" (get-all-tick-params))
(set! view->dc (make-view->dc left right top bottom))
;(printf "~v~n" (get-all-tick-params))
(append (append* (map (λ (params) (send/apply pd get-text-corners params))
(get-all-label-params)))
(append* (map (λ (params) (send/apply pd get-tick-endpoints (rest params)))
@ -545,32 +568,32 @@
(send pd set-minor-pen)
(when (plot-x-axis?)
(send pd draw-line
(plot->dc/no-axis-trans (vector x-min x-axis-y z-min))
(plot->dc/no-axis-trans (vector x-max x-axis-y z-min))))
(norm->dc (vector -0.5 x-axis-norm-y -0.5))
(norm->dc (vector 0.5 x-axis-norm-y -0.5))))
(when (plot-x-far-axis?)
(send pd draw-line
(plot->dc/no-axis-trans (vector x-min x-far-axis-y z-min))
(plot->dc/no-axis-trans (vector x-max x-far-axis-y z-min))))
(norm->dc (vector -0.5 x-far-axis-norm-y -0.5))
(norm->dc (vector 0.5 x-far-axis-norm-y -0.5))))
(when (plot-y-axis?)
(send pd draw-line
(plot->dc/no-axis-trans (vector y-axis-x y-min z-min))
(plot->dc/no-axis-trans (vector y-axis-x y-max z-min))))
(norm->dc (vector y-axis-norm-x -0.5 -0.5))
(norm->dc (vector y-axis-norm-x 0.5 -0.5))))
(when (plot-y-far-axis?)
(send pd draw-line
(plot->dc/no-axis-trans (vector y-far-axis-x y-min z-min))
(plot->dc/no-axis-trans (vector y-far-axis-x y-max z-min))))))
(norm->dc (vector y-far-axis-norm-x -0.5 -0.5))
(norm->dc (vector y-far-axis-norm-x 0.5 -0.5))))))
(define (draw-front-axes)
(when (plot-decorations?)
(send pd set-minor-pen)
(when (plot-z-axis?)
(send pd draw-line
(plot->dc/no-axis-trans (vector z-axis-x z-axis-y z-min))
(plot->dc/no-axis-trans (vector z-axis-x z-axis-y z-max))))
(norm->dc (vector z-axis-norm-x z-axis-norm-y -0.5))
(norm->dc (vector z-axis-norm-x z-axis-norm-y 0.5))))
(when (plot-z-far-axis?)
(send pd draw-line
(plot->dc/no-axis-trans (vector z-far-axis-x z-far-axis-y z-min))
(plot->dc/no-axis-trans (vector z-far-axis-x z-far-axis-y z-max))))))
(norm->dc (vector z-far-axis-norm-x z-far-axis-norm-y -0.5))
(norm->dc (vector z-far-axis-norm-x z-far-axis-norm-y 0.5))))))
(define (draw-ticks tick-params)
(for ([params (in-list tick-params)])
@ -590,18 +613,19 @@
(define (add-shapes! shapes) (set! render-list (append shapes render-list)))
(define (draw-shapes ss)
(define s+cs (map (λ (s) (cons s (plot->view/no-rho (shape-center s)))) ss))
(define s+cs (map (λ (s) (cons s (norm->view/no-rho (shape-center s)))) ss))
(for ([s+c (in-list (depth-sort (reverse s+cs)))])
(match-define (cons s c) s+c)
(draw-shape s (rotate/rho c))))
(define (draw-polygon alpha center vs norm pen-color pen-width pen-style brush-color brush-style)
(define-values (diff spec) (get-light-values center (norm->view norm)))
(define (draw-polygon alpha center vs normal
pen-color pen-width pen-style brush-color brush-style)
(define-values (diff spec) (get-light-values center (norm->view normal)))
(let ([pen-color (map (λ (v) (+ (* v diff) spec)) pen-color)]
[brush-color (map (λ (v) (+ (* v diff) spec)) brush-color)])
(send pd set-pen pen-color pen-width pen-style)
(send pd set-brush brush-color brush-style)
(send pd draw-polygon (map (λ (v) (plot->dc v)) vs))))
(send pd draw-polygon (map (λ (v) (norm->dc v)) vs))))
(define (draw-shape s center)
(send pd set-alpha (shape-alpha s))
@ -609,19 +633,20 @@
;; shapes
[(shapes alpha _ ss) (draw-shapes ss)]
;; polygon
[(polygon alpha _ vs norm pen-color pen-width pen-style brush-color brush-style)
(draw-polygon alpha center vs norm pen-color pen-width pen-style brush-color brush-style)]
[(polygon alpha _ vs normal pen-color pen-width pen-style brush-color brush-style)
(draw-polygon alpha center vs normal pen-color pen-width pen-style brush-color brush-style)]
;; rectangle
[(rectangle alpha _ r pen-color pen-width pen-style brush-color brush-style)
(for ([face (in-list (rect-visible-faces r theta))])
(match face
[(list norm vs ...) (draw-polygon alpha center vs norm pen-color pen-width pen-style
brush-color brush-style)]
[(list normal vs ...) (draw-polygon alpha center vs
normal pen-color pen-width pen-style
brush-color brush-style)]
[_ (void)]))]
;; line
[(line alpha _ v1 v2 pen-color pen-width pen-style)
(send pd set-pen pen-color pen-width pen-style)
(send pd draw-line (plot->dc v1) (plot->dc v2))]
(send pd draw-line (norm->dc v1) (norm->dc v2))]
;; text
[(text alpha _ anchor angle str font-size font-family color)
(send pd set-font font-size font-family)
@ -639,21 +664,14 @@
;; arrow glyph
[(arrow-glyph alpha _ v1 v2 pen-color pen-width pen-style)
(send pd set-pen pen-color pen-width pen-style)
(send pd draw-arrow (plot->dc v1) (plot->dc v2))]
(send pd draw-arrow (norm->dc v1) (norm->dc v2))]
[_ (error 'draw-shapes "shape not implemented: ~e" s)]))
;; Use a special view transform for the light so that the light angle is always the same
;; regardless of theta (but rotates rho). This also doesn't do any axis transforms, which could
;; fail; e.g. log transform when the light is at a negative position.
(define light-transform-matrix (m3* (m3-rotate-x rho) scale-matrix))
(define (plot->light-view v) (m3-apply light-transform-matrix (center v)))
;; Light position, in normalized view coordinates: 5 units up, ~3 units back and to the left
;; (simulates non-noon daylight conditions)
(define light (vector-map exact->inexact
(plot->light-view (vector (- x-min (* 2 x-size))
(- y-min (* 2 y-size))
(+ z-max (* 5 z-size))))))
(define light (m3-apply rotate-rho-matrix (vector (- -0.5 2.0)
(- -0.5 2.0)
(+ 0.5 5.0))))
;; View direction, in normalized view coordinates: many graph widths backward
(define view-dir (vector 0.0 -50.0 0.0))
@ -663,17 +681,17 @@
(define get-light-values
(cond
[(not (or diffuse-light? specular-light?)) (λ (v norm) (values 1.0 0.0))]
[(not (or diffuse-light? specular-light?)) (λ (v normal) (values 1.0 0.0))]
[else
(λ (v norm)
(λ (v normal)
; common lighting values
(define light-dir (vnormalize (v- light v)))
; diffuse lighting: typical Lambertian surface model
(define diff (if diffuse-light? (abs (vdot norm light-dir)) 1.0))
(define diff (if diffuse-light? (abs (vdot normal light-dir)) 1.0))
; specular highlighting: Blinn-Phong model
(define spec (cond [specular-light?
(define lv (v* (v+ light-dir view-dir) 0.5))
(define cos-angle (/ (abs (vdot norm lv)) (vmag lv)))
(define cos-angle (/ (abs (vdot normal lv)) (vmag lv)))
(* 32.0 (expt cos-angle 10.0))]
[else 0.0]))
; put it all together
@ -793,12 +811,13 @@
(values v1 v2))])
(unless (and v1 v2) (return (void)))
(cond [identity-transforms?
(add-shape! (line alpha c v1 v2
(add-shape! (line alpha (plot->norm c) (plot->norm v1) (plot->norm v2)
pen-color pen-width pen-style))]
[else
(define vs (subdivide-line plot->dc v1 v2))
(for ([v1 (in-list vs)] [v2 (in-list (rest vs))])
(add-shape! (line alpha c v1 v2 pen-color pen-width pen-style)))]))))
(add-shape! (line alpha (plot->norm c) (plot->norm v1) (plot->norm v2)
pen-color pen-width pen-style)))]))))
(define/public (put-lines vs)
(for ([vs (vregular-sublists vs)])
@ -811,7 +830,7 @@
(when (or (empty? vs) (not (and (andmap vregular? vs) (vregular? c))))
(return lst))
(define norm (vnormal (map plot->norm vs)))
(define normal (vnormal (map plot->norm vs)))
(let* ([vs (if clipping?
(clip-polygon vs clip-x-min clip-x-max
clip-y-min clip-y-max
@ -819,7 +838,8 @@
vs)]
[vs (if identity-transforms? vs (subdivide-polygon plot->dc vs))])
(when (empty? vs) (return lst))
(cons (polygon alpha c vs norm pen-color pen-width pen-style brush-color brush-style)
(cons (polygon alpha (plot->norm c) (map plot->norm vs)
normal pen-color pen-width pen-style brush-color brush-style)
lst))))
(define/public (put-polygon vs [c (vcenter vs)])
@ -830,33 +850,43 @@
#:when (not (empty? vs)))
(add-polygon lst vs (vcenter vs))))
(when (not (empty? lst))
(add-shape! (shapes alpha c lst))))
(add-shape! (shapes alpha (plot->norm c) lst))))
(define/public (put-rect r [c (rect-center r)])
(when (rect-regular? r)
(let ([r (rect-meet r bounds-rect)])
(add-shape! (rectangle alpha c r pen-color pen-width pen-style brush-color brush-style)))))
(match-define (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max)) r)
(match-let ([(vector x-min y-min z-min) (plot->norm (vector x-min y-min z-min))]
[(vector x-max y-max z-max) (plot->norm (vector x-max y-max z-max))])
(add-shape! (rectangle alpha (plot->norm c)
(vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max))
pen-color pen-width pen-style brush-color brush-style))))))
(define/public (put-text str v [anchor 'center] [angle 0])
(when (and (vregular? v) (in-bounds? v))
(add-shape! (text alpha v anchor angle str font-size font-family text-foreground))))
(add-shape! (text alpha (plot->norm v) anchor angle str
font-size font-family text-foreground))))
(define/public (put-glyphs vs symbol size)
(for ([v (in-list vs)])
(when (and (vregular? v) (in-bounds? v))
(add-shape!
(glyph alpha v symbol size pen-color pen-width pen-style brush-color brush-style)))))
(glyph alpha (plot->norm v) symbol size
pen-color pen-width pen-style brush-color brush-style)))))
(define/public (put-arrow v1 v2 [c (v* (v+ v1 v2) 1/2)])
(when (and (vregular? v1) (vregular? v2) (in-bounds? v1))
(cond [(in-bounds? v2)
(add-shape!
(arrow-glyph alpha c v1 v2 (->brush-color (plot-background)) (+ 2 pen-width) 'solid))
(arrow-glyph alpha (plot->norm c) (plot->norm v1) (plot->norm v2)
(->brush-color (plot-background)) (+ 2 pen-width) 'solid))
(add-shape!
(arrow-glyph alpha c v1 v2 pen-color pen-width pen-style))]
(arrow-glyph alpha (plot->norm c) (plot->norm v1) (plot->norm v2)
pen-color pen-width pen-style))]
[else (put-line v1 v2)])))
(define/public (put-tick v radius angle)
(when (and (vregular? v) (in-bounds? v))
(add-shape! (tick-glyph alpha v radius angle pen-color pen-width pen-style))))
(add-shape! (tick-glyph alpha (plot->norm v) radius angle
pen-color pen-width pen-style))))
)) ; end class

5
collects/plot/plot3d/plot.rkt
View File

@ -36,14 +36,12 @@
(define (get-bounds-rect renderer-list x-min x-max y-min y-max z-min z-max)
(define given-bounds-rect (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max)))
(define plot-bounds-rect (bounds-fixpoint renderer-list given-bounds-rect))
(when (or (not (rect-regular? plot-bounds-rect))
(rect-zero-area? plot-bounds-rect))
(match-define (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max)) plot-bounds-rect)
(error 'plot "could not determine sensible plot bounds; got x ∈ [~a,~a], y ∈ [~a,~a], z ∈ [~a,~a]"
x-min x-max y-min y-max z-min z-max))
(rect-inexact->exact plot-bounds-rect))
plot-bounds-rect)
(define (get-ticks renderer-list bounds-rect)
(define-values (all-x-ticks all-x-far-ticks all-y-ticks all-y-far-ticks all-z-ticks all-z-far-ticks)
@ -56,7 +54,6 @@
(define ticks-fun (plot-element-ticks-fun r))
(cond [ticks-fun (ticks-fun bounds-rect)]
[else (values empty empty empty empty empty empty)])))
(values (remove-duplicates (append* all-x-ticks))
(remove-duplicates (append* all-x-far-ticks))
(remove-duplicates (append* all-y-ticks))

13
collects/plot/plot3d/surface.rkt
View File

@ -11,17 +11,20 @@
(define ((surface3d-render-proc f samples color style line-color line-width line-style alpha label)
area)
(match-define (vector (ivl x-min x-max) (ivl y-min y-max) z-ivl) (send area get-bounds-rect))
(match-define (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max))
(send area get-bounds-rect))
(define sample (f x-min x-max (animated-samples samples)
y-min y-max (animated-samples samples)))
(send area put-alpha alpha)
(send area put-brush color style)
(send area put-pen line-color line-width line-style)
(for-2d-sample
(xa xb ya yb z1 z2 z3 z4) sample
(send area put-polygon
(list (vector xa ya z1) (vector xb ya z2) (vector xb yb z3) (vector xa yb z4))))
(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
(send area put-polygon
(list (vector xa ya z1) (vector xb ya z2) (vector xb yb z3) (vector xa yb z4)))))
(cond [label (rectangle-legend-entry label color style line-color line-width line-style)]
[else empty]))

56
collects/plot/tests/extreme-bounds-tests.rkt Normal file
View File

@ -0,0 +1,56 @@
#lang racket
(require plot plot/utils)
(plot-x-label #f)
(plot-y-label #f)
(define stops (list (* 2 (inexact->exact -max.0))
(inexact->exact -max.0)
(inexact->exact -min.0)
0
(inexact->exact +min.0)
(inexact->exact +max.0)
(* 2 (inexact->exact +max.0))))
(define (extreme-real->string x)
(real->plot-label x (digits-for-range 0 (abs x))))
(define num-2d 0)
(time
(for ([x-min (in-list stops)]
[x-max (in-list (rest stops))]
#:when #t
[y-min (in-list stops)]
[y-max (in-list (rest stops))])
(displayln
(list (format "[~a,~a] × [~a,~a]"
(extreme-real->string x-min) (extreme-real->string x-max)
(extreme-real->string y-min) (extreme-real->string y-max))
(time (plot (lines (list (vector x-min y-min) (vector x-max y-max)))))))
(set! num-2d (+ 1 num-2d))))
(printf "Total 2D plots: ~a~n" num-2d)
(printf "-----------------------------------------~n")
(define num-3d 0)
(time
(for ([x-min (in-list stops)]
[x-max (in-list (rest stops))]
#:when #t
[y-min (in-list stops)]
[y-max (in-list (rest stops))]
#:when #t
[z-min (in-list stops)]
[z-max (in-list (rest stops))])
(displayln
(list (format "[~a,~a] × [~a,~a] × [~a,~a]"
(extreme-real->string x-min) (extreme-real->string x-max)
(extreme-real->string y-min) (extreme-real->string y-max)
(extreme-real->string z-min) (extreme-real->string z-max))
(time (plot3d (lines3d (list (vector x-min y-min z-min) (vector x-max y-max z-max)))))))
(set! num-3d (+ 1 num-3d))))
(printf "Total plots: ~a~n" num-3d)
(printf "-----------------------------------------~n")

6
collects/plot/tests/isosurface-tests.rkt
View File

@ -55,3 +55,9 @@
#:color "navajowhite" #:width 2)
(lines3d '(#(0 0 2) #(0 0 -2)) #:color "navajowhite" #:width 2))
#:title "A Seashell" #:x-label #f #:y-label #f #:angle 210 #:altitude 30))
(define flepsilon (flnext 0.0))
(time
(plot3d (isosurface3d (λ (x y z) (+ (sqr x) (sqr y) (sqr z))) (sqr (inexact->exact flepsilon))
(- flepsilon) flepsilon (- flepsilon) flepsilon (- flepsilon) flepsilon)))

4
collects/plot/tests/plot2d-tests.rkt
View File

@ -18,6 +18,10 @@
(time (plot (function values 0 1000)))
(time
(define flepsilon (ordinal->flonum 1))
(plot (lines (list (vector 0 0) (vector flepsilon flepsilon)))))
(parameterize ([plot-x-transform log-transform]
[plot-x-ticks (log-ticks #:base 4)]
[plot-y-transform log-transform]

9
collects/plot/tests/plot3d-tests.rkt
View File

@ -10,6 +10,15 @@
(time
(plot3d (points3d empty) #:x-min -1 #:x-max 1 #:y-min -1 #:y-max 1 #:z-min -1 #:z-max 1))
(define flepsilon (ordinal->flonum 1))
(time
(plot3d (lines3d (list (vector 0 0 0) (vector flepsilon flepsilon flepsilon)))))
(time
(plot3d (contour-intervals3d (λ (x y) (- (sqr x) (sqr y)))
(- flepsilon) flepsilon (- flepsilon) flepsilon)))
(parameterize ([plot-background "black"]
[plot-foreground "white"]
[plot-background-alpha 1/2]