簡介

"平均顏色是使用平滑迭代計數的小數部分在平均總和之間進行插值的著色函式族。" ^[1]

這裡，條紋是指垂直於迭代帶（逃逸時間的水平集）的線。這些"徑向線... 遵循外部射線，外部射線在 M 的理論研究和全純動力學中非常重要。" (A Cheritat)^[2]

"基本上是計算角度的一種廉價方法。TIA 對所有迭代中的角度進行平均以獲得平滑的結果。因此，需要進行一些初始化和每次迭代的計算來進行求和。它必須在新增下一個值之前儲存先前迭代中的總和，以便您可以在迭代帶之間插值以獲得連續函式（與通常對連續勢進行插值的方式相同）。這種插值是在迭代退出後完成的。" ^[3] (xenodreambuie - Jux 的作者)

"... 它可能顯示出接近外部射線的曲線，但您仍然沒有 SAC 返回的數字與真實外部引數之間的明確關係。" 沃爾夫·容

名稱

SAM = 條紋平均方法
SAC = 條紋平均著色^[4]
徑向線

歷史

Jussi Härkönen 在 2007 年引入了一種名為條紋平均的新著色方式。它基於三角不等式平均著色的行為。

示例

影像

fractalforums stripe-average-mandelbrot-2 by Zephitmaal

以灰色顯示的曼德勃羅集，影像和原始碼 (processing 和 c)
曼德勃羅集縮放 (wake 1/3) 及 c 原始碼
巴西利卡朱利亞集
不連通的朱利亞集：內爆花椰菜

影片

youtube: 曼德勃羅集著色演算法 by sfract

演算法

對於引數平面矩形中的每個點 c，執行以下操作：^[5]

計算點 z0 的軌道。結果是複數序列 {z0,z1,z2,...}
將函式 t 應用於該序列以獲得新的實數序列：{ t0=t(z0), t1=t(z1), t2=t(z2), ...}.
計算該實數序列的平均值 A ( = 算術平均數)
將實數 A 對映到顏色 ( 著色函式)

軌道

截斷軌道:

${\mathcal {O}}_{n}(z)=\{z,z_{1},z_{2},\ldots z_{n}\}$

跳過軌道是截斷軌道的子集

從序列中跳過前 k+1 個專案
不要使用前 k+1 個專案來計算平均值

${\mathcal {O}}_{n,k}(z)=\{z_{n-k},z_{n-k+1},\ldots z_{n}\}$

加數函式

加數函式 t

將複數 $z$ 對映到實數 $t_{n}$
在跳過的軌道元素的元素上是連續的

t_{n}=t(z_{n})={\frac {sin(s*arg(z_{n}))}{2}}+{\frac {1}{2}}

t\ \colon \mathbb {C} \to \mathbb {R}

其中

s 是條紋密度

平均

全部

截斷軌道的平均值

從截斷軌道的所有點計算平均值

A_{n}=A({\mathcal {O}}_{n})={\frac {1}{n}}\sum _{i=0}^{n}t_{i}={\frac {t_{0}+t_{1}+\cdots +t_{n}}{n}}

如果計算所有點（從 i=1 到 i=n）的平均值，那麼

條紋將可見（好）
逃逸時間的等值線也將可見（不好，不連續）

解決方案是

插值
跳過的截斷軌道

跳過

跳過軌道的平均值

從計算平均值中跳過前 (k+1) 個專案

A_{n,k}={\frac {1}{n-k}}\sum _{i=n-k}^{n}t_{i}={\frac {t_{n-k}+t_{n-k+1}+\cdots +t_{n}}{n-k}}

請注意

$k<n$

$n-k=k+1$

它消除了一些不連續性，但沒有消除逃逸時間的等值線

插值

逃逸時間的等值線仍然可見（不連續）。可以使用插值消除它們。^[6]

插值

在 2 個值之間
使用
- 線性函式 = 線性插值
- 樣條 = 樣條插值
- 多項式

