From 56f70fb4f2b3af001b6b41dc7763b51608a693cb Mon Sep 17 00:00:00 2001
From: Neil Toronto <ntoronto@racket-lang.org>
Date: Fri, 4 Nov 2011 23:53:22 -0600
Subject: [PATCH] Exposed marching squares and cubes algorithms through
 plot/utils, speed improvements, fixes

---
 collects/plot/common/marching-cubes.rkt   | 170 ++++++++++------------
 collects/plot/common/marching-squares.rkt |  48 ++----
 collects/plot/common/marching-utils.rkt   |  32 ++++
 collects/plot/common/math.rkt             |  28 ++--
 collects/plot/doc.rkt                     |   8 +-
 collects/plot/plot3d/isosurface.rkt       |  45 ++----
 collects/plot/plot3d/plot-area.rkt        |   2 +-
 collects/plot/scribblings/todo.scrbl      |   3 +
 collects/plot/tests/isosurface-tests.rkt  |   8 +-
 collects/plot/utils.rkt                   |   2 +
 10 files changed, 164 insertions(+), 182 deletions(-)
 create mode 100644 collects/plot/common/marching-utils.rkt

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