Drawing Koch's curve

Greetings,

I’ve been trying to draw a Koch’s curve for some time now and this is the best result so far, but i still have problems with calculating the third point of a triangle with given 2.
If I set number of iterations to 1,2 or 3 the program draws it correctly but if iteration count rises on some places program fails to calculate the prober third point. As it can bee seen if the code is compiled and run on some places the third point is set on the inner part not on the outer and it looks like the distance is miscalculated as well.

I really couldn’t understand why is this happening, hope some of you can.

Regards.


#include <GL/glut.h>
#include <math.h>
#include <list>
#include <iostream>
#define PI 3.14159/180

typedef struct _point {
  float x;
  float y;
}point;

using namespace std;

/* Deklaracije callback funkcija. */
static void on_display(void);

int main(int argc, char **argv)
{
    /* Inicijalizuje se GLUT. */
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_RGB | /*GLUT_DEPTH | */ GLUT_DOUBLE);

    /* Kreira se prozor. */
    glutInitWindowSize(800, 800);
    glutInitWindowPosition(100, 100);
    glutCreateWindow(argv[0]);
    glEnable (GL_LINE_SMOOTH);
    glEnable (GL_BLEND);
    glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
    glHint (GL_LINE_SMOOTH_HINT, GL_NICEST);


    /* Registruju se callback funkcije. */
    glutDisplayFunc(on_display);

    /* Obavlja se OpenGL inicijalizacija. */
    glClearColor(0.75, 0.75, 0.75, 0);
    glEnable(GL_DEPTH_TEST);

    /* Program ulazi u glavnu petlju. */
    glutMainLoop();

    return 0;
}

static void on_display(void)
{
    /* Brise se prethodni sadrzaj prozora. */
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  
    list<point> l;
    list<point>::iterator i;
    point niz[3];
    point p1,p2;
    
    p1.x = -0.8;
    p1.y = 0;
    p2.x = 0.8;
    p2.y = 0;
    
    l.push_back(p1);
    l.push_back(p2);

    int br_iteracija = 4;

    while(br_iteracija > 0)
    {
        unsigned z = l.size();
        int p = 1;
        int j = 0;
        while(p < z)
        {
            if(j==0)
            {
              i = l.begin();
	      point pom1 = *i;
              i++;
	      point pom3 = *i;
              niz[0].x = ((pom3.x-pom1.x)/3) + pom1.x;
              niz[0].y = ((pom3.y-pom1.y)/3) + pom1.y;
	      niz[2].x = (2*((pom3.x-pom1.x))/3) + pom1.x;
	      niz[2].y = (2*((pom3.y-pom1.y))/3) + pom1.y;
	      
	      float r = sqrt(pow(niz[0].x-niz[2].x,2)+pow(niz[0].y-niz[2].y,2));
	      float altitude = (r*sqrt(3.0))/2;
	      float dx = niz[0].x - niz[2].x;
	      
	      point p;
	      p.x = ((pom3.x-pom1.x)/2) + pom1.x;
	      p.y = ((pom3.y-pom1.y)/2) + pom1.y;
	      
	      float slope = -dx / (niz[0].y-niz[2].y+0.000001);
	      
	      niz[1].x = sqrt(fabs(altitude*altitude - dx*dx)) / slope + p.x;
	      niz[1].y = slope * (niz[1].x - p.x) + p.y;
	      
              l.insert(i,niz,niz+3);
            }
            else
            {
                point pom1 = *i;
                i++;
		point pom3 = *i;
                niz[0].x = ((pom3.x-pom1.x)/3) + pom1.x;
		niz[0].y = ((pom3.y-pom1.y)/3) + pom1.y;
		niz[2].x = (2*((pom3.x-pom1.x))/3) + pom1.x;
		niz[2].y = (2*((pom3.y-pom1.y))/3) + pom1.y;
		
		float r = sqrt(pow(niz[0].x-niz[2].x,2)+pow(niz[0].y-niz[2].y,2));
		float altitude = (r*sqrt(3.0))/2;
		float dx = niz[0].x - niz[2].x;
		
		point p;
		p.x = ((pom3.x-pom1.x)/2) + pom1.x;
		p.y = ((pom3.y-pom1.y)/2) + pom1.y;
		
		float slope = -dx / (niz[0].y-niz[2].y+0.000001);
		
		niz[1].x = sqrt(fabs(altitude*altitude - dx*dx)) / slope + p.x;
		niz[1].y = slope * (niz[1].x - p.x) + p.y;
		
                l.insert(i,niz,niz+3);
            }
            p++;
            j++;
        }
        br_iteracija--;
    }

    glColor3f(0,0,1);
    //glPointSize(2);
    list<point>::iterator it;
    it = l.end();
    it--;
    
      for(i = l.begin();i!=it;i++) {
	glBegin(GL_LINES);
	point tmp = *i;
	i++;
	point tmp2 = *i;
	glVertex2f(tmp.x,tmp.y);
	glVertex2f(tmp2.x,tmp2.y);
	glEnd();
	i--;
      }

    /* Nova slika se salje na ekran. */
    glutSwapBuffers();
}


