scaling 2D gray map by bi-linear interpolation

bi-linear.clj

; 2012.03.23 境界の処理がイマイチだったので全面書き換え
; clojure.math.numeric-tower (https://github.com/clojure/math.numeric-tower)
; が必要.
; clojure 1.3 向け

(use '[clojure.math.numeric-tower :only (floor)])

; 長さ w のベクタの, 長さ h のベクタとして与えられる
; 解像度 w * h の値マップ src を双線形補完を用いて
; 解像度 W * H に拡大縮小し,
; 長さ W のベクタの, 長さ H のベクタとして返す.
(defn bi-linear [src W H]
  (let [w (count (first src))
        h (count src)
        srcij (fn [i j] (nth (nth src (mod j h)) (mod i w)))
        abc (fn [a b] (- (* 2 b) a))
        ijij ; i1 j1 i2 j2 を与えると
             ; v(i1, j1) から v(i2 j2) に引いた直線の
             ; (2*i2 - i1, 2*j2 - j1) での値を求む
             (fn [i1 j1 i2 j2] (abc (srcij i1 j1) (srcij i2 j2)))
        exsrc (reduce
               (fn [acc j]
                 (assoc acc j
                   (reduce
                     (fn [acc2 i]
                       (assoc acc2 i
                         (cond (< j 0)
                                 (cond (< i 0) (ijij  1  1  0  0)
                                       (< i w) (ijij  i  1  i  0)
                                       :else   (ijij -2  1 -1  0))
                               (< j h)
                                 (cond (< i 0) (ijij  1  j  0  j)
                                       (< i w) (srcij i  j)
                                       :else   (ijij -2  j -1  j))
                               :else
                                 (cond (< i 0) (ijij  1 -2  0 -1)
                                       (< i w) (ijij  i -2  i -1)
                                       :else   (ijij -2 -2 -1 -1))
                         )
                       )
                     )
                     {}
                     (range -1 (+ w 1)))))
               {}
               (range -1 (+ h 1)))
        ]
    (apply vector
      (for [J (range H)]
        (apply vector 
          (for [I (range W)]
            (let [;方形領域 a (0 <= x <= w, 0 <= y <= h) を
                  ;方形領域 A (0 <= X <= W, 0 <= Y <= H) に拡縮する場合
                  ;                2*I + 1  2*J + 1
                  ;A において座標 (-------, -------) である点の
                  ;                   2*W     2*H
                  ;a における座標 (x, y)
                  x (/ (* w (+ (* 2.0 I) 1.0)) (* 2.0 W))
                  y (/ (* h (+ (* 2.0 J) 1.0)) (* 2.0 H))
                  ; 元のデータ src のインデクス i, j の値は
                  ; 方形領域 a の座標 (i + 0.5, j + 0.5) の値を
                  ; 表していると考える.
                  ; 補完に使う元の空間の四つの点うち左上の点の
                  ; 拡張した元データにおけるインデクス (xi, yj)
                  xi (long (floor (- x 0.5)))
                  yj (long (floor (- y 0.5)))
                  ; 四つの点の座標 (x0, y0) (x1, y0) (x0, y1) (x1, y1)
                  x0 (+ xi 0.5)
                  y0 (+ yj 0.5)
                  x1 (+ xi 1.5)
                  y1 (+ yj 1.5)
                  ; 四つの点の値 v(x0, y0), v(x1, y0), v(x0, y1), v(x1, y1)
                  vx0y0 ((exsrc      yj )      xi )
                  vx1y0 ((exsrc      yj ) (inc xi))
                  vx0y1 ((exsrc (inc yj))      xi )
                  vx1y1 ((exsrc (inc yj)) (inc xi))
                  ; y0 における x 方向の値直線の x における値 v(x, y0)
                  ;        vx1y0 - vx0y0
                  ; vxy0 = ------------- (x - x0) + vx0y0
                  ;           x1 - x0
                  vxy0 (+ (/ (* (- vx1y0 vx0y0) (- x x0)) (- x1 x0)) vx0y0)
                  ; y1 における x 方向の値直線の x における値 v(x, y1)
                  ;        vx1y1 - vx0y1
                  ; vxy1 = ------------- (x - x0) + vx0y1
                  ;           x1 - x0
                  vxy1 (+ (/ (* (- vx1y1 vx0y1) (- x x0)) (- x1 x0)) vx0y1)
                  ; x における y 方向の値直線の y における値 v(x, y)
                  ;       vxy1 - vxy0
                  ; vxy = ----------- (y - y0) + vxy0
                  ;         y1 - y0
                  vxy (+ (/ (* (- vxy1 vxy0) (- y y0)) (- y1 y0)) vxy0)
                  ]
              vxy
              )))))
  ))