sandbox/Antoonvh/koch.c

The Fractal Dimension of the Koch Snowflake

This page presents an investigation on the Koch snowflake. You may want to watch this youtube video for a detailed, 20 minute introduction to the topic of fractal curves. Since we have the tools, let us use an adaptive quad-tree grid and the ppm2gif converter.

#include "grid/quadtree.h"
#include "utils.h"

int main(){

The Koch Curve

First, we need to obtain the snowflake curve. We determine it upto 8 fractal levels of refinement, resulting in 348200.000 points on the curve. This determines the maximum level of refinement that we can use in our grid based analysis.

  int jmax=9;
  int n=4;
  int nm=((n-1)*pow(4,jmax-1))+1;
  double xn[nm],yn[nm];
  double xk[nm];
  double yk[nm];
  xk[0]=0;
  xk[1]=1;
  xk[2]=0.5;
  xk[3]=0;
  yk[0]=0;
  yk[1]=0;
  yk[2]=-sqrt(3.)/2;
  yk[3]=0;
  fprintf(ferr,"\nyk[n-1]=%g\n\n",yk[n-1]);
  for (int j=1;j<jmax;j++){//For every fractal level
    for (int k=0;k<(n-1);k++){//For every existing point add three new points. 
      xn[(k*4)]=xk[(k)];
      xn[(k*4)+1]=(2*xk[k]+xk[k+1])/3;
      xn[(k*4)+2]=((xk[k]+xk[k+1])/2)+((yk[k]-yk[k+1])*sqrt(3.)/6);      
      xn[(k*4)+3]=(xk[k]+2*xk[k+1])/3;
      yn[(k*4)]=yk[k];
      yn[(k*4)+1]=(2*yk[k]+yk[k+1])/3;
      yn[(k*4)+2]=((yk[k]+yk[k+1])/2)-((xk[k]-xk[k+1])*sqrt(3.)/6);
      yn[(k*4)+3]=(yk[k]+2*yk[k+1])/3;
    }
    yn[4*(n-1)]=yk[n-1];
    xn[4*(n-1)]=xk[n-1];
    FILE * fpkoch = fopen("koch.dat","w");
    for (int l=0;l<(4*(n-1)+1);l++){// Copy yn and xn into yk and xk. 
      yk[l]=yn[l];
      xk[l]=xn[l];
    }
    if (j==jmax-1){//Output the points of the curve after the last interation. 
      for (int l=0;l<(4*(n-1)+1);l++){
        fprintf(fpkoch,"%g %g\n",xk[l],yk[l]);
      }
    }
    n=(n-1)*4+1;
  }

Let us check-out the obtained snowflake.

Looks O.K. Note that the curve is defined at a much finer resolution than is displayed in the plot.

The Algorithm

To find the fractal dimension of this curve we iteratively refine the grid and log how many cells are located on the curve.

  FILE * fpit = fopen("cells.dat","w");
  X0=-0.5;
  Y0=-1.2;
  L0=2.0;
  init_grid(32);
  
  scalar c[];
  for (int m=0;m<7;m++){
    foreach()
      c[]=0;
    for(int i=0;i<n;i+=64/pow(2,m)){
      Point point = locate(xn[i],yn[i]);
      c[]+=1.0;
    }
    int gcc=0;
    int gcn=0;
    foreach(reduction(+:gcc) reduction(+:gcn)){
      gcn++;
      if(c[]>0.5)
	gcc++;
    }

For the refinment we (once again) have faith in the wavelet based adaption algorithm to refine (and coarsen) the grid.

    boundary({c});
    while(adapt_wavelet({c},(double[]){0.1},(7+m),4,{c}).nf);
    
    fprintf(fpit,"%d %d %d\n",m,gcc,gcn);
    fflush(fpit);
    int i=m;
    scalar lev[];
    foreach(){
      lev[]=level;
    }
    static FILE * fp1 = popen ("ppm2gif --delay 200 > f.gif", "w");
    output_ppm (lev, fp1,512,min=4,max=15);
  }
}

Results

The adaptive grid that is used to refine the cells in the neighborhood of the curve is visualized below. More red colors represent a higher level of refinement.

The various stages of refinement

The various stages of refinement

Again the displayed resolution (512×512) is much less than the final resulution of our analysis (that corresponds to 4096×4096).

Now we check if we can find the fractal dimension of this curve by plotting the number of grid cells against the refinement iteration.

Apparently the fractal dimension is 1.26 and not 1. The latter value would be typical for a ‘regular’ curve. As was explained in the aforementioned youtube video and also noted in another source, the fractal dimension can be expressed analytically as: ln(4)ln(3)1.26. That number is rather close to the one found with the present method.

Finally, we check if the total employed number of gridcells scales with the fractal dimension of this “problem”.

It does appear to be the case. Well done adapt_wavelet() function!

Follow-up

There are two other, related cases to adress the following questions: