View of /lisp/curveto.l

Parent Directory

Revision Log

Initial revision

;
; アレイの定義
;
(defun array-rank (a)
  (cond ((vectorp a) (1+ (array-rank (vref a 0))))
	(t 0)))
;
(defun array-dimension (a dim)
  (lets ((tmparray a))
	(do ((i 0 (1+ i)))
	    ((= i dim))
	    (setq tmparray (vref tmparray 0)))
	(vector-length tmparray)))
;
(defun array-dimensions (a)
  (lets ((dimlist))
	(do ((tmparray a (vref tmparray 0)))
	    ((not (vectorp tmparray)))
	    (push (vector-length tmparray) dimlist))
	(reverse dimlist)))
;
;(defun make-array (dimlist)
;  (lets ((tmpvec))
;	(do ((rev (reverse dimlist) (cdr rev)))
;	    ((atom rev))
;	    (print rev)
;	    (setq tmpvec (vector (car rev) `(,tmpvec)))
;	    (print-vec tmpvec)
;	    )
;	tmpvec))
;



;
; ベクトルの加算
;
(defun vplus (v1 v2)
  (cond ((= (vector-length v1) (vector-length v2))
	 (lets ((l (vector-length v1))(v (vector l)))
	       (do ((i 0 (1+ i)))
		   ((= i l))
		   (vset v i (plus (vref v1 i)(vref v2 i))))
	       v))
	(t nil)))
;
; ベクトルの減算
;
(defun vminus (v1 v2)
  (cond ((= (vector-length v1)(vector-length v2))
	 (lets ((l (vector-length v1))(v (vector l)))
	       (do ((i 0 (1+ i)))
		   ((= i l))
		   (vset v i (difference (vref v1 i)(vref v2 i))))
		   v))
	(t nil)))
;
; ベクトルの加算（自動拡張あり）
;
(defun vpluse (v1 v2)
  (lets ((l1 (vector-length v1))(l2 (vector-length v2))
	 (v (vector (max l1 l2))))
	(do ((i 0 (1+ i)))
	    ((= i (min l1 l2)))
	    (vset v i (plus (vref v1 i)(vref v2 i))))
	(cond ((> l1 l2)
	       (do ((i l2 (1+ i)))
		   ((= i l1))
		   (vset v i (vref v1 i))))
	      (t (do ((i l1 (1+ i)))
		     ((= i l2))
		     (vset v i (vref v2 i)))))
	v))
;
; ベクトルの減算（自動拡張あり）
;
(defun vminuse (v1 v2)
  (lets ((l1 (vector-length v1))(l2 (vector-length v2))
	 (v (vector (max l1 l2))))
	(do ((i 0 (1+ i)))
	    ((= i (min l1 l2)))
	    (vset v i (difference (vref v1 i)(vref v2 i))))
	(cond ((> l1 l2)
	       (do ((i l2 (1+ i)))
		   ((= i l1))
		   (vset v i (vref v1 i))))
	      (t (do ((i l1 (1+ i)))
		     ((= i l2))
		     (vset v i (minus (vref v2 i))))))
	v))
;
; ベクトルのスカラー倍
;
(defun vtimes (v k)
  (mapvector v #'(lambda (e)(times e k))))
;
; ベクトルの要素ごとの積
;
(defun vproduct (v1 v2)
  (cond ((= (vector-length v1)(vector-length v2))
	 (lets ((l (vector-length v1))(v (vector l)))
	       (do ((i 0 (1+ i)))
		   ((= i l))
		   (vset v i (times (vref v1 i)(vref v2 i))))
	       v))
	(t nil)))
;
; ベクトルの内積
;
(defun vip (v1 v2)
  (cond ((= (vector-length v1)(vector-length v2))
	 (do ((i 0 (1+ i))
	      (l (vector-length v1))
	      (sum 0 (plus sum (times (vref v1 i)(vref v2 i)))))
	     ((= i l) sum)))
	(t nil)))
;
; 行列（＝ベクトルのベクトル）のある要素へのアクセス
;
(defun matrix (v i j)
  (cond ((and (<= 0 i)(> (vector-length v) i))
	 (cond ((and (<= 0 j)(> (vector-length (vref v i)) j))
		(vref (vref v i) j))
	       (t nil)))
	(t nil)))
;
; ４元連立一次方程式を解く（ピボットなし）
;
(defun solve4 (a1 a2 a3 a4 b)
  (lets ((a11 (float(car a1)))(a12 (float(cadr a1)))(a13 (float(caddr a1)))(a14 (float(cadddr a1)))
	 (a21 (float(car a2)))(a22 (float(cadr a2)))(a23 (float(caddr a2)))(a24 (float(cadddr a2)))
	 (a31 (float(car a3)))(a32 (float(cadr a3)))(a33 (float(caddr a3)))(a34 (float(cadddr a3)))
	 (a41 (float(car a4)))(a42 (float(cadr a4)))(a43 (float(caddr a4)))(a44 (float(cadddr a4)))
	 (b1 (float(car b)))(b2 (float(cadr b)))(b3 (float(caddr b)))(b4 (float(cadddr b)))
; 1st
	 (a12 (//$ a12 a11))
	 (a13 (//$ a13 a11))
	 (a14 (//$ a14 a11))
	 (b1 (//$ b1 a11))
	 (a22 (-$ a22 (*$ a21 a12)))
	 (a23 (-$ a23 (*$ a21 a13)))
	 (a24 (-$ a24 (*$ a21 a14)))
	 (b2 (-$ b2 (*$ a21 b1)))
	 (a32 (-$ a32 (*$ a31 a12)))
	 (a33 (-$ a33 (*$ a31 a13)))
	 (a34 (-$ a34 (*$ a31 a14)))
	 (b3 (-$ b3 (*$ a31 b1)))
	 (a42 (-$ a42 (*$ a41 a12)))
	 (a43 (-$ a43 (*$ a41 a13)))
	 (a44 (-$ a44 (*$ a41 a14)))
	 (b4 (-$ b4 (*$ a41 b1)))
; 2nd
	 (a23 (//$ a23 a22))
	 (a24 (//$ a24 a22))
	 (b2 (//$ b2 a22))
	 (a13 (-$ a13 (*$ a12 a23)))
	 (a14 (-$ a14 (*$ a12 a24)))
	 (b1 (-$ b1 (*$ a12 b2)))
	 (a33 (-$ a33 (*$ a32 a23)))
	 (a34 (-$ a34 (*$ a32 a24)))
	 (b3 (-$ b3 (*$ a32 b2)))
	 (a43 (-$ a43 (*$ a42 a23)))
	 (a44 (-$ a44 (*$ a42 a24)))
	 (b4 (-$ b4 (*$ a42 b2)))
; 3rd
	 (a34 (//$ a34 a33))
	 (b3 (//$ b3 a33))
	 (a14 (-$ a14 (*$ a13 a34)))
	 (b1 (-$ b1 (*$ a13 b3)))
	 (a24 (-$ a24 (*$ a23 a34)))
	 (b2 (-$ b2 (*$ a23 b3)))
	 (a44 (-$ a44 (*$ a43 a34)))
	 (b4 (-$ b4 (*$ a43 b3)))
; 4th
	 (b4 (//$ b4 a44))
	 (b1 (-$ b1 (*$ a14 b4)))
	 (b2 (-$ b2 (*$ a24 b4)))
	 (b3 (-$ b3 (*$ a34 b4))))
	(list b1 b2 b3 b4)))
;
; ４元連立一次方程式を解く（ピボットあり）
;
(defun solve4p (a1 a2 a3 a4 b)
  (lets ((m1 (vector 5 `(,(float (car a1)),(float (cadr a1)),(float (caddr a1)),(float (cadddr a1)),(float (car b)))))
	 (m2 (vector 5 `(,(float (car a2)),(float (cadr a2)),(float (caddr a2)),(float (cadddr a2)),(float (cadr b)))))
	 (m3 (vector 5 `(,(float (car a3)),(float (cadr a3)),(float (caddr a3)),(float (cadddr a3)),(float (caddr b)))))
	 (m4 (vector 5 `(,(float (car a4)),(float (cadr a4)),(float (caddr a4)),(float (cadddr a4)),(float (cadddr b)))))
	 (m (vector 4 (list m1 m2 m3 m4))))
	(do ((i 0 (1+ i))
	     (v (vector 5))
	     (p -1))
	    ((>= i 4))
	    (do ((ii (1+ i) (1+ ii)))
		((>= ii 4))
		(cond ((>$ (abs (vref (vref m ii) i))(abs (vref (vref m i) i)))
		       (setq p ii))))
	    (cond ((> p i)
;		   (format "p =/c is larger than i =/c, Swapping!!/n" p i)
		   (setq v (vref m i))
		   (vset m i (vref m p))
		   (vset m p v)))
	    (vset m i (vtimes (vref m i)(//$ 1.0 (vref (vref m i) i))))
	    (do ((ii 0 (1+ ii)))
		((>= ii i))
		(vset m ii
		      (vminus (vref m ii) (vtimes (vref m i)(vref (vref m ii) i)))))
	    (do ((ii (1+ i) (1+ ii)))
		((>= ii 4))
		(vset m ii
		      (vminus (vref m ii) (vtimes (vref m i)(vref (vref m ii) i)))))
;	    (do ((l 0 (1+ l)))
;		((>= l 4))
;		(print-vec (vref m l)))
	    )
	(list (vref (vref m 0) 4)(vref (vref m 1) 4)(vref (vref m 2) 4)(vref (vref m 3) 4))))
;
; 幅のあるベジェ曲線（ピボットなし）
;
(defun curveto (z1 z2 z3 z4 w)
  (lets ((l (times2 0.25 (plus2 (plus2 z1 z2) (plus2 z2 z3))))
	 (r (times2 0.25 (plus2 (plus2 z2 z3) (plus2 z3 z4))))
	 (v1 (diff2 z2 z1))
	 (v2 (diff2 z3 z4))
	 (zi1 (plus2 z1 (times2 w (norm2 (rot90 v1)))))
	 (zo1 (plus2 z1 (times2 w (norm2 (rot270 v1)))))
	 (zi4 (plus2 z4 (times2 w (norm2 (rot270 v2)))))
	 (zo4 (plus2 z4 (times2 w (norm2 (rot90 v2)))))
	 (c (times2 0.5 (plus2 l r)))
	 (vc (diff2 r l))
	 (ci (plus2 c (times2 w (norm2 (rot90 vc)))))
	 (co (plus2 c (times2 w (norm2 (rot270 vc)))))
	 (l1 (list (times 2 (car v1)) (car v2) (times -4 (car vc)) 0))
	 (l2 (list (times 2 (cadr v1)) (cadr v2) (times -4 (cadr vc)) 0))
	 (l3 (list (car v1) (times 2 (car v2)) 0 (times -4 (car vc))))
	 (l4 (list (cadr v1) (times 2 (cadr v2)) 0 (times -4 (cadr vc))))
	 (bi (list (plus (times -3 (car zi1))(times -1 (car zi4))(times 4 (car ci)))
		   (plus (times -3 (cadr zi1))(times -1 (cadr zi4))(times 4 (cadr ci)))
		   (plus (times -1 (car zi1))(times -3 (car zi4))(times 4 (car ci)))
		   (plus (times -1 (cadr zi1))(times -3 (cadr zi4))(times 4 (cadr ci)))))
	 (bo (list (plus (times -3 (car zo1))(times -1 (car zo4))(times 4 (car co)))
		   (plus (times -3 (cadr zo1))(times -1 (cadr zo4))(times 4 (cadr co)))
		   (plus (times -1 (car zo1))(times -3 (car zo4))(times 4 (car co)))
		   (plus (times -1 (cadr zo1))(times -3 (cadr zo4))(times 4 (cadr co)))))
	 (ansi (solve4 l1 l2 l3 l4 bi))
	 (anso (solve4 l1 l2 l3 l4 bo)))
	(list (list zi1 (list (plus (car zi1)(times (car ansi)(car v1)))(plus (cadr zi1)(times (car ansi)(cadr v1))))
		    (list (plus (car zi4)(times (cadr ansi)(car v2)))(plus (cadr zi4)(times (cadr ansi)(cadr v2)))) zi4)
	      (list zo1 (list (plus (car zo1)(times (car anso)(car v1)))(plus (cadr zo1)(times (car anso)(cadr v1))))
		    (list (plus (car zo4)(times (cadr anso)(car v2)))(plus (cadr zo4)(times (cadr anso)(cadr v2)))) zo4))))
;
; 幅のあるベジェ曲線(2)（ピボットあり）
;
(defun curvetop (z1 z2 z3 z4 w)
  (lets ((l (times2 0.25 (plus2 (plus2 z1 z2) (plus2 z2 z3))))
	 (r (times2 0.25 (plus2 (plus2 z2 z3) (plus2 z3 z4))))
	 (v1 (diff2 z2 z1))
	 (v2 (diff2 z3 z4))
	 (zi1 (plus2 z1 (times2 w (norm2 (rot90 v1)))))
	 (zo1 (plus2 z1 (times2 w (norm2 (rot270 v1)))))
	 (zi4 (plus2 z4 (times2 w (norm2 (rot270 v2)))))
	 (zo4 (plus2 z4 (times2 w (norm2 (rot90 v2)))))
	 (c (times2 0.5 (plus2 l r)))
	 (vc (diff2 r l))
	 (ci (plus2 c (times2 w (norm2 (rot90 vc)))))
	 (co (plus2 c (times2 w (norm2 (rot270 vc)))))
	 (l1 (list (times -1 (car v1)) (car v2) (times -4 (car vc)) 0))
	 (l2 (list (times -1 (cadr v1)) (cadr v2) (times -4 (cadr vc)) 0))
	 (l3 (list (times 2 (car v1)) (car v2) 0 (times 4 (car vc))))
	 (l4 (list (times 2 (cadr v1)) (cadr v2) 0 (times 4 (cadr vc))))
	 (bi (list (plus (times 2 (car zi1))(times -2 (car zi4)))
		   (plus (times 2 (cadr zi1))(times -2 (cadr zi4)))
		   (plus (times -3 (car zi1))(times -1 (car zi4))(times 4 (car ci)))
		   (plus (times -3 (cadr zi1))(times -1 (cadr zi4))(times 4 (cadr ci)))))
	 (bo (list (plus (times 2 (car zo1))(times -2 (car zo4)))
		   (plus (times 2 (cadr zo1))(times -2 (cadr zo4)))
		   (plus (times -3 (car zo1))(times -1 (car zo4))(times 4 (car co)))
		   (plus (times -3 (cadr zo1))(times -1 (cadr zo4))(times 4 (cadr co)))))
	 (ansi (solve4p l1 l2 l3 l4 bi))
	 (anso (solve4p l1 l2 l3 l4 bo)))
	(list (list zi1 (list (plus (car zi1)(times (car ansi)(car v1)))(plus (cadr zi1)(times (car ansi)(cadr v1))))
		    (list (plus (car zi4)(times (cadr ansi)(car v2)))(plus (cadr zi4)(times (cadr ansi)(cadr v2)))) zi4)
	      (list zo1 (list (plus (car zo1)(times (car anso)(car v1)))(plus (cadr zo1)(times (car anso)(cadr v1))))
		    (list (plus (car zo4)(times (cadr anso)(car v2)))(plus (cadr zo4)(times (cadr anso)(cadr v2)))) zo4))))
;
; 簡単な出力
;
;(drawcurves (curvetop z1 z2 z3 z3 width)) として使う
(defun drawcurves (lines (psfile 'testcurve.ps))
  (lets ((standard-output (outopen (stream psfile)))
	 (zi (car lines))(zo (cadr lines))
	 (zi1 (car zi))(zi2 (cadr zi))(zi3 (caddr zi))(zi4 (cadddr zi))
	 (zo1 (car zo))(zo2 (cadr zo))(zo3 (caddr zo))(zo4 (cadddr zo))
	 (date (date-time)))
	(format "%!/nnewpath/n")
	(format "//X {moveto currentlinewidth 1 setlinewidth 5 5 rmoveto -10 -10 rlineto 10 0 rmoveto -10 10 rlineto stroke setlinewidth} def/n")
	(format "//Times-Roman findfont 20 scalefont setfont/n")
	(format "50 50 moveto (/c    /c-/c-/c /c:/c) show/n"
		psfile (substring date 0 2)(substring date 2 4)
		(substring date 4 6)(substring date 6 8)
		(substring date 8 10))
	(format "100 100 translate/n")
	(format "2 setlinewidth/n")
	(format "/c /c moveto " (fix (car zi1))(fix (cadr zi1)))
	(format "/c /c /c /c /c /c curveto stroke/n" (fix (car zi2))(fix (cadr zi2))(fix (car zi3))(fix (cadr zi3))(fix (car zi4))(fix (cadr zi4)))
	(format "/c /c moveto " (fix (car zo1))(fix (cadr zo1)))
	(format "/c /c /c /c /c /c curveto stroke/n" (fix (car zo2))(fix (cadr zo2))(fix (car zo3))(fix (cadr zo3))(fix (car zo4))(fix (cadr zo4)))
	(format "/c /c X /c /c X /c /c X /c /c X/n" (fix (car zi1))(fix (cadr zi1))(fix (car zi2))(fix (cadr zi2))(fix (car zi3))(fix (cadr zi3))(fix (car zi4))(fix (cadr zi4)))
	(format "/c /c X /c /c X /c /c X /c /c X/n" (fix (car zo1))(fix (cadr zo1))(fix (car zo2))(fix (cadr zo2))(fix (car zo3))(fix (cadr zo3))(fix (car zo4))(fix (cadr zo4)))
	(format "showpage/n")))
;
;
;