Delaunay diagram sheep

core.clj

(ns sheep.core
  (:require [quil.core :as q])
  (:require [quil.middleware :as m])
  (:require [sheep.delaunay :as d]))

(def next-color!
  (let [current (atom 50)
        hue (atom 150)]
    (fn [] (let [b (swap! current (fn [_] (+ 70 (rand-int 30))))
                 hue (swap! hue (fn [_] (+ 150 (rand-int 40))))]
            (q/color hue 30 b)))))

(defn inner-p [n center radius]
  (let [start-x (- (:x center) radius)
        start-y (- (:y center) radius)
        width (- (+ (:x center) radius) start-x)
        height (- (+ (:y center) radius) start-y)]
    (vec (loop [points [] len n]
      (if (= len 0)
        points
        (let [p (d/->Point (+ (rand-int width) start-x) (+ (rand-int height) start-y))
              dist (d/distance center p)]
          (if (< dist radius)
            (recur (conj points p) (dec len))
            (recur points len))))))))

(defn circle-points [center radius length]
  (let [{center-x :x center-y :y} center
        radian-unit (/ q/TWO-PI length)]
    (for [i (range length)]
      (let [radian (* i radian-unit)]
      (d/->Point (+ (* (q/cos radian) radius) center-x)
                 (+ (* (q/sin radian) radius) center-y))))))

(defn draw-tri [tri color]
  (q/fill color)
  (q/stroke color)
  (let [{p1 :p1 p2 :p2 p3 :p3} tri
        {x1 :x y1 :y} p1, {x2 :x y2 :y} p2, {x3 :x y3 :y} p3]
    (q/triangle x1 y1 x2 y2 x3 y3)))

(defn draw-foots []
  (q/fill (q/color 80))
  (q/no-stroke)
  (let [width 60, length 130]
    (q/with-translation [320 520]
      (q/with-rotation [(* 30 (/ q/PI 360))]
        (q/rect 0 0 width length 15)
        (q/rect 80 0 width length 15)))
    (q/with-translation [550 530]
      (q/with-rotation [(- (* 30 (/ q/PI 360)))]
        (q/rect 0 -10 width length 15)
        (q/rect 80 0 width length 15)
        ))))

(defn draw-face []
  (q/color-mode :rgb 255)
  (q/no-stroke)
  (q/with-translation [300 400]
    ; 顔
    (q/fill (q/color 80))
    (q/ellipse 0 0 250 200)
    ; 目
    (q/fill (q/color 245))
    (q/ellipse 70 -30 20 20)
    ; 鼻
    (q/triangle -60 20 -36 20 -48 30)
    ; 口
    (q/stroke-weight 3)
    (q/stroke-cap :round)
    (q/stroke (q/color 245))
    (q/curve -80 20 -70 50 0 50 10 20))
  ;角
  (q/stroke-weight 3)
  (q/stroke (q/color 188 189 85))
  (q/with-translation [510 350]
    (q/with-rotation [(+ q/PI (* 30 (/ q/PI 180)))]
    (loop [radian 0 offset 0 weight 1]
      (when (>= radian (- q/TWO-PI))
        (let [x1 (* (q/cos radian) offset), y1 (* (q/sin radian) offset)
              x2 (* (q/cos radian) (+ offset weight)), y2 (* (q/sin radian) (+ offset weight))]
        (q/line x1 y1 x2 y2))
        (recur (- radian (/ q/TWO-PI 360)) (+ offset 0.2) (+ weight 0.2)))))))

(defn setup []
  (q/background (q/color 255))
  (q/smooth)
  (draw-foots)
  (q/color-mode :hsb 400 100 100)
  (let [c-points (circle-points (d/->Point 500 300) 300 15)
        points (vec (concat c-points (inner-p 15 (d/->Point 500 300) 270)))
        tris (d/delaunay-diagram points)]
    (doseq [tri tris] (draw-tri tri (next-color!)) ))
  (draw-face))

(defn key-pressed [s event]
  (when (= (:raw-key event) \s)
    (q/save "delauny.jpg"))
  s)

