返回列表 回复 发帖
精解代码:
aug01-hilbert-spline2 {
init:
  float x=0.0
  float y=0.0
  float r=0.0
  float rmin=1e20
  float u=0.0
  float v=0.0
  int iter=0
  int nlast=0
  int n=0
  int fac3=0
  int fac4=0
  int power4=0
  zh=(0,0)
;
; set up endpoints for the lines
;
  float h=0.5
  float ooh=1/h
  float ooomh=1/(1-h)
  float hoomh=h/(1-h)
  z0a=@enterexit
  z0b=flip(@enterexit)
  z1=@centerll
  z2=@centerul
  z3=@centerur
  z4=@centerlr
  z5=1+flip(@enterexit)
loop:
final:
  iter=0
  zh=#z*0.25+(0.5,0.5)
  while(iter<@niter)
    iter=iter+1
    x=real(zh)
    y=imag(zh)
;
; lower left sub-square:  shrink, flip horizontally, rotate by -90 degrees
;
    if((x<h)&&(y<h))
      nlast=3
      u=ooh*y
      v=ooh*x
;
; upper left sub-square:  just shrink
;
    elseif((x<h)&&(y>=h))
      nlast=2
      u=ooh*x
      v=ooomh*y-hoomh
;
; upper right sub-square:  shrink & flip horizontally
;
    elseif((x>=h)&&(y>=h))
      nlast=1
      u=ooomh-ooomh*x
      v=ooomh*y-hoomh
;
; lower right sub-square:  shrink, rotate by 90 degrees
;
    elseif((x>=h)&&(y<h))
      nlast=4
      u=ooh*y
      v=ooomh-ooomh*x
    else
      nlast=-1
      u=x
      v=y
    endif
    zh=u+flip(v)
    n=n*4+nlast-1
  endwhile
;
; determine how to start line 0-1
;
  z0=z0a

  iter=0
  fac3=2
  fac4=1
  power4=4
  while(iter<@niter)
    iter=iter+1
    if(iter%2==1)
      if(n%power4==fac3)
        z0=z0b

      endif
    else
      if(n%power4==fac3)
        z0=z0a

      endif
    endif
    power4=4*power4
    fac4=4*fac4
    fac3=fac3+2*fac4
  endwhile
;
; spline from enter to lower left sub-square
;
  zz0=z0
  zz2=z1


  zz1=(zz0+zz2)/2

  c=zz0-zh
  b=2*(zz1-zz0)

    root1=-c/b


  x=real(root1)
  y=imag(root1)
  if(x<0.0)
    r=cabs(root1)
  elseif(x>1.0)
    r=cabs(root1-1)
  else
    r=abs(y)
  endif

  if(r<rmin)
    rmin=r
  endif



;
; spline from lower left to upper left sub-square
;
  zz0=z1
  zz2=z2


    zz1=(zz0+zz2)/2

  c=zz0-zh
  b=2*(zz1-zz0)


    root1=-c/b


  x=real(root1)
  y=imag(root1)
  if(x<0.0)
    r=cabs(root1)
  elseif(x>1.0)
    r=cabs(root1-1)
  else
    r=abs(y)
  endif

  if(r<rmin)
    rmin=r
  endif



;
; spline from upper left to upper right sub-square
;
  zz0=z2
  zz2=z3


    zz1=(zz0+zz2)/2

    c=zz0-zh
  b=2*(zz1-zz0)


    root1=-c/b

  x=real(root1)
  y=imag(root1)
  if(x<0.0)
    r=cabs(root1)
  elseif(x>1.0)
    r=cabs(root1-1)
  else
    r=abs(y)
  endif

  if(r<rmin)
    rmin=r
  endif



; spline from upper right to lower right sub-square
;
  zz0=z3
  zz2=z4

    zz1=(zz0+zz2)/2

  c=zz0-zh
  b=2*(zz1-zz0)

    root1=-c/b


  x=real(root1)
  y=imag(root1)
  if(x<0.0)
    r=cabs(root1)
  elseif(x>1.0)
    r=cabs(root1-1)
  else
    r=abs(y)
  endif

  if(r<rmin)
    rmin=r
  endif



; spline from lower right sub-square to exit
;
  zz0=z4
  zz2=z5

    zz1=(zz0+zz2)/2

  c=zz0-zh
  b=2*(zz1-zz0)


    root1=-c/b

  x=real(root1)
  y=imag(root1)
  if(x<0.0)
    r=cabs(root1)
  elseif(x>1.0)
    r=cabs(root1-1)
  else
    r=abs(y)
  endif

  if(r<rmin)
    rmin=r
  endif



;
  #index=rmin
default:
  title="Hilbert spline 2"
  param niter
    caption="iterations"
    default=0
    min=0
  endparam
  param centerll
    caption="lower left center"
    default=(0.25,0.25)
    hint="Center of the lower left sub-square. Make both coordinates \
    between 0 & 1; use (0.25,0.25) for standard Hilbert curve."
  endparam
  param centerul
    caption="upper left center"
    default=(0.25,0.75)
    hint="Center of the upper left sub-square. Make both coordinates \
    between 0 & 1; use (0.25,0.75) for standard Hilbert curve."
  endparam
  param centerur
    caption="upper right center"
    default=(0.75,0.75)
    hint="Center of the upper right sub-square. Make both coordinates \
    between 0 & 1; use (0.75,0.75) for standard Hilbert curve."
  endparam
  param centerlr
    caption="lower right center"
    default=(0.75,0.25)
    hint="Center of the lower right sub-square. Make both coordinates \
    between 0 & 1; use (0.75,0.25) for standard Hilbert curve."
  endparam
  param enterexit
    caption="enter/exit"
    default=0.25
    min=0.0
    max=1.0
    hint="Where the curve enters and exits the block of 4 sub-squares. \
      Between 0 & 1; use 0.25 for standard Hilbert curve."
  endparam


}
未命名.jpg
2013-6-13 22:31
未命名.jpg
2013-6-13 22:33

