This section implements the find-contours2d function, which accepts
a 2D ndarray and a scalar. It computes a list of point sequences
representing all interpolated level crossings for the given value.
These contour lines can then be further processed or visualized, as
shown in the example below…
The following demo utilizes thi.ng/geom, thi.ng/math and thi.ng/color libs to create SVG based contour maps of generated dummy data. To try it out yourself, simply add the following dependency (in addition to this (ndarray) library) - geom itself depends on math & color so they don’t need to be specified in your project:
[thi.ng/geom "0.0.815"]These images were generated with the demo code below and show the impact of different matrix resolutions on the precision and quality of the resulting visualizations.
Note: Click on images to view in bigger resolution
 |  |  |
| 32 x 32 | 64 x 64 | 128 x 128 |
| 8916 vertices | 19172 vertices | 39702 vertices |
 |  |  |
| 32 x 32 | 64 x 64 | 128 x 128 |
| 7020 vertices | 14652 vertices | 29874 vertices |
(require
'[thi.ng.ndarray.core :as nd]
'[thi.ng.ndarray.contours :as contours]
'[thi.ng.geom.core :as g]
'[thi.ng.geom.core.vector :as v]
'[thi.ng.geom.svg.core :as svg]
'[thi.ng.math.core :as m]
'[thi.ng.math.simplexnoise :as n]
'[thi.ng.color.gradients :as grad])
(def res 128)
(def width 640.0)
(def scale (/ width (- res 2)))
(def clipped (- width (* 2.0 scale)))
(def n-scale 0.03)
(def num-contours 60)
(defn contour->svg
"Takes a single seq of contour coordinates and converts it into an
SVG polygon (hiccup format)."
[contour]
(-> (map #(-> % v/vec2 v/yx (g/scale scale)) contour)
(svg/polygon)))
(defn noise-matrix
"Creates a new 2D matrix of size res and populates it with simplex
noise, then sets border cells to 1.0 and returns matrix"
[res ns]
(let [mat (nd/ndarray :float32 (float-array (* res res)) [res res])]
(dorun
(for [[y x] (nd/position-seq mat)]
(nd/set-at mat x y (+ 0.5 (* 0.5 (n/noise2 (+ 101 (* x ns)) (* y ns)))))))
(contours/set-border-2d mat 1)))
(defn circle-matrix
"Creates new 2D matrix of size res and populates it w/ modulated
normalized distance values from center. `spikes` arg is number of
oscillations used to modulate. `amp` is modulation strength. Sets
border cells to 1.0 and returns matrix."
[res spikes amp]
(let [mat (nd/ndarray :float32 (float-array (* res res)) [res res])
c (/ res 2.0)
dmax (* m/SQRT2 0.5 res)]
(dorun
(for [[y x] (nd/position-seq mat)
:let [dx (- c x)
dy (- c y)
t (Math/atan2 dy dx)
d (Math/sqrt (* (+ 0.5 (* amp (Math/sin (* t spikes))))
(+ (* dx dx) (* dy dy))))]]
(nd/set-at mat x y (/ d dmax))))
(contours/set-border-2d mat 1)))
(def palette
(apply grad/cosine-gradient num-contours (:orange-blue grad/cosine-schemes)))
;; choose one...
(def mat (noise-matrix res n-scale))
;;(def mat (circle-matrix res 8 0.25))
(->> (m/norm-range num-contours)
(rest)
(map
#(svg/group
{:stroke (palette (int (* % (dec num-contours))))
:fill "none"}
(map contour->svg (contours/find-contours-2d mat %))))
(svg/svg
{:width width
:height width
:viewBox (format "%1.2f %1.2f %1.2f %1.2f" scale scale clipped clipped)})
(svg/serialize)
(spit "iso.svg"))Loosely based on Marching Squares/Cubes implementations by Paul Bourke (C) & Murphy Stein (Java):
- http://paulbourke.net/geometry/polygonise/
- https://github.com/murphydactyl/JavaKinectFingerTracker/blob/master/imageprocessor/FindIsolines.java
(def edge-index-2d
[nil [2 0] [1 0] [1 0]
[0 0] nil [0 0] [0 0]
[3 0] [2 0] nil [1 0]
[3 0] [2 0] [3 0] nil])
(def next-edges-2d
[[-1 0] [0 1] [1 0] [0 -1]])
(defn set-border-2d
[mat x]
(let [[h w] (nd/shape mat)
h' (dec h)
w' (dec w)
l (nd/pick mat nil 0)
r (nd/pick mat nil w')
t (nd/pick mat 0 nil)
b (nd/pick mat h' nil)]
(loop [i w']
(when (>= i 0)
(nd/set-at t i x)
(nd/set-at b i x)
(recur (dec i))))
(loop [i h']
(when (>= i 0)
(nd/set-at l i x)
(nd/set-at r i x)
(recur (dec i))))
mat))
(defn encode-crossings-2d
[src isoval]
(let [out (nd/ndarray :int8 (#?(:clj byte-array :cljs a/int8) (nd/size src)) (nd/shape src))
iso? (fn [y x m] (if (< (nd/get-at src y x) isoval) m 0))]
(loop [pos (nd/position-seq (nd/truncate-h src -1 -1))]
(if pos
(let [[y x] (first pos)
x' (inc x)
y' (inc y)]
(nd/set-at
out y x
(-> (iso? y x 0x08)
(bit-or (iso? y x' 0x04))
(bit-or (iso? y' x' 0x02))
(bit-or (iso? y' x 0x01))))
(recur (next pos)))
out))))
(defn mean-cell-value-2d
[src y x]
(* (+ (+ (nd/get-at src y x) (nd/get-at src y (inc x)))
(+ (nd/get-at src (inc y) x) (nd/get-at src (inc y) (inc x))))
0.25))
(defn process-saddle5
[src y x iso from]
(if (> (mean-cell-value-2d src y x) iso)
(if (== 3 from) [2 0x04] [0 0x01])
(if (== 3 from) [0 0x0d] [2 0x07])))
(defn process-saddle10
[src y x iso from]
(if (> (mean-cell-value-2d src y x) iso)
(if (== 0 from) [3 0x02] [1 0x08])
(if (== 2 from) [3 0x0b] [1 0x0e])))
(defn mix2d
[src y1 x1 y2 x2 iso]
(let [a (nd/get-at src y1 x1)
b (nd/get-at src y2 x2)]
(if (== a b) 0 (/ (- a iso) (- a b)))))
(defn contour-vertex-2d
[src y x to iso]
(let [x' (inc x) y' (inc y)]
(case (int to)
0 [y (+ x (mix2d src y x y x' iso))]
1 [(+ y (mix2d src y x' y' x' iso)) x']
2 [y' (+ x (mix2d src y' x y' x' iso))]
3 [(+ y (mix2d src y x y' x iso)) x]
nil)))
(defn find-contours-2d
[src isoval]
(let [[h' w'] (nd/shape src)
h' (dec h')
w' (dec w')
coded (encode-crossings-2d src isoval)
contours (volatile! (transient []))]
(loop [pos (nd/position-seq coded)
curr (transient [])
to nil
p nil]
(if pos
(let [from to
[y x] (if p p (first pos))]
(if (or (>= x w') (>= y h'))
(recur (next pos) curr to nil)
(let [id (nd/get-at coded y x)
[to clear] (case (int id)
5 (process-saddle5 src y x isoval from)
10 (process-saddle10 src y x isoval from)
(edge-index-2d (int id)))
curr (if (and (nil? from) to (pos? (count curr)))
(do (vswap! contours conj! (persistent! curr))
(transient []))
curr)]
(when clear
(nd/set-at coded y x clear))
(if (and to (>= to 0))
(let [vertex (contour-vertex-2d src y x to isoval)
[oy ox] (next-edges-2d to)]
(recur (next pos) (conj! curr vertex) to [(+ y oy) (+ x ox)]))
(recur (next pos) curr to nil)))))
(persistent! (conj! @contours (persistent! curr)))))))This functionality will be ported from the thi.ng/geom-voxel module…
This section contains somewhat less high level operations to find
level crossings in 1D, 2D and 3D ndarrays. Unlike the contour
extractions above, which procude a logical sequence of connected
points/segments/facets in the grid, these functions here merely
produce a sequence of (potentially) unconnected cell positions where
thresholds are crossed and are intended for more analytical use cases.
The functions all take an ndarray and contour level value and assume
the array to be in major shape order (the default order), i.e. in 2D
row-major (YX), in 3D slice-row-major (ZYX). The functions return
interpolated grid positions where the given contour level is crossed
between cells. If the array doesn’t conform to this ordering, use the
nd/transpose method to create a properly ordered view before using
the functions below.
(level-crossings1d (nd/ndarray :float32 [0 0 1 0]) 4 0.25)
;; (1.25 2.75)
(let [a [0 0 0
0 1 0
0 0 0]
a (nd/ndarray :float32 a [3 3])]
{:x (level-crossings2d-x a 0.25)
:y (level-crossings2d-y a 0.25)
:all (level-crossings2d a 0.25)})
;; {:x ([1 0.25] [1 1.75])
;; :y ([0.25 1] [1.75 1])
;; :all ([1 0.25] [1 1.75] [0.25 1] [1.75 1])}
(let [a (nd/ndarray :float32 (float-array 27) [3 3 3])]
(nd/set-at a 1 1 1 1)
{:x (level-crossings3d-x a 0.25)
:y (level-crossings3d-y a 0.25)
:z (level-crossings3d-z a 0.25)
:all (level-crossings3d a 0.25)})
;; {:x ((1 1 0.25) (1 1 1.75))
;; :y ((1 0.25 1) (1 1.75 1))
;; :z ([0.25 1 1] [1.75 1 1])
;; :all ((1 1 0.25) (1 1 1.75) (1 0.25 1) (1 1.75 1) [0.25 1 1] [1.75 1 1])}(defn level-crossing
[offset a b level]
(let [da (- a level)
db (- b level)]
(if-not (= (>= da 0.0) (>= db 0.0))
(+ offset (+ 0.5 (* 0.5 (/ (+ da db) (- da db))))))))
(defn level-crossings1d
[mat shape level]
(for [x (range (dec (if (number? shape) shape (first shape))))
:let [x' (level-crossing x (nd/get-at mat x) (nd/get-at mat (inc x)) level)]
:when x']
x'))
(defn level-crossings2d-x
([mat level]
(level-crossings2d-x mat (nd/shape mat) level))
([mat [sy sx] level]
(mapcat
(fn [y] (map #(vector y %) (level-crossings1d (nd/pick mat y nil) sx level)))
(range sy))))
(defn level-crossings2d-y
([mat level]
(level-crossings2d-y mat (nd/shape mat) level))
([mat [sy sx] level]
(mapcat
(fn [x] (map #(vector % x) (level-crossings1d (nd/pick mat nil x) sy level)))
(range sx))))
(defn level-crossings2d
([mat level]
(level-crossings2d mat (nd/shape mat) level))
([mat shape level]
(concat
(level-crossings2d-x mat shape level)
(level-crossings2d-y mat shape level))))
(defn level-crossings3d-x
([mat level]
(level-crossings3d-x mat (nd/shape mat) level))
([mat [sz sy sx] level]
(mapcat
(fn [z] (map #(cons z %) (level-crossings2d-x (nd/pick mat z nil nil) [sy sx] level)))
(range sz))))
(defn level-crossings3d-y
([mat level]
(level-crossings3d-y mat (nd/shape mat) level))
([mat [sz sy sx] level]
(mapcat
(fn [z] (map #(cons z %) (level-crossings2d-y (nd/pick mat z nil nil) [sy sx] level)))
(range sz))))
(defn level-crossings3d-z
([mat level]
(level-crossings3d-z mat (nd/shape mat) level))
([mat [sz sy sx] level]
(mapcat
(fn [x] (map #(conj % x) (level-crossings2d-y (nd/pick mat nil nil x) [sz sy] level)))
(range sx))))
(defn level-crossings3d
([mat level]
(level-crossings3d mat (nd/shape mat) level))
([mat shape level]
(concat
(level-crossings3d-x mat shape level)
(level-crossings3d-y mat shape level)
(level-crossings3d-z mat shape level))))(ns thi.ng.ndarray.contours
(:require
[thi.ng.ndarray.core :as nd]
#?(:cljs [thi.ng.typedarrays.core :as a])))
<<level-crossings>>
<<contours2d>>