(defn start-sketch []
  (q/sketch
    :title "sheep"
    :size [1000 800]
    :setup setup
    :key-pressed key-pressed
    :middleware [m/fun-mode]))

(defn -main [& args]
  (start-sketch))

delaunay.clj

(ns sheep.delaunay)

(def sqrt3 "3の平方根" (Math/sqrt 3))
(defn abs [n] (max n (- n)))
(defn pow2 "数値を2乗する" [b] (Math/pow b 2))
(defn rm [coll item]
  (remove #(= item %) coll))

; -Shape --------------------------------------
(defprotocol Shape (contains-p? [this p]))

; -- Point ------------------------------------
(defrecord Point [x y])
(defn distance [p1 p2]
  "2点間の距離を求める"
  (let [{x1 :x y1 :y} p1, {x2 :x y2 :y} p2]
    (Math/sqrt (+ (pow2 (- x1 x2)) (pow2 (- y1 y2))))))

(defrecord Rectangle [point width height])

;-- Edge --------------------------------------
(defrecord Edge [p1 p2])
(defn side? [edge p]
  "点が線分のどちら側にあるか :left :right :online"
  (let [{p1 :p1 p2 :p2} edge
        {x1 :x y1 :y} p1, {x2 :x y2 :y} p2, {x :x y :y} p
        n (+ (* x (- y1 y2)) (* x1 (- y2 y)) (* x2 (- y y1)))]
    (cond
      (> n 0) :left
      (< n 0) :right
      :else   :online)))

;-- Triangle ----------------------------------
(defn get-points [tri] [(:p1 tri) (:p2 tri) (:p3 tri)])
(defn get-edges [tri] 
  (if tri
    (let [[p1 p2 p3] (get-points tri)]
      [(->Edge p1 p2) (->Edge p2 p3) (->Edge p3 p1)])
    []))

(defrecord Triangle [p1 p2 p3]
  Shape
  (contains-p? [this p]
    (let [[ab bc ca] (get-edges this)
          side-ab (side? ab p), side-bc (side? bc p), side-ca (side? ca p)]
      (or (some #(= :online %) [side-ab side-bc side-ca]) (= side-ab side-bc side-ca)))))

(defn new-triangle [[x1 y1] [x2 y2] [x3 y3]]
  (->Triangle (->Point x1 y1) (->Point x2 y2) (->Point x3 y3)))

(defn equal-tri? [tri1 tri2]
  (let [tri1-points (get-points tri1)
        {p1 :p1 p2 :p2 p3 :p3} tri2]
    (and (some #(= p1 %) tri1-points)
         (some #(= p2 %) tri1-points)
         (some #(= p3 %) tri1-points))))

(defn contains-edge? [tri edge]
  "三角形に辺が含まれるか"
  (let [tri-points (get-points tri)] 
    (and (some #(= (:p1 edge) %) tri-points)
         (some #(= (:p2 edge) %) tri-points))))

(defn get-another-p [tri edge]
  "辺の端点以外の頂点を求める"
  (let [points (get-points tri)
        equal-p1 #(= (:p1 edge) %)
        equal-p2 #(= (:p2 edge) %)]
    (some #(if (not (or (equal-p1 %) (equal-p2 %))) % nil) points)))


;-- Circle ------------------------------------
(defrecord Circle [center radius]
  Shape
  (contains-p? [this point]
    (let [{center :center radius :radius} this]
      (<= (distance center point) radius))))

        

;----------------------------------------------
(defn circumscribed-rect [points]
  "点のコレクションを全て内包する四角形を求める"
  (let [min-x (:x (apply min-key :x points))
        max-x (:x (apply max-key :x points))
        min-y (:y (apply min-key :y points))
        max-y (:y (apply max-key :y points))]
    (->Rectangle (->Point min-x min-y) (- max-x min-x) (- max-y min-y))))

(defn circumscribed-triangle [rect]
  "四角形に外接する三角形を作成する"
  (let [p (:point rect)
        width (:width rect)
        height (:height rect)
        center-x (+ (/ width 2) (:x p))
        center-y (+ (/ height 2) (:y p))
        center (->Point center-x center-y)
        radius (distance p center)
        p1 (->Point center-x (- center-y (* radius 2)))
        p2 (->Point (- center-x (* sqrt3 radius)) (+ center-y radius))
        p3 (->Point (+ center-x (* sqrt3 radius)) (+ center-y radius))]
        (->Triangle p1 p2 p3)))

(defn circumscribed-circle [tri]
  "三角形に外接する円を作成する"
  (let [{p1 :p1 p2 :p2 p3 :p3} tri
        {x1 :x y1 :y} p1, {x2 :x y2 :y} p2, {x3 :x y3 :y} p3
        ;中心点
        x1p2 (pow2 x1), y1p2 (pow2 y1), x2p2 (pow2 x2), y2p2 (pow2 y2), x3p2 (pow2 x3), y3p2 (pow2 y3)
        a (+ x2p2 (- x1p2) y2p2 (- y1p2))
        b (+ x3p2 (- x1p2) y3p2 (- y1p2))
        c (* 2 (- (* (- x2 x1) (- y3 y1)) (* (- y2 y1) (- x3 x1))))
        x (/ (+ (* (- y3 y1) a) (* (- y1 y2) b)) c)
        y (/ (+ (* (- x1 x3) a) (* (- x2 x1) b)) c)
        center (->Point x y)
        ;半径
        radius (distance p1 center)]
    (->Circle center radius)))

(defn flip [tri1 tri2 edge]
  "２個の三角形の共有する辺をフリップする
    引数：三角形1, 三角形2, 共有辺
    戻り値：新しい三角形1, 新しい三角形2, 辺の集合, 共有辺"
  (let [ap1 (get-another-p tri1 edge)
        ap2 (get-another-p tri2 edge)
        edges [(->Edge ap1 (:p1 edge)), (->Edge ap1 (:p2 edge))
               (->Edge ap2 (:p1 edge)), (->Edge ap2 (:p2 edge))]
        share-edge (->Edge ap1 ap2)
        {p1 :p1 p2 :p2} edge]
    [(->Triangle ap1 ap2 p1), (->Triangle ap1 ap2 p2), edges, share-edge]))


(defn delaunay-one [tris p]
  "ドロネー三角分割済みの三角形の集合に新たな点を一個追加する"
  (let [tri (some #(if (contains-p? % p) % nil) tris) ;点Pを包含する三角形
        [a b c] (get-points tri)
        tris (-> tris
                 (rm tri)
                 (conj (->Triangle a b p))
                 (conj (->Triangle b c p))
                 (conj (->Triangle c a p))) ;点P追加による新たな三角形を追加
        edges (vec (concat (get-edges tri) [(->Edge a p) (->Edge b p) (->Edge c p)]))];フリップ判定対象の辺
    (loop [edges edges tris tris]
      (if (= 0 (count edges))
        tris
        (let [[edge, edges] [(peek edges), (pop edges)]
              edge-tris (filter #(contains-edge? % edge) tris)]
          (if (not (= (count edge-tris) 2))
            (recur edges tris)
            (let [[tri1 tri2] edge-tris
                  another-p (get-another-p tri2 edge)
                  c-circle (circumscribed-circle tri1)]
              (if (not (contains-p? c-circle another-p))
                (recur edges tris)
                (let [[new-tri1 new-tri2 new-edges new-share-edge] (flip tri1 tri2 edge)
                      tris (-> tris (rm tri1) (rm tri2) (conj new-tri1) (conj new-tri2))
                      edges (vec (concat edges new-edges))]
                  (recur edges tris))))))))))

(defn delaunay-diagram [points]
  "点の集合をドロネー三角分割する"
  (let [super-tri (-> points (circumscribed-rect) (circumscribed-triangle))
        super-points (get-points super-tri)
        all-tris  (loop [points points tris [super-tri]]
                    (if (= 0 (count points))
                      tris 
                      (let [[point rest-points] [(peek points) (pop points)]
                            new-tris (delaunay-one tris point)]
                        (recur rest-points new-tris))))]
    (->> all-tris
        (filter (fn [tri] 
                  (let [tmp-points (get-points tri)
                        points-set (set (concat tmp-points super-points))] 
                    (= 6 (count points-set))))))))