牛顿分形的LISP程序

highflybird 周四, 01/21/2021 - 21:17

牛顿分形简介：

在复数域上使用牛顿迭代生成分形图像，对函数公式\(F(z) = z^3 – 1\)在复数域上面有三个根，一个是1,另外两个是复数:

$$
\left \{
\begin{array}{ll}
Z_1=1 \\
Z_2=\frac {-1+\sqrt{3}i}{2} \\
Z_3=\frac {-1-\sqrt{3}i}{2}
\end{array}
\right.
$$

根据计算出来根的值不同转换为RGB三种不同的颜色，根据迭代次数的多少设置颜色值的大小，即颜色强度。

牛顿分形

下面是用LISP代码的实现。

;;;************
;;;牛顿迭代分形
;;;************
(defun C:Nt (/ con eps  OLDCMD oldsnap t0  gg ita)
  (vl-load-com)
  (setq *APP (vlax-get-acad-object))
  (setq colObj (vla-getinterfaceobject *APP "AutoCAD.AcCmColor.16"))
  (setq OLDCMD (getvar "CMDECHO"))
  (setq oldsnap (getvar "OSMODE"))
  (princ "\n请选择三种不同的颜色:")
  (princ "\n选择颜色1:\n")
  (setq c1 (acad_truecolordlg 1))
  (princ "\n选择颜色2:\n")
  (setq c2 (acad_truecolordlg 3))
  (princ "\n选择颜色3:\n")
  (setq c3 (acad_truecolordlg 5))
  (initget 7)
  (setq gg (getint "\n输入颜色梯度:"))
  (initget 7)
  (setq sol (getint "\n输入分辨率:"))
  (initget 7)
  (setq EPS (getreal "\n输入迭代误差:"))     ;迭代误差
  (initget 7)
  (setq ita (getint "\n最大迭代次数:"))
  (if (and c1 c2 c3 gg sol EPS ita)
    (progn   
      (setvar "CMDECHO" 0)
      (setvar "OSMODE" 0)
      (setq t0 (getvar "TDUSRTIMER"))
      (NewTon_fractal sol sol EPS c1 c2 c3 gg ita)
      (princ "\n画牛顿分形用时")
      (princ (* (- (getvar "TDUSRTIMER") t0) 86400))
      (princ "秒")
      (command ".ZOOM" "E")
      (setvar "OSMODE" oldsnap)
      (setvar "CMDECHO" OLDCMD)
    )
    (alert "n你没有输入有效选择!")
  )
  (gc)
  (princ)
)
 
 
 
(defun NewTon_fractal (solx soly EPS c1 c2 c3 grad ita /
  col1 col2 col3 hsl1 hsl2 hsl3 con xx yy nx ny
  cm newx newy x y k)
  (setq col1 (trans->RGB c1))  
  (setq HSL1 (RGB->HSL (car col1) (cadr col1) (caddr col1)))
  (setq col2 (trans->RGB c2))  
  (setq HSL2 (RGB->HSL (car col2) (cadr col2) (caddr col2)))
  (setq col3 (trans->RGB c3))  
  (setq HSL3 (RGB->HSL (car col3) (cadr col3) (caddr col3)))
  (setq CON (/ (sqrt 3.0) 2))
  (setq xx (/ solX 4.0))
  (setq yy (/ solY 4.0))
  (setq nx (/ solx -2))
  (repeat solx
    (setq ny (/ soly -2))
    (repeat soly
      (setq x (/ nx xx))
      (setq y (/ ny yy))
      (setq k 0)
      (while (and (< k ita) (/= (dist x y) 0))
        (setq cm   (* 3 (dist x y) (dist x y))
              newx (+ (* 2 x (/ 1.0 3)) (/ (- (* x x) (* y y)) cm))
              newy (- (* 2 y (/ 1.0 3)) (/ (* 2 x y) cm))
              x    newx
              y    newy
        )
        (cond
          ( (< (dist (1- x) y) EPS)
            (putcolor hsl1 nx ny k)
            (setq k ita)
          )
          ( (< (dist (+ x 0.5) (- y con)) EPS)
            (putcolor hsl2 nx ny k)
            (setq k ita)
          )  
          ( (< (dist (+ x 0.5) (+ y con)) EPS)
            (putcolor hsl3 nx ny k)
            (setq k ita)
          )
        )
        (setq k (1+ k))
      )
      (setq ny (1+ ny))
    )
    (setq nx (1+ nx))
  )
)
 
;;;==========================
;;; 距离平方函数             
;;;==========================
(defun dist (a b)
  (+ (* a a) (* b b))
)
 
;;;==========================
;;; 颜色转换函数             
;;;==========================
(defun trans->RGB (color / ccc) 
  (if (setq ccc (cdr (assoc 420 color)))
    (Number->RGB ccc)
    (Index->RGB colObj (cdr (assoc 62 color)))
  )
)
 
;;;==========================
;;;用entmake方法画像素点函数 
;;;==========================
(defun putcolor (HSL nx ny k / nH ncolor)
  (setq nH (rem (+ (car HSL) (* grad k)) 361))
  (setq ncolor (hsl->rgb nH (cadr hsl) (caddr hsl)))
  (setq ncolor (rgb->number ncolor))
  (putpixel nx ny ncolor)
)
 
;;;==========================
;;;用entmake方法画像素点函数 
;;;==========================
(defun putpixel (a b c)
  (entmake
    (list
      '(0 . "LWPOLYLINE")    
      '(100 . "AcDbEntity")
      '(100 . "AcDbPolyline")
      '(90 . 2)
      '(43 . 1.0)
      (cons 420 c)
      (cons 10 (list a b))
      (cons 10 (list (1+ a) b))
    )
  )
)
 
;;;==========================
;;;HSL值转RGB值                 
;;;返回RGB值的列表           
;;;==========================
 
;;;==========================
;;;Hue转RGB                  
;;;==========================     
(defun Hue->rgb (v1 v2 vHue / vH)
  (cond
    ( (< vHue 0) (setq vH (1+ vHue)))
    ( (> vHue 1) (setq vH (1- vHue)))
    (t (setq vH vHue))
  )
  (cond
    ( (< (* 6 vH) 1) (+ v1 (* (- v2 v1) 6 vH)))
    ( (< (* 2 vH) 1) v2)
    ( (< (* 3 vH) 2) (+ v1 (* (- v2 v1) 6 (- 0.66666667 vH))))
    (t v1)
  ) 
)
 
;;;==========================
;;;HSL转RGB                  
;;;==========================
(defun Hsl->rgb (Hue Saturation Light / h s l r g b var2 var1)
  (setq h (/ Hue 360.0)
        s (/ Saturation 100.0)
        l (/ Light 100.0)
  )
  (if (= s 0)
    (setq r (* l 255)
          g (* l 255)
          b (* l 255)
    )
    (setq var2 (if (< l 0.5)
                 (* l (1+ s))
                 (+ l s (* s l -1))
               )
          var1 (- (* 2 l) var2)
          r (* 255 (Hue->RGB var1 var2 (+ h 0.33333333)))
          g (* 255 (Hue->RGB var1 var2 h))
          b (* 255 (Hue->RGB var1 var2 (- h 0.33333333)))
    )
  )
  (list (fix r) (fix g) (fix b))
)
 
;;;==========================
;;;RGB值转HSL值              
;;;返回HSL值的列表           
;;;==========================
(defun RGB->HSL(R G B / 
  var_R var_G var_B var_min var_max
  del_max del_R del_G del_B H L S)
  (setq var_R (/ R 255.0))
  (setq var_G (/ G 255.0))
  (setq var_B (/ B 255.0))
  (setq var_min (min var_R var_G var_B))
  (setq var_max (max var_R var_G var_B))
  (setq del_max (- var_max var_min))
  (setq L (/ (+ var_max var_min) 2))
  (if (= del_max 0)
    (setq H 0 S 0)
    (progn 
      (setq S (if (< L 0.5)
                (/ del_max (+ var_max var_min))
                (/ del_max (- 2 var_max var_min))
              )
            del_R (/ (+ (/ (- var_max var_R) 6)  (/ del_Max 2 )) del_max)
            del_G (/ (+ (/ (- var_max var_G) 6)  (/ del_Max 2 )) del_max)
            del_B (/ (+ (/ (- var_max var_B) 6)  (/ del_Max 2 )) del_max)
      )
      (cond
        ( (= var_R var_max)
          (setq H (- del_B del_G))
        )
        ( (= var_G var_max)
          (setq H (+ (/ 1.0 3) del_R (- del_B)))
        )
        ( (= var_B var_max)
          (setq H (+ (/ 2.0 3) del_G (- del_R)))
        )
      )
      (cond
        ( (< H 0) (setq  H (1+ H)))
        ( (> H 1) (setq  H (1- H)))
      )
    )
  )
  (setq h (* 360 H)
        S (* 100 S)
        l (* 100 l) 
  )
  (list (fix H) (fix S) (fix L))
)
 
;;;==========================
;;;把truecolordlg            
;;;420构成的数值返           
;;;回RGB列表.                
;;;==========================
(defun Number->RGB (C)
  (list (lsh C -16)
        (lsh (lsh C 16) -24)
        (lsh (lsh C 24) -24)
  )
)
 
;;;==========================
;;;把truecolordlg            
;;;420构成的数值返           
;;;回RGB列表.                
;;;==========================
(defun RGB->Number (lst)
  (+ (lsh (car lst) 16) (lsh (cadr lst) 8) (caddr lst))
)
 
;;;==========================
;;;RGB转化成索引号           
;;;==========================
(defun RGB->Index (colorObj colRGB / )
  (vla-setRGB colorobj (car colRGB) (cadr colRGB) (caddr colRGB)) 
  (vla-get-ColorIndex colorobj)
)
 
;;;==========================
;;;索引号转化成RGB           
;;;==========================
(defun Index->RGB (colorobj ci / )
  (vla-put-ColorIndex  colorobj ci)
  (list (vla-get-red   colorobj)
        (vla-get-green colorobj)
        (vla-get-blue  colorobj)
  )
)

标签

牛顿 分形
编程