简介

"平均颜色是使用平滑迭代计数的小数部分在平均总和之间进行插值的着色函数族。" ^[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 [email protected]: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://wikibooks.cn/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]