步驟

構建一個平滑的迭代計數 real escape time u = 實數（參見 J Harkonen 論文第 21 頁的公式 3.11）
取其小數部分^[7] d

$u_{n}=u(z_{n})=n+1+{\frac {1}{ln(2)}}ln{\frac {ln(ER)}{ln(abs(z_{n}))}}$

d_{n}=\operatorname {frac} (u_{n})\ =\ u_{n}-\lfloor u_{n}\rfloor \ for\ u_{n}\geq 0

小數（分數）部分 d 為

在先前等值線 $B_{n-1}$ 的邊界處為零
在下一個等級集邊界 $B_{n}$ 的鄰域內收斂到 1

“因此它可以用作插值係數。”

線性

計算線性平均值 $A^{l}$ 在兩個值之間插值

一個是從“完整”系列（截斷並跳過）獲得的 = Avg = $A_{n}$
另一個是從除最後一個元素外的所有 t 個元素的平均值獲得的 = prevAvg = $A_{n-1}$

  $A^{l}=dA_{n}+(1-d)A_{n-1}$

參見 J Harkonen 論文第 28 頁的公式 4.2

或 Syntopia 的 GLSL 程式碼

float mix = frac*avg+(1.0-frac)*prevAvg;

以下是 Syntopia 的 GLSL 程式碼的重要片段

// iteration
for ( i = 0; i < Iterations; i++) {	
	z = complexMul(z,z) + c;	
	if (i>=Skip) {		
		count++;	
		lastAdded = 0.5+0.5*sin(StripeDensity*atan(z.y,z.x));		
		avg +=  lastAdded;
	}
	z2 = dot(z,z);
	if (z2> EscapeRadius2 && i>Skip) break;
}


prevAvg = (avg -lastAdded)/(count-1.0); 

avg = avg/count;

frac =1.0+(log2(log(EscapeRadius2)/log(z2)));	

float mix = frac*avg+(1.0-frac)*prevAvg;

樣條曲線

使用 Catmull-Rom 樣條曲線計算樣條曲線平均值 $A^{s}$

多項式 $H_{i}(d)$

${\begin{array}{lcl}H_{0}={\tfrac {1}{2}}(-d^{2}+d^{3})\\H_{1}={\tfrac {1}{2}}(d+4d^{2}-3d^{3})\\H_{2}={\tfrac {1}{2}}(2-5d^{2}+3d^{3})\\H_{3}={\tfrac {1}{2}}(-d+2d^{2}-d^{3})\end{array}}$

那麼

$A^{s}=H_{0}A_{n}+H_{1}A_{n-1}+H_{2}A_{n-2}+H_{3}A_{n-3}$

參見 J Harkonen 論文第 28 頁的公式 4.3

程式碼

C

M_LN2 在 math.h 中定義^[8]

#define M_LN2 0.69314718055994530942 /* log_e 2 */ The natural logarithm of the 2.^[9]

/* 
   c program: console, 1-file
   samm.c
   
   samm = Stripe Average Method  
   with :
   - skipping first (i_skip+1) points from average
   - linear interpolation 
   
   en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/stripeAC
   
  
   
   --------------------------------
   1. draws Mandelbrot set for complex quadratic polynomial 
   Fc(z)=z*z +c
   using samm = Stripe Average Method  
   -------------------------------         
   2. technique of creating ppm file is  based on the code of Claudio Rocchini
   http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg
   create 24 bit color graphic file ,  portable pixmap file = PPM 
   see http://en.wikipedia.org/wiki/Portable_pixmap
   to see the file use external application ( graphic viewer)
   -----
   it is example  for : 
   https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set
 
   -------------
   compile : 

 
 
   gcc samm.c -lm -Wall
 
 
   ./a.out
   
   
   -------- git --------
   
   
   cd existing_folder
   git init
   git remote add origin git@gitlab.com:adammajewski/mandelbrot_wiki_ACh.git
   git add samm.c
   git commit -m "samm + linear interpolation "
   git push -u origin master



 
 
 
*/
#include <stdio.h>
#include <math.h>
#include <complex.h> // https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c
 
 