tutorial on doing the koch curve on the gpu.
It should be pretty simple to port to your code

I dont see much similarity, the only thing i’m failing is calculating third point

Ok here is the comented code:


#include <GL/glut.h>
#include <math.h>
#include <list>

/* Structure that is used to hold (x,y) coordinate of a Von Koch's curve */
typedef struct _point {
  float x;
  float y;
}point;

using namespace std;

/* Declaration of callback function */
static void on_display(void);

int main(int argc, char **argv)
{
    /* GLUT initialization */
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);

    /* Window creation */
    glutInitWindowSize(800, 800);
    glutInitWindowPosition(100, 100);
    glutCreateWindow("Von Koch's curve");

    /* Enabling blending and antialiasing for lines */
    glEnable (GL_LINE_SMOOTH);
    glEnable (GL_BLEND);
    glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
    glHint (GL_LINE_SMOOTH_HINT, GL_NICEST);

    /* Registration of callback function */
    glutDisplayFunc(on_display);

    /* OpenGL initialization, bacground color set to gray */
    glClearColor(0.75, 0.75, 0.75, 0);

    /* Program enters the main loop */
    glutMainLoop();

    return 0;
}

static void on_display(void)
{
    /* Clearing the previous window content */
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    /* pointList will contain points of Koch's curve at the end of iterations */
    list<point> pointList;
    list<point>::iterator i;

    /* interpolationArry will contain for every 2 points calculated 3 new points which represent a triangle */
    point interpolationArray[3];
    /* The program will start from 2 initial points and then for every 2 points it will interpolate 3 new point */
    point initialPoint1,initialPoint2;

    /* Setting the coordinates */
    initialPoint1.x = -0.8;
    initialPoint1.y = 0;
    initialPoint2.x = 0.8;
    initialPoint2.y = 0;

    /* Inserting initial points */
    pointList.push_back(initialPoint1);
    pointList.push_back(initialPoint2);

    /* Here the problem starts, for iteratioCount 1,2 the program works proerly but if it is switched to 3 or 4 ...
       it begins to miscalculate the third point of triangle */
    int iterationCount = 3;

    /* For number of iterations, for each iteration and every 2 points we interpolate 3 new points
      for example we have points 1 2 -> 1 3 4 5 2 -> 1 6 7 8 3 9 10 11 4 12 13 14 5 15 16 17 2 -> ... */
    while(iterationCount > 0)
    {
	/* Getting the current list size */
        unsigned pointListSize = pointList.size();
        int counter = 1;
        int j = 0;
	/* If in current iteration list conatins n points it is needed to interpolate n-1 times */
        while(counter < pointListSize)
        {
	    /* If it is a first iteration of interpolation interator must be set to point at second element */
            if(j==0)
            {
              i = pointList.begin();
	      /* Get the first point */
	      point pom1 = *i;
              i++;
	      /* Get the second point */
	      point pom3 = *i;
	      /* For the given 2 points calculate 3 new points
		 first of 3 new points is 1/3 of the length of defined line given by 2 points we took */
              interpolationArray[0].x = ((pom3.x-pom1.x)/3) + pom1.x;
              interpolationArray[0].y = ((pom3.y-pom1.y)/3) + pom1.y;
	      /* Second one is 2/3 of length */
	      interpolationArray[2].x = (2*((pom3.x-pom1.x))/3) + pom1.x;
	      interpolationArray[2].y = (2*((pom3.y-pom1.y))/3) + pom1.y;

	      /* Now it is needed to find the third point which will define a triangle (isoscale) */
	      /* Calculate the side length */
	      // THE PROBLEM IS PROBABLY HERE SOMETHING IS MISCALCULATED, AND IN ITERATION COUNT > 3
	      // IT STARTS TO MISCALCULATE THIRD POINT
	      float r = sqrt(pow(interpolationArray[0].x-interpolationArray[2].x,2)+pow(interpolationArray[0].y-interpolationArray[2].y,2));
	      /* Calculate the altitude */
	      float altitude = (r*sqrt(3.0))/2;
	      /* Calculate dx */
	      float dx = interpolationArray[0].x - interpolationArray[2].x;

	      /* Calculate midPoint from the previously calculated 2 */
	      point midPoint;
	      midPoint.x = ((pom3.x-pom1.x)/2) + pom1.x;
	      midPoint.y = ((pom3.y-pom1.y)/2) + pom1.y;

	      /* Calculate slope, +0,000001 to avoid deviding by 0 in some cases */
	      float slope = -dx / (interpolationArray[0].y-interpolationArray[2].y+0.000001);

	      /* Finaly obtain the third point */
	      interpolationArray[1].x = sqrt(fabs(altitude*altitude - dx*dx)) / slope + midPoint.x;
	      interpolationArray[1].y = slope * (interpolationArray[1].x - midPoint.x) + midPoint.y;

	      /* Interpolate 3 new points between current 2 */
              pointList.insert(i,interpolationArray,interpolationArray+3);
            }
            /* Otherwise iterator is set on proper location from the previous iteration */
            else
            {
		/* This is the same procedure just a special case when the iteratior is set on
		   proper location from the previous iteration */
                point pom1 = *i;
                i++;
		point pom3 = *i;
                interpolationArray[0].x = ((pom3.x-pom1.x)/3) + pom1.x;
		interpolationArray[0].y = ((pom3.y-pom1.y)/3) + pom1.y;
		interpolationArray[2].x = (2*((pom3.x-pom1.x))/3) + pom1.x;
		interpolationArray[2].y = (2*((pom3.y-pom1.y))/3) + pom1.y;
		
		float r = sqrt(pow(interpolationArray[0].x-interpolationArray[2].x,2)+pow(interpolationArray[0].y-interpolationArray[2].y,2));
		float altitude = (r*sqrt(3.0))/2;
		float dx = interpolationArray[0].x - interpolationArray[2].x;
		
		point midPoint;
		midPoint.x = ((pom3.x-pom1.x)/2) + pom1.x;
		midPoint.y = ((pom3.y-pom1.y)/2) + pom1.y;
		
		float slope = -dx / (interpolationArray[0].y-interpolationArray[2].y+0.000001);
		
		interpolationArray[1].x = sqrt(fabs(altitude*altitude - dx*dx)) / slope + midPoint.x;
		interpolationArray[1].y = slope * (interpolationArray[1].x - midPoint.x) + midPoint.y;
		
                pointList.insert(i,interpolationArray,interpolationArray+3);
            }
            counter++;
            j++;
        }
        iterationCount--;
    }

  
    glColor3f(0,0,1);
    list<point>::iterator it;
    it = pointList.end();
    it--;

    /* Interate through list and draw lines from i,i+1 point respectevly */
    for(i = pointList.begin();i!=it;i++) {
      glBegin(GL_LINES);
	point tmp = *i;
	i++;
	point tmp2 = *i;
	glVertex2f(tmp.x,tmp.y);
	glVertex2f(tmp2.x,tmp2.y);
      glEnd();
      i--;
    }

    /* Send new picture on display */
    glutSwapBuffers();
}


Help people, I’m really stuck here :frowning:

Sorry… if you could write up all the math, maybe with some diagrams, that would hielp a ton, but otherwise I’m not sure how to start looking at your code. (well, at least not quickly or without a headache)

There is a way with which you can draw generalized koch curves (base-motif fractals) MUCH easier than this. I’ll describe to you my way:

first, you have two things: A base, which is a series of line segments anywhere in the coordinate system. For any two adjacent line segments on the base, we call them A and B. The second thing is a motif, which is also a series of line segments, with its first point at (0,0), and its second point at (1,0). We will call a point on this motif, (x,y).
in the case of the koch curve:
base:

motif:

So, FOR every two adjacent points in the base (line AB), take the motif, stretch it, rotate it, and place it inbetween the two points. In other words, the interval on the X axis [0,1], is transformed to A and B.

Now, imagine if you will, a point on the motif being the value given by the addition of two parametric lines. One horizontal, from 0 to 1:
finalx=(ending_x-starting_x)*x+starting_x
finaly=0

and another vertical, from 0 to 1
finalx=0
finaly=(ending_y-starting_y)*y+starting_y

the solution, quite simply, is finalx=x and finaly=y, but the idea of a parametric line like this can help us.

If you have a point on the X axis, and you want to put it between A and B, you just use x as the coefficient in the parametric equation of a line between A and B:
finalx=(B.x-A.x)*x+A.x
finaly=(B.y-A.y)*x+A.y

