Fractals/Iterations in the complex plane/stripeAC

From Wikibooks, open books for an open world
Jump to navigation Jump to search

introduction[edit | edit source]

"Average colorings are a family of coloring functions that use the decimal part of the smooth iteration count to interpolate between average sums." [1]

Here stripes are lines perpendicular to the iteration bands ( level sets of escape time). These " radial strands ... follows the external rays which are of great importance in the theoretical study of M and in holomorphic dynamics." (A Cheritat)[2]


" is basically a cheap way to calculate an angle. TIA is averaging the angle over all iterations to get a smooth result. So there is some initialization and some calculations per iteration to do the sum. It has to store the sum at the previous iteration before adding the next one so you can interpolate between them to get a continuous function between iteration bands (in the same way that continuous potential is usually interpolated). This interpolation is done after the iterations bail out."[3] (xenodreambuie - author of Jux)

" ... it probably shows curves close to external rays, but nevertheless you do not have an explicit relation between the number returned by SAC and the true external argument." Wolf Jung

name[edit | edit source]

  • SAM = Stripe Average Method
  • SAC = Stripe Average Coloring[4]
  • radial strands

history[edit | edit source]

A new coloring called the Stripe Average was introduced by Jussi Härkönen in 2007. It is based on the behavior of the Triangle Inequality Average coloring.

examples[edit | edit source]

Images[edit | edit source]

fractalforums stripe-average-mandelbrot-2 by Zephitmaal

video[edit | edit source]


algorithm[edit | edit source]

For every point c of rectangle from parameter plane do:[5]

  • Calculate the orbit of the point z0. The result is series of complex numbers {z0,z1,z2,...}
  • Apply function t to this series to get a new series of real numbers: { t0=t(z0), t1=t(z1), t2=t(z2), ...}.
  • Take the average A ( = arithmetic mean) of this series of real numbers
  • map real number A to a color ( coloring function)

orbit[edit | edit source]

Truncated orbit:


Skipped orbit is a subset of truncated orbit:

  • skip a first k+1 of items from the sequence
  • do not use some first k+1 items to compute average


addend function[edit | edit source]

The addend function t

  • maps complex number to real number
  • is continuous on the elements of elements in skipped orbit

where :

  • s is a stripe density

average[edit | edit source]

full[edit | edit source]

Average of full truncated orbit:

  • compute average from all points of the truncated orbit



If one compute average of all points ( from i=1 to i=n) then:

  • stripes will be visible ( good)
  • level sets of escape time will be also visible (bad, discontinuity )

The solution is :

  • the interpolation
  • skipped truncated orbit

skipped[edit | edit source]

Average of skipped orbit:

  • skip first (k+1) items from computing average


note that:


It removes some discontinuities but not level sets of escape time

interpolated[edit | edit source]

Level sets of escape time are still visible ( discontinuity). One can remove them using interpolation.[6]

interpolation:

  • between 2 values
  • using:
    • linear function = linear interpolation
    • splines = spline interpolation
    • polynomial


steps:

  • Construct a smooth iteration count u = real number ( see equation 3.11 from page 21 in J Harkonen thesis)
  • take it's fractional part[7] d


Decimal (frational) part d is:

  • zero at boundary of previous level set
  • converges to 1 in the neighborhood of next level set boundary

"Thus it can be used as the interpolation coefficient."

linear[edit | edit source]

To compute linear average interpolate between two values:

  • one obtained from the "full" series ( truncated and skipped) = Avg =
  • one obtained from the average of all t elements except the last one = prevAvg =


 

see equation 4.2 from page 28 in J Harkonen thesis

or in GLSL code by Syntopia:

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

Here is important fragment of GLSL code by Syntopia:

// 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;
spline[edit | edit source]

Compute spline average using Catmull-Rom splines:

Polynomials :


then


see equation 4.3 from page 28 in J Harkonen thesis

Code[edit | edit source]

C[edit | edit source]

  • M_LN2 is defined in 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://en.wikibooks.org/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[edit | edit source]

Ultra Fractal

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[edit | edit source]

// 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[edit | edit source]

// 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");
}

Programs[edit | edit source]

References[edit | edit source]

  1. Jussi Harkonen : On Fractal Coloring Techniques
  2. A Cheritat wiki : Mandelbrot_set - Radial_strands
  3. fractalforums : triangle-inequality-average-coloring
  4. fractalforum : master's-thesis-on-fractal-coloring-techniques
  5. fractalforums : stripe-average-coloring
  6. wikipedia : Interpolation
  7. wijkipedia : Fractional part
  8. stackoverflow question: how-to-use-m-ln2-from-math-h ?
  9. M_LN2