#define M_PI        3.14159265358979323846    /* pi */
 
/* screen ( integer) coordinate */
int iX,iY;
const int iXmax = 1200; 
const int iYmax = 1000;

/* world ( double) coordinate = parameter plane*/
double Cx,Cy;
const double CxMin=-2.5;
const double CxMax=1.5;
const double CyMin=-5.0/3.0;
const double CyMax= 5.0/3.0;
/* */
double PixelWidth; //=(CxMax-CxMin)/iXmax;
double PixelHeight; // =(CyMax-CyMin)/iYmax;

/* color component ( R or G or B) is coded from 0 to 255 */
/* it is 24 bit color RGB file */
const int MaxColorComponentValue=255; 
FILE * fp;
char *filename="samm_i.ppm"; // https://commons.wikimedia.org/wiki/File:Mandelbrot_set_-_Stripe_Average_Coloring.png
char *comment="# ";/* comment should start with # */
        
static unsigned char color[3]; // 24-bit rgb color

unsigned char s = 5; // stripe density

const int IterationMax=1000; //  N in wiki 
int i_skip = 1; // exclude (i_skip+1) elements from average

/* bail-out value for the bailout test for exaping points
   radius of circle centered ad the origin, exterior of such circle is a target set  */
const double EscapeRadius=10000; // = ER big !!!!
double lnER; // ln(ER)
        
        
        
        
double complex give_c(int iX, int iY){
  double Cx,Cy;
  Cy=CyMin + iY*PixelHeight;
  if (fabs(Cy)< PixelHeight/2) Cy=0.0; /* Main antenna */
  Cx=CxMin + iX*PixelWidth;
   
  return Cx+Cy*I;
 
 
}
 


// the addend function
// input : complex number z
// output : double number t 
double Give_t(double complex z){

  return 0.5+0.5*sin(s*carg(z));

}

/*
  input :
  - complex number
  - intege
  output = average
 
*/
double Give_Arg(double complex C , int iMax)
{
  int i=0; // iteration 
   
   
  double complex Z= 0.0; // initial value for iteration Z0
  double A = 0.0; // A(n)
  double prevA = 0.0; // A(n-1)
  double R; // =radius = cabs(Z)
  double d; // smooth iteration count
   
    
  // iteration = computing the orbit
  for(i=0;i<iMax;i++)
    {  
      Z=Z*Z+C; // https://wikibook.tw/wiki/Fractals/Iterations_in_the_complex_plane/qpolynomials
      
      if (i>i_skip) A += Give_t(Z); // 
      
      R = cabs(Z);
      if(R > EscapeRadius) break; // exterior of M set
   
      prevA = A; // save value for interpolation
        
    } // for(i=0
   
   
  if (i == iMax) 
    A = -1.0; // interior 
  else { // exterior
    // computing interpolated average
    A /= (i - i_skip) ; // A(n)
    prevA /= (i - i_skip - 1) ; // A(n-1) 
    // smooth iteration count
    d = i + 1 + log(lnER/log(R))/M_LN2;
    d = d - (int)d; // only fractional part = interpolation coefficient
    // linear interpolation
    A = d*A + (1.0-d)*prevA;
        
  }
    
  return A;  
}
 

/* 
 
 input = complex number
 output 
  - color array as Formal parameters
  - int = return code
*/
int compute_color(complex double c, unsigned char color[3]){
 
  double arg;
  unsigned char b;
   
  // compute 
  arg = Give_Arg( c, IterationMax);
     
  // 
  if (arg < 0.0)
    { /*  interior of Mandelbrot set = inside_color =  */
      color[0]=191; // 
      color[1]=191;
      color[2]=191;                           
    }
  else // exterior of Mandelbrot set = CPM  
     { 
       // gray gradient
      b = (unsigned char) (255 - 255*arg );
     
      color[0]= b;  /* Red*/
      color[1]= b;  /* Green */ 
      color[2]= b;  /* Blue */
    };
    
  return 0;
}
 
 
 