In matrix form, where A and B are column vector matrices:
final=(B-A)*x+A

Now, that only works for a point where the y value is 0. If you wanted to include the y value, how could you write a parametric equation similar to the above?

well… the vector (B-A) is B’s position relative to A, so A is the origin, and B is out there somewhere. so, If you apply a rotation matrix to (B-A), you’re rotating B about A. if you rotate B 90 degrees to the left about A, and then use the parametric line from A to B’, then you can say, if x=0:
[0 -1] *(B-A)*y+A = final
[1 -0]

So, basically, transform your coordinates so that A is the origin to get (B-A), then get a point from a parametric equation between (0,0) and A-B, with x as the coefficient. then rotate (B-A) 90 degrees, and use y as the parametric coefficient of a line from (0,0) to the rotated point. Then, add these two points together, and transform the whole point back to A (by adding A to the point).

In matrix form, before simplification:

aaand the beautiful one after simplification:

(long-winded explanation to a simple equation :'D)

soooo now we can write a function:


2d_point transform(2d_point A, 2d_point B, 2d_point motifcoord)
{
 2d_point returnvalue;
 returnvalue.x=motifcoord.x*(B.x-A.x)-motifcoord.y*(B.y-A.y)+A.x;
 returnvalue.y=motifcoord.y*(B.x-A.x)+motifcoord.x*(B.y-A.y)+A.y;
 return returnvalue;
}

remember, the motif coordinate would be a point on a series of line segments starting at (0,0) and ending at (1,0). The motif for a koch curve would be:
{{0,0},{1/3,0},{.5,sin(60)*1/3},{2/3,0},{1,0}}
aaand so “motifcoord” would ALWAYS be one of the values in that array.

The algorithm running would look like this:

call the array containing the base’s coordinates “base”
call the array containing the motif’s coordinates “motif”
so base[0] is the first vertex of base, and motif[0] is the first vertex of the motif.

(1) create a list containing the final results (i’ll call it “myList”)
(2) copy the contents of the base into myList.
-(3) for each pair of points n and n+1 (where n+1 is at maximum the last point on the list):
–(a) add a point in myList at point n (same as saying “transform(point_n,point_nplusone,motif[0]”)
–(b) add a point in myList: transform(point_n,point_nplusone,motif[1]);
–(c) add transform(point_n,point_nplusone,motif[2]);
–…
–(i) add transform(point_n,point_nplusone,motif);
-(4) increase n, [i]taking special care: if the most recent point in myList is “point_nplusone”, and you increase n, it will now be “point_n”, that means that if you add “point_n” and “point_nplusone” all the time, then you will get duplicate vertices, and you’ll waste computing time. Only add point_nplusone if it is the last vertex.

-(5) repeat step (3) and upwards until you reach the desired iteration of the koch curve.

The algorithm is much, much simpler than the one you have. The one you have looks like it uses slope/intercept lines? Not nice things, those.

I have more about it including images I’ve rendered here:
http://www.neurofuzzydev.com/papers/basemotiffractals.php
this includes a snippet (in BASIC) of an implementation of the algorithm.

[I edited this post a bunch becaues some images didn’t want to show up]