1
$\begingroup$

hey, please check the edits at the bottom of the post, thanks :)

In these days I tried to write a program that computes the dimension of a fractal (given an image) in Java (more specifically in Processing 3).

I'm using a simple box-counting method (just scale the image and count the boxes the fractal "touches" every time).

In this super-mega-sketch I just :

  • Take a bunch of convolutions to get an image of the border of the fractal, using the Sobel operator here's the Wikipedia link

  • I take only the really, really "bright" pixels

  • Scale the image a bunch of times and store in the memory how many boxes are "touched" (I don't really count the number of boxes the fractal touches, I just scale the image and count the white pixels) and the scaling factor, by which you scale the image. I store each "box" and each scaling factor in arrays.

****Warning Critical Part****

  • take the natural logarithm of every item in the Array in order to compute a linear regression. The slope of the line should be the dimension of the fractal, since:

eq.1 - y = ax^(d)

//where y is the number of boxes touched, and x is the scaling factor

eq.2 - ln(y) = ln(a) + d*ln(x)

//If you where to plot the log of the number of boxes touched in function of the log of the scaling factor, you should get a linear relation y = a + dx (or y = ln(a) + dx , I don't really know cause I'm stupid :( ).

The linear regression doesn't really work tough, the resulting dimension of the Britain is always 0.99983 (of course it doesn't make sense, the dimension of a fractal can't be less than one [or can it? ow0]).

I also tried to compute the dimension of the fractal through the median slope of each point in the array, just use this formula (m = (y_a - y_b)/(x_a-x_b)). Of course I get nonsense again, for the Britain I get a dimension of 1.99983 (same as before + 1);

I guess I made some mistake taking the natural log of both arrays (check the sketch below). Since Desmos says ln(y) = ln(a) + d*ln(x) is quadratic... Check this

There might also be something wrong with my data normalization, though I don't thing so...

Though I spent my all day trying to figure what's wrong, and here I am... asking for help.

Please let me know if you find anything

Here's my sketch, thanks in advance!!!!!

//Variables//Edge Detection//Variables//
PImage img,imgGray,imgBlur,sobelX,sobelY,sobelMagnitude;
boolean start=false;
int matrixsize = 3;

final int BLURS = 0;
final color FILTER = color(192,192,192);

float[][] gaussianBlurKernel={
                             { 0.0625, 0.125, 0.0625 },
                             { 0.125, 0.25, 0.125 },
                             { 0.0625, 0.123, 0.0625 }
                             };
float[][]sobelOperatorX={
                        { -1, 0, 1 },
                        { -2, 0, 2 },
                        { -1, 0, 1 }
                        };
float[][]sobelOperatorY={
                        { -1, -2, -1 },
                        { 0, 0, 0 },
                        { 1, 2, 1 }
                        };

//Variables//Nearest Neightbour and Box Counting Algorythm//Variables//
final int DATASET_SIZE = 128;

PImage fractal;
float[] dataset_boxes = new float[DATASET_SIZE];
float[] dataset_scaling_factors = new float[DATASET_SIZE];

//Varibles//Computing the Fractal Dimension//Variables//
float m = 1;
float q = 0;

float[] log_dataset_scaling_factors = new float[DATASET_SIZE];
float[] log_dataset_boxes = new float[DATASET_SIZE];

float c; //e^q = c
float d; //y = c + d*x
float[] D = new float[DATASET_SIZE];

void setup(){
  //fullScreen(2);
  size(8192,8192);
  selectInput("Select the image of a fractal","imageAssigner");
  while(!start){
    println();
  }
  noSmooth();
}

void imageAssigner(File selection){
  String h = selection.getAbsolutePath();
  img = loadImage(selection.getAbsolutePath());
  if (selection == null) {
    println("Window was closed or the user hit cancel.");
  } else {
    println("User selected " + selection.getAbsolutePath());
  }
  start=true;
}

void draw(){
  //Start//Edge Detection//Start//

  img.filter(GRAY);
  img.save("1-imgGray.jpg");

  imgGray = loadImage("1-imgGray.jpg");

  int xstart = 0;
  int ystart = 0;
  int xend = imgGray.width;
  int yend = imgGray.height;

  for(int i = 0; i<BLURS; i++){
    for (int x = xstart; x < xend; x++) {
      for (int y = ystart; y < yend; y++ ) {
        color c = convolution(x, y, gaussianBlurKernel, matrixsize, imgGray);
        int loc = x +y*imgGray.width;
        imgGray.pixels[loc] = c;
      }
    }
  }

  imgGray.save("2-imgBlur.jpg");

  sobelX=loadImage("2-imgBlur.jpg");
  image(sobelX, 0, 0);

  loadPixels();
  for (int x = xstart; x < xend; x++) {
    for (int y = ystart; y < yend; y++ ) {
      color c = convolution(x, y, sobelOperatorX, matrixsize, sobelX);
      int loc = x + y*width;
      pixels[loc] = c;
    }
  }
  updatePixels();

  sobelX = get(0,0,sobelX.width,sobelX.height);
  sobelX.save("3-sobelX.jpg");

  background(0);
  sobelY=loadImage("2-imgBlur.jpg");
  image(sobelY, 0, 0);

  loadPixels();
  for (int x = xstart; x < xend; x++) {
    for (int y = ystart; y < yend; y++ ) {
      color c = convolution(x, y, sobelOperatorY, matrixsize, sobelY);
      int loc = x + y*width;
      pixels[loc] = c;
    }
  }
  updatePixels();

  sobelY = get(0,0,sobelY.width,sobelY.height);
  sobelY.save("4-sobelY.jpg");

  //background(0);
  sobelX = loadImage("3-sobelX.jpg");
  sobelY = loadImage("4-sobelY.jpg");
  sobelMagnitude = loadImage("3-sobelX.jpg");

  for (int x = xstart; x < xend; x++) {
    for (int y = ystart; y < yend; y++ ) {
      color c = 0;
      int loc = x + y*sobelMagnitude.width;
      sobelMagnitude.pixels[loc] = c;
    }
  }

  //background(0);

  loadPixels();
  for (int x = xstart; x < xend; x++) {
    for (int y = ystart; y < yend; y++ ) {
      int loc = x + y*sobelX.width;

      if(sobelX.pixels[loc]<=FILTER)//metà sarebbe -8388608
        sobelX.pixels[loc] = #000000;
      else
        sobelX.pixels[loc] = #FFFFFF;

      if(sobelY.pixels[loc]<=FILTER)//metà sarebbe -8388608
        sobelY.pixels[loc] = #000000;
      else
        sobelY.pixels[loc] = #FFFFFF;

      //pixels[loc] = Math.round(sqrt(pow(sobelX.pixels[loc],2) + pow(sobelY.pixels[loc],2)));
      //pixels[loc] = sobelX.pixels[loc] + sobelY.pixels[loc];

      sobelMagnitude.pixels[loc] = sobelX.pixels[loc] + sobelY.pixels[loc];
    }
  }
  updatePixels();

  sobelX.save("6-sobelX-debug.jpg");
  sobelY.save("7-sobelY-debug.jpg");
  sobelMagnitude.save("5-LivingOnTheEdge.jpg");  

  //End//Edge Detection//End//

  //Start//Scaling and Box Counting Algorythm//Start//

  fractal=loadImage("5-LivingOnTheEdge.jpg");


  for(int i=0; i<DATASET_SIZE; i++){
    Scale(i,fractal,DATASET_SIZE);
    dataset_boxes[i] = countBoxes();
    println("[",i,"] = ",countBoxes());
    //String string = "s"+str(i)+"-ScaledFractal.jpg";
    //save(string);
  }
  save("8-ScaledFractal.jpg");
  //End//Scaling//End//

  //Start//Computing the Fractal Dimension//Start//
  pushMatrix();
  rotate(PI/2);
  translate(-width/2,height/2);
  for(int i = 0; i < DATASET_SIZE; i++){
    log_dataset_scaling_factors[i] = log(dataset_scaling_factors[i]);
    log_dataset_boxes[i] = log(dataset_boxes[i]);
  }
  background(51);  
  Normalize();
  for(int i = 0; i < dataset_scaling_factors.length; i++){
    float x = map(dataset_scaling_factors[i], 0, 1, 0, width);
    float y = map(dataset_boxes[i], 1, 0, 0, height);
    fill(255);
    stroke(255);
    ellipse(x,y,8,8);
  }

  linearRegression(dataset_scaling_factors,dataset_boxes);
  drawLine();

  println("the dimension of the fractal is : ",m);

  for(int i = 1; i < D.length; i++){
    D[i] = (log_dataset_boxes[i] - log_dataset_boxes[i-1])/(log_dataset_scaling_factors[i] - log_dataset_scaling_factors[i-1]);
  }

  for(int i = 0; i < D.length; i++){
    d += D[i];
  }

  d /= D.length-1;

  println("maybe the dimension of the fractal is ",d);

  popMatrix();
  exit();
   //End//Computing The Fractal Dimension//End//
}

//functions//Edge Detection//functions//

color convolution(int x, int y, float[][] matrix, int matrixsize, PImage img){
  float rtotal = 0.0;
  float gtotal = 0.0;
  float btotal = 0.0;
  int offset = matrixsize / 2;
  for (int i = 0; i < matrixsize; i++){
    for (int j= 0; j < matrixsize; j++){
      int xloc = x+i-offset;
      int yloc = y+j-offset;
      int loc = xloc + img.width*yloc;
      loc = constrain(loc,0,img.pixels.length-1);
      rtotal += (red(img.pixels[loc]) * matrix[i][j]);
      gtotal += (green(img.pixels[loc]) * matrix[i][j]);
      btotal += (blue(img.pixels[loc]) * matrix[i][j]);
    }
  }
  rtotal = constrain(rtotal, 0, 255);
  gtotal = constrain(gtotal, 0, 255);
  btotal = constrain(btotal, 0, 255);
  return color(rtotal, gtotal, btotal);
}

//functions//Scaling and Box Counting Algorythm//functions//

void Scale(int iteration, PImage image,int dataset){
  float scalingFactor;

  if(width - image.width > height - image.height){
    float lengthDifference = height - image.height;
    float increment = iteration*lengthDifference/dataset;
    float supposedLenght = image.height + increment;
    scalingFactor = supposedLenght/image.height;
    println(supposedLenght/image.height);
  }
  else if(width - image.width < height - image.height){
    int lengthDifference = width - image.width;
    int increment = iteration*lengthDifference/dataset;
    int supposedLenght = image.width + increment;
    scalingFactor = supposedLenght/image.width;
  }
  else {
    scalingFactor = 0.9;
  }

  dataset_scaling_factors[iteration] = scalingFactor;
  pushMatrix();
  scale(scalingFactor);
  background(0);
  image(image,0,0);
  popMatrix();
};

int countBoxes(){
  loadPixels();
  int counter = 0;
  for(int i = 0; i < width; i++){
    for(int j = 0; j < height; j++){
      int loc = i + j*width;
      if(pixels[loc] == #FFFFFF)
        counter++;
    }
  }
  updatePixels();
  return counter;
}

//Functions//Computing the Fractal Dimension//Functions//
void drawLine(){
  float x1 = 0;
  float y1 = m * x1 + q;
  float x2 = 1;
  float y2 = m * x2 + q;

  x1 = map(x1, 0, 1, 0, width);
  y1 = map(y1, 0, 1, height, 0);
  x2 = map(x2, 0, 1, 0, width);
  y2 = map(y2, 0, 1, height, 0);

  stroke(255,0,255);
  line(x1,y1,x2,y2);
}

void linearRegression(float[] dataX, float[] dataY){
  float xsum = 0;
  float ysum = 0;
  for(int i = 0; i < dataX.length; i++){
    xsum += dataX[i];
    ysum += dataY[i];
  }
  float xmean = xsum / dataX.length;
  float ymean = ysum / dataX.length;

  float num = 0; //numerator
  float den = 0; //denominator

  for(int i = 0; i < dataX.length; i++){
    float x = dataX[i];
    float y = dataY[i];

    num += (x - xmean) * (y - ymean);
    den += (x - xmean) * (x - xmean);
  }
  m = num / den;
  q = ymean - m * xmean;
}

void Normalize(){
  float minX = dataset_scaling_factors[0];
  float maxX = dataset_scaling_factors[dataset_scaling_factors.length - 1];
  float minY = dataset_boxes[0];
  float maxY = dataset_boxes[dataset_boxes.length - 1];

  for(int i = 0; i < dataset_scaling_factors.length; i++){
    dataset_scaling_factors[i] = map(dataset_scaling_factors[i],minX,maxX,0,1);
  }
  for(int i = 0; i < dataset_boxes.length; i++){
    dataset_boxes[i] = map(dataset_boxes[i],minY,maxY,0,1);
  }
}

EDIT Yes... Normalizing is a very very bad idea, if you use the map function and don't use the same minimum and maximum in the first two parameters. It's worthless anyway... Remove it if you are trying to get my code work.

EDIT Ok, I checked and checked... the linear regression works properly, what doesn't work is the whole thing before... I don't thing scaling the image and counting the white pixels works as I thought...

So, I am completely open to suggestions :

  • which language would you use? Wolfram Mathematica? MATLAB? please let me know.

  • which method would you use to compute the dimension of a non-self-similar fractal?

$\endgroup$

closed as off-topic by Xander Henderson, Paul Frost, Christopher, ancientmathematician, Rafa Budría Sep 3 '18 at 19:12

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question is not about mathematics, within the scope defined in the help center." – Xander Henderson, Paul Frost, Christopher, ancientmathematician, Rafa Budría
If this question can be reworded to fit the rules in the help center, please edit the question.

Browse other questions tagged or ask your own question.