void setup(){
 
  //
  PixelWidth=(CxMax-CxMin)/iXmax;
  PixelHeight=(CyMax-CyMin)/iYmax;
  
  lnER = log(EscapeRadius); // ln(ER) 
        
         
  /*create new file,give it a name and open it in binary mode  */
  fp= fopen(filename,"wb"); /* b -  binary mode */
  /*write ASCII header to the file*/
  fprintf(fp,"P6\n %s\n %d\n %d\n %d\n",comment,iXmax,iYmax,MaxColorComponentValue);
 
}
 



void info(){

  double distortion;
  // widt/height
  double PixelsAspectRatio = (double)iXmax/iYmax;  // https://en.wikipedia.org/wiki/Aspect_ratio_(image) 
  double WorldAspectRatio = (CxMax-CxMin)/(CyMax-CyMin);
  printf("PixelsAspectRatio = %.16f \n", PixelsAspectRatio );
  printf("WorldAspectRatio = %.16f \n", WorldAspectRatio );
  distortion = PixelsAspectRatio - WorldAspectRatio;
  printf("distortion = %.16f ( it should be zero !)\n", distortion );
  printf("\n");
  printf("bailut value = Escape Radius = %.0f \n", EscapeRadius);
  printf("IterationMax = %d \n", IterationMax);
  printf("i_skip = %d = number of skipped elements ( including t0 )= %d \n", i_skip, i_skip+1);

  // file  
  printf("file %s saved. It is called .... in  A Cheritat wiki\n", filename);

}
 



 
void close(){
 
  fclose(fp);
  info(); 

}
 
 
 
 
 
// ************************************* main ************************* 
int main()
{
        
  complex double c;
        
        
 
  setup();      
        
        
  printf(" render = compute and write image data bytes to the file \n");
 
  for(iY=0;iY<iYmax;iY++)
    for(iX=0;iX<iXmax;iX++)
      { // compute pixel coordinate        
	c = give_c(iX, iY);  
	/* compute  pixel color (24 bit = 3 bytes) */
	compute_color(c,color);         
	/*write color to the file*/
	fwrite(color,1,3,fp);
      }
        
  
  
  close();
  
         
  return 0;
}

Ultra Fractal

uf
jh.ucl

