racket/collects/plot/plot3d/isosurface.rkt

276 lines
12 KiB
Racket
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#lang racket/base
(require racket/class racket/match racket/list racket/flonum racket/contract racket/math
plot/utils
"../common/contract-doc.rkt")
(provide (all-defined-out))
;; ===================================================================================================
;; Surfaces of constant value (isosurfaces)
(define ((isosurface3d-render-proc
f d 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) (ivl z-min z-max))
(send area get-bounds-rect))
(match-define (3d-sample xs ys zs dsss d-min d-max)
(f x-min x-max (animated-samples samples)
y-min y-max (animated-samples samples)
z-min z-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 ([za (in-list zs)]
[zb (in-list (rest zs))]
[dss0 (in-vector dsss)]
[dss1 (in-vector dsss 1)]
#:when #t
[ya (in-list ys)]
[yb (in-list (rest ys))]
[ds00 (in-vector dss0)]
[ds01 (in-vector dss0 1)]
[ds10 (in-vector dss1)]
[ds11 (in-vector dss1 1)]
#:when #t
[xa (in-list xs)]
[xb (in-list (rest xs))]
[d1 (in-vector ds00)]
[d2 (in-vector ds00 1)]
[d3 (in-vector ds01 1)]
[d4 (in-vector ds01)]
[d5 (in-vector ds10)]
[d6 (in-vector ds10 1)]
[d7 (in-vector ds11 1)]
[d8 (in-vector ds11)])
(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)]
[else empty]))
(defproc (isosurface3d [f (real? real? real? . -> . real?)] [d real?]
[x-min (or/c regular-real? #f) #f]
[x-max (or/c regular-real? #f) #f]
[y-min (or/c regular-real? #f) #f]
[y-max (or/c regular-real? #f) #f]
[z-min (or/c regular-real? #f) #f]
[z-max (or/c regular-real? #f) #f]
[#:samples samples (and/c exact-integer? (>=/c 2)) (plot3d-samples)]
[#:color color plot-color/c (surface-color)]
[#:style style plot-brush-style/c (surface-style)]
[#:line-color line-color plot-color/c (surface-line-color)]
[#:line-width line-width (>=/c 0) (surface-line-width)]
[#:line-style line-style plot-pen-style/c (surface-line-style)]
[#:alpha alpha (real-in 0 1) (surface-alpha)]
[#:label label (or/c string? #f) #f]
) renderer3d?
(define g (3d-function->sampler f))
(renderer3d (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max)) #f default-ticks-fun
(isosurface3d-render-proc
g d samples color style line-color line-width line-style alpha label)))
;; ===================================================================================================
;; Nested isosurfaces
(define ((isosurfaces3d-render-proc f rd-min rd-max levels samples colors styles
line-colors line-widths line-styles alphas label)
area)
(match-define (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max))
(send area get-bounds-rect))
(match-define (3d-sample xs ys zs dsss fd-min fd-max)
(f x-min x-max (animated-samples samples)
y-min y-max (animated-samples samples)
z-min z-max (animated-samples samples)))
(define d-min (if rd-min rd-min fd-min))
(define d-max (if rd-max rd-max fd-max))
(cond
[(not (and d-min d-max)) empty]
[else
(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)]
[color (in-cycle colors)]
[style (in-cycle styles)]
[alpha (in-cycle alphas)]
[line-color (in-cycle line-colors)]
[line-width (in-cycle line-widths)]
[line-style (in-cycle line-styles)])
(send area put-alpha alpha)
(send area put-brush color style)
(send area put-pen line-color line-width line-style)
(for ([za (in-list zs)]
[zb (in-list (rest zs))]
[dss0 (in-vector dsss)]
[dss1 (in-vector dsss 1)]
#:when #t
[ya (in-list ys)]
[yb (in-list (rest ys))]
[ds00 (in-vector dss0)]
[ds01 (in-vector dss0 1)]
[ds10 (in-vector dss1)]
[ds11 (in-vector dss1 1)]
#:when #t
[xa (in-list xs)]
[xb (in-list (rest xs))]
[d1 (in-vector ds00)]
[d2 (in-vector ds00 1)]
[d3 (in-vector ds01 1)]
[d4 (in-vector ds01)]
[d5 (in-vector ds10)]
[d6 (in-vector ds10 1)]
[d7 (in-vector ds11 1)]
[d8 (in-vector ds11)])
(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-entries
label ds colors styles line-colors line-widths line-styles)]
[else empty]))]))
(defproc (isosurfaces3d
[f (real? real? real? . -> . real?)]
[x-min (or/c regular-real? #f) #f] [x-max (or/c regular-real? #f) #f]
[y-min (or/c regular-real? #f) #f] [y-max (or/c regular-real? #f) #f]
[z-min (or/c regular-real? #f) #f] [z-max (or/c regular-real? #f) #f]
[#:d-min d-min (or/c regular-real? #f) #f] [#:d-max d-max (or/c regular-real? #f) #f]
[#:samples samples (and/c exact-integer? (>=/c 2)) (plot3d-samples)]
[#:levels levels (or/c 'auto exact-positive-integer? (listof real?)) (isosurface-levels)]
[#:colors colors (plot-colors/c (listof real?)) (isosurface-colors)]
[#:styles styles (plot-brush-styles/c (listof real?)) (isosurface-styles)]
[#:line-colors line-colors (plot-colors/c (listof real?)) (isosurface-line-colors)]
[#:line-widths line-widths (pen-widths/c (listof real?)) (isosurface-line-widths)]
[#:line-styles line-styles (plot-pen-styles/c (listof real?)) (isosurface-line-styles)]
[#:alphas alphas (alphas/c (listof real?)) (isosurface-alphas)]
[#:label label (or/c string? #f) #f]
) renderer3d?
(define g (3d-function->sampler f))
(renderer3d (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max)) #f default-ticks-fun
(isosurfaces3d-render-proc g d-min d-max levels samples colors styles
line-colors line-widths line-styles alphas
label)))
;; ===================================================================================================
(define ((polar3d-render-proc f g 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) (ivl z-min z-max))
(send area get-bounds-rect))
(match-define (3d-sample xs ys zs dsss d-min d-max)
(g x-min x-max (animated-samples samples)
y-min y-max (animated-samples samples)
z-min z-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 ([za (in-list zs)]
[zb (in-list (rest zs))]
[dss0 (in-vector dsss)]
[dss1 (in-vector dsss 1)]
#:when #t
[ya (in-list ys)]
[yb (in-list (rest ys))]
[ds00 (in-vector dss0)]
[ds01 (in-vector dss0 1)]
[ds10 (in-vector dss1)]
[ds11 (in-vector dss1 1)]
#:when #t
[xa (in-list xs)]
[xb (in-list (rest xs))]
[d1 (in-vector ds00)]
[d2 (in-vector ds00 1)]
[d3 (in-vector ds01 1)]
[d4 (in-vector ds01)]
[d5 (in-vector ds10)]
[d6 (in-vector ds10 1)]
[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 xa xb ya yb za zb 0.0 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 [(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)]))
(cond [label (rectangle-legend-entry
label color style line-color line-width line-style)]
[else empty]))
(define 2pi (* 2 pi))
(define ((2d-polar->3d-function f) x y z)
(let ([x (exact->inexact x)]
[y (exact->inexact y)]
[z (exact->inexact z)])
(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 (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] [#:x-max x-max (or/c regular-real? #f) #f]
[#:y-min y-min (or/c regular-real? #f) #f] [#:y-max y-max (or/c regular-real? #f) #f]
[#:z-min z-min (or/c regular-real? #f) #f] [#:z-max z-max (or/c regular-real? #f) #f]
[#:samples samples (and/c exact-integer? (>=/c 2)) (plot3d-samples)]
[#:color color plot-color/c (surface-color)]
[#:style style plot-brush-style/c (surface-style)]
[#:line-color line-color plot-color/c (surface-line-color)]
[#:line-width line-width (>=/c 0) (surface-line-width)]
[#:line-style line-style plot-pen-style/c (surface-line-style)]
[#: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.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))
(cond [(empty? rvs) (renderer3d #f #f #f #f)]
[else
(match-define (list (vector rxs rys rzs) ...) rvs)
(let ([x-min (if x-min x-min (apply min* rxs))]
[x-max (if x-max x-max (apply max* rxs))]
[y-min (if y-min y-min (apply min* rys))]
[y-max (if y-max y-max (apply max* rys))]
[z-min (if z-min z-min (apply min* rzs))]
[z-max (if z-max z-max (apply max* rzs))])
(define new-f (2d-polar->3d-function f))
(define g (3d-function->sampler new-f))
(renderer3d (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max)) #f
default-ticks-fun
(polar3d-render-proc new-f g samples color style
line-color line-width line-style alpha label)))]))