Hilbert spline 2.gsp (45.55 KB)
本人调色始终差炎候,UF中的色彩真艳:
Fractal2.JPG
2013-6-13 22:34
调一幅:
未命名.jpg
2013-6-13 22:37
将原M集的数据粘到#13中,将M集的迭代终点坐标加载到原希氏数据中,只需花几分钟时长。扫得一图:
未命名.jpg
2013-6-10 11:23
未命名.jpg
2013-6-13 22:39
M集作陷阱的UF代码,精简,整理,编辑,几处代码动的手脚比较大,发现图形没什么两样:
juliatrap(BOTH) {
; By Samuel Monnier, 1999
init:
   count = 0
  float stdz = 0.0
  float logp = 1/log(2)
  float logb = log(log(1e20))
    float dz = 0.0
loop:
  trz = exp(flip(-pi/180*@rot))*(#z-@center)/@size
  x = 1/sqrt(@r)*real(trz)
  y = sqrt(@r)*imag(trz)
    trz = x + flip(y)
  zz = trz
  float i = 0
  iter = @niter


    start = zz
    zz = @seed


  while i < @niter
    i = i + 1
    zz = zz^2 + start
    if |zz| > 1e20
      iter = i
      i = @niter
    endif

  endwhile
  dz = real(iter + logp*logb - logp*log(log(cabs(zz))))
    if dz>stdz
      stdz = dz


    endif


final:
  #index = .1*stdz
default:
  title = "Julia Trap"
  helpfile = "sam-help/juliatrap.htm"
  param seed
    caption = "Seed (for 'custom')"
    default = (0,0)
  endparam
  param center
    caption = "Center"
    default = (1,1)
  endparam
  param rot
    caption = "Rotation"
    default = -30.0
  endparam
  param size
    caption = "Size"
    default = 0.386750
  endparam
  param r
    caption = "Ratio Width/Heigh"
    default = 1.0
  endparam
  param @niter
    caption = "Julia Iterations"
    default = 100
  endparam
  param freq
    caption = "Trap Frequency"
    default = 0.0
  endparam
}
Fractal2.JPG
2013-6-13 22:42

红色部分是我动过手脚比较大的部分,原代码是:
if first
      stdz = dz
      first = false
    elseif dz > stdz
      stdz = dz
    endif
这段代码如读天书,所以改成红色部分后,效果一样。红色部分是什么意思?
没考虑红色部分,用画板则得效果如下:
未命名.JPG
2013-6-13 22:46

等势圈遮住了M集部分,图看起来不舒服。在UF中拿掉这部分代码,则与GSP一致!这段代码是什么意思?如何处理?dz中的iter,其实就是充当陷阱M集的et。
8# changxde


常老师,把迭代次数设为30,放大看看有没有断头现象。我在电脑上试试,次数不大时显示正常,但当次数超过30时就出现问题。
返回列表