Stripes {
; Jussi Härkönen, 07-01-02
; See Fibers and Things by Ron Barnett for a similar coloring
; for convergent fractals.
;
; TIA originally developed by Kerry Mitchell
; Curvature originally developed by Damien Jones
;
; 07-02-15 Added the interpolation mode parameter
; 07-03-07 Added the "None" interpolation mode
; 07-03-30 Added Skip iterations, seed (TIA only) and Mandel version (TIA only)
;          parameters
;
; 16-10-30 Stripes + Perlin routines added by jam.
;
; from jh.ucl
; http://formulas.ultrafractal.com/cgi/formuladb?view;file=jh.ucl;type=.txt

global:
  ; Precalculate the attenuation factor for attenuated sum
  float attenuationFactor = 1 - exp(-@attenuation)
  float lp = log(log(@bailout))
  float twopi = 2 * #pi


if (iteration > @skipIter)
  if (@coloring == "Stripes")
    ; Angle dependent component  in [-1,1]
    float TkAng = sin(@angularFrequency * atan2(z) + @angularOffset)

    ; Radius-dependent component in [-1,1]
    float TkRad = sin(@circularFrequency * log((cabs(z))) + @circularOffset)

    ; The weighted sum of radius and angle dependent terms is in [-1,1]
    ; Scale and translate it to [0,1]
    Tk = 0.5 + 0.5 * ((1 - @circularWeight) * TkAng + @circularWeight * TkRad)

  elseif (@coloring == "Curvature")
    ; Skip two first iterations in order to get two nonzero previous values

    if (iteration > @skipIter + 2)
      q = (z - zPrev) / (zPrev - zPrev2)
      Tk = (abs(atan2(q)) / (pi))
    else
      Tk = 0
    endif
  else ; "TIA"
    ; Skip first iteration
    if (iteration > 1 + @skipIter)
      ; |z(k-1)^n|
      if (@usePixel)
        ; Use pixel value
        znAbs = cabs(z - #pixel)
      else
        ; Use seed specified by the user
        znAbs = cabs(z - @seed)
      endif

      ; Bounds:
      ; m(k) = ||z(k-1)|^n - |c||
      float lowBound = abs(znAbs - cAbs)
      ; M(k) = |z(k-1)|^n + |c|
      float upBound = znAbs + cAbs

      ; T(k)
      Tk = (cabs(z) - lowBound) / (upBound - lowBound)
    endif
  endif

  ; Update previous sums
  sum3 = sum2
  sum2 = sum1
  sum1 = sum


  ; Calculate S(k)
  if (@attenuate)
    sum = attenuationFactor*sum + Tk
  else
    sum = sum + Tk
  endif

endif

; ...
;
#index = tempcolor

default:
  title = "Stripes"

  param coloring
    caption = "Coloring mode"
    default = 0
    enum = "Stripes" "TIA" "Curvature"
    hint = "Specifies the coloring algorithm."
  endparam

  param usePixel
    caption = "Mandel version"
    default = false
    visible = @coloring == "TIA"
    hint = "Specifies whether to use the algorithm adapted \
            to Julia or Mandelbrot sets"
  endparam

  param seed
    caption = "Seed"
    default = (0.5,0.25)
    visible = (@coloring == "TIA") && (@usePixel == false)
    hint = "Seed of the Julia set."
  endparam

  param interpolation
    caption = "Interpolation"
    default = 0
    enum = "Linear" "Smooth" "None"
    hint = "Specifies the interpolation method."
  endparam

  float param circularWeight
    caption = "Circle weight"
    default = 0.0
    visible = (@coloring == "Stripes")
    min = 0
    max = 1
    hint = "Weight of circular and radial components. Weight = 0 shows only \
            stripes, whereas weight = 1 shows only circles."
  endparam

  float param circularFrequency
    caption = "Circle frequency"
    default = 5.0
    visible = (@coloring == "Stripes")
    hint = "Specifies the density of circles. Use integer values for smooth results."
  endparam

  float param circularOffset
    caption = "Circle offset"
    default = 0.0
    visible = (@coloring == "Stripes")
    hint = "Circle offset given in radians."
  endparam

  float param angularFrequency
    caption = "Stripe density"
    default = 5.0
    visible = (@coloring == "Stripes")
    hint = "Specifies the density of stripes. Use integer values for smooth results."
  endparam

  float param angularOffset
    caption = "Stripe offset"
    default = 0.0
    visible = (@coloring == "Stripes")
    hint = "Stripe offset given in radians."
  endparam

  int param skipIter
    caption = "Skip iterations"
    default = 0
    hint = "Excludes a given number of iterations from the end of the orbit."
  endparam

  ; Fractional iteration count parameters
  float param bailout
    caption = "Bailout"
    default = 1e20
    min = 1
    hint = "Bail-out value used by the fractal formula. Use large bailouts for \
      best results."
  endparam

 ; Attenuation parameters
  bool param attenuate
    caption = "Moving average"
    default = false
    hint = "Use moving average instead of average. This may make the coloring \
      structure appear clearer (especially with TIA), but may also make busy \
      areas to appear flat."
  endparam

  float param attenuation
    caption = "Filter factor"
    default = 2
    visible = @attenuate
    min = 0
    hint = "Specifies the moving average strength. Values close to 0 cause the last \
      terms to be weighted creating results similar to Cilia. Large values \
      make the sum look more like the usual average."
  endparam
}

Fragmentarium

// http://www.fractalforums.com/general-discussion/stripe-average-coloring/
// code by Syntopia
#include "2D.frag"
#group Mandelbrot

// Number of iterations
uniform int Iterations; slider[10,200,5000]
uniform float R; slider[0,0,1]
uniform float G; slider[0,0.4,1]
uniform float B; slider[0,0.7,1]
uniform bool Julia; checkbox[false]
uniform float JuliaX; slider[-2,-0.6,2]
uniform float JuliaY; slider[-2,1.3,2]
vec2 c2 = vec2(JuliaX,JuliaY);

void init() {}

vec2 complexMul(vec2 a, vec2 b) {	
	return vec2( a.x*b.x - a.y*b.y,a.x*b.y + a.y * b.x);	
}

vec2 mapCenter = vec2(0.5,0.5);
float mapRadius =0.4;
uniform bool ShowMap; checkbox[true]
uniform float MapZoom; slider[0.01,2.1,6]

vec3 getMapColor2D(vec2 c) {
	vec2 p = (aaCoord-mapCenter)/(mapRadius);
	p*=MapZoom; p.x/=pixelSize.x/pixelSize.y;
	if (abs(p.x)<2.0*pixelSize.y*MapZoom) return vec3(0.0,0.0,0.0);
	if (abs(p.y)<2.0*pixelSize.x*MapZoom) return vec3(0.0,0.0,0.0);
	p +=vec2(JuliaX, JuliaY) ;
	vec2 z = vec2(0.0,0.0);
	int i = 0;
	for (i = 0; i < Iterations; i++) {	
		z = complexMul(z,z) +p;	
		if (dot(z,z)> 200.0) break;	
	}
	
	if (i < Iterations) {	
		float co = float( i) + 1.0 - log2(.5*log2(dot(z,z)));	
		co = sqrt(co/256.0);	
		return vec3( .5+.5*cos(6.2831*co),.5+.5*cos(6.2831*co),.5+.5*cos(6.2831*co) );	
	} else {
		return vec3(0.0);	
	}
	
}

// Skip initial iterations in coloring
uniform int Skip; slider[0,1,100]
// Scale color function
uniform float Multiplier; slider[-10,0,10]
uniform float StripeDensity; slider[-10,1,10]
// To test continous coloring
uniform float Test; slider[0,1,1]
uniform float EscapeRadius2; slider[0,1000,100000]

vec3 getColor2D(vec2 c) {
	if (ShowMap && Julia) {	
		vec2 w = (aaCoord-mapCenter);	
		w.y/=(pixelSize.y/pixelSize.x);	
		if (length(w)<mapRadius) return getMapColor2D(c);	
		if (length(w)<mapRadius+0.01) return vec3(0.0,0.0,0.0);	
	}
	
	vec2 z = Julia ? c : vec2(0.0,0.0);
	int i = 0;
	float count = 0.0;
	float avg = 0.0; // our average
	float lastAdded = 0.0;
	float z2 = 0.0; // last squared length
	for ( i = 0; i < Iterations; i++) {	
		z = complexMul(z,z) + (Julia ? c2 : c);	
		if (i>=Skip) {		
			count++;	
			lastAdded = 0.5+0.5*sin(StripeDensity*atan(z.y,z.x));		
			avg +=  lastAdded;
		}
		z2 = dot(z,z);
		if (z2> EscapeRadius2 && i>Skip) break;
	}
	float prevAvg = (avg -lastAdded)/(count-1.0);
	avg = avg/count;
	float frac =1.0+(log2(log(EscapeRadius2)/log(z2)));	
	frac*=Test;
	float mix = frac*avg+(1.0-frac)*prevAvg;
	if (i < Iterations) {		
		float co = mix*pow(10.0,Multiplier);
		co = clamp(co,0.0,10000.0);
		return vec3( .5+.5*cos(6.2831*co+R),.5+.5*cos(6.2831*co + G),.5+.5*cos(6.2831*co +B) );		
	} else {		
		return vec3(0.0);		
	}	
}

#preset Default
Center = -0.587525,0.297888
Zoom = 1.79585
AntiAliasScale = 1
AntiAlias = 1
Iterations = 278
R = 0
G = 0.4
B = 0.7
Julia = false
JuliaX = -0.6
JuliaY = 1.3
ShowMap = true
MapZoom = 2.1
Skip = 6
Multiplier = -0.1098
StripeDensity = 1.5384
Test = 1
EscapeRadius2 = 74468
#endpreset

#preset Julia
Center = -0.302544,-0.043626
Zoom = 4.45019
AntiAliasScale = 1
Iterations = 464
R = 0.58824
G = 0.3728
B = 0.27737
Julia = true
JuliaX = -1.26472
JuliaY = -0.05884
ShowMap = false
MapZoom = 1.74267
Skip = 4
Test = 1
EscapeRadius2 = 91489
Multiplier = 0.4424
StripeDensity = 2.5
AntiAlias = 2
#endpreset

Processing

// http://www.fractalforums.com/general-discussion/stripe-average-coloring/
// code by asimes{{typo help inline|reason=similar to amises|date=October 2022}}
// result : https://commons.wikimedia.org/wiki/File:Mandelbrot_set_-_Stripe_Average_Coloring.png

float xmin = -2.5; 
float ymin = -2.0; 
float wh = 4;
float stripes = 5.0;
int maxIterations = 1000;

void setup() {
  size(10000, 10000, P2D);
}

void draw() {
  loadPixels();
  float xmax = xmin+wh;
  float ymax = ymin+wh;
  float dx = (xmax-xmin)/width;
  float dy = (ymax-ymin)/height;
  float x = xmin;
  for (int i = 0; i < width; i++) {
    float y = ymin;
    for (int j = 0;  j < height; j++) {
      float zr = x;
      float zi = y;
      float lastzr = x;
      float lastzi = y;
      float orbitCount = 0;
      int n = 0;
      while (n < maxIterations) {
        float zrr = zr*zr;
        float zii = zi*zi;
        float twori = 2*zr*zi;
        zr = zrr-zii+x;
        zi = twori+y;
        if (zrr+zii > 10000) break;
        orbitCount += 0.5+0.5*sin(stripes*atan2(zi, zr));
        lastzr = zr;
        lastzi = zi;
        n++;
      }
      if (n == maxIterations) pixels[i+j*width] = 0;
      else {
        float lastOrbit = 0.5+0.5*sin(stripes*atan2(lastzi, lastzr));
        float smallCount = orbitCount-lastOrbit;
        orbitCount /= n;
        smallCount /= n-1;
        float frac = -1+log(2.0*log(10000))/log(2)-log(0.5*log(lastzr*lastzr+lastzi*lastzi))/log(2);      

        float mix = frac*orbitCount+(1-frac)*smallCount;
        float orbitColor = mix*255;
        pixels[i+j*width] = color(orbitColor);
      }
      y += dy;
    }
    x += dx;
  }
  updatePixels();
  noLoop();
}

void mousePressed(){
  save("a10000.png");
}

程式

參考文獻

[1] Jussi Harkonen : 關於分形著色技術

[2] A Cheritat wiki : 曼德爾布羅集 - 徑向線

[3] ractalforums : 三角不等式平均著色

[4] ractalforum : 關於分形著色技術的碩士論文

[5] ractalforums : 條紋平均著色

[6] 維基百科 : 插值

[7] 維基百科 : 小數部分

[8] stackoverflow 問題：如何使用 M_LN2 從 math.h？

[9] M_LN2

[1]

[2]

[3]

[4]

[5]

[8]

[9]