Implementing Robert Hodgin's "Meander": trouble with curvature and 'modified bitangents'

p5.js Coding Questions
1

Hi!

I’ve been trying to create a very basic version of Robert Hodgin (aka Flight 404) 's Meander program but I am having some difficulty and I’m not sure what the problem is.

The sketch I am working on is here: https://www.openprocessing.org/sketch/904851.

I have been using his write up and this how I have implemented the steps:

  1. Implement the “bitangent”. (I would normally refer to this as a normal rather than a bitangent–is this right? I’ll stick to his wording though so as not to complicate this.)

    His implementation: http://roberthodgin.com/v2/wp-content/uploads/2020/05/process_2_bitangent.png
    My implemenation: https://imgur.com/a/MDIDcAc

    This seems to be working fine. This is the code I have been using for it:

//m=Δy/Δx
let tangentGradient = (currentPoint.y - previousPoint.y) / (currentPoint.x - previousPoint.x);

//convert from gradient to angle
let tangentAngle = atan(tangentGradient) + ((currentPoint.x - previousPoint.x) < 0 ? -PI : 0);
// if (tangentAngle > TAU / 4 && tangentAngle < TAU * 3 / 4) tangentAngle -= PI;

//rotate 90 degrees
let bitangentAngle = tangentAngle - HALF_PI;

let bitangent = createVector(1, 0).rotate(bitangentAngle);

  1. Modify the bitangent by scaling according to the curvature at that point. This seemed to work for me too. I used this code to find the curvature and to find out whether it is negative or positive I checked to see if the change in angle was positive or negative.

    This is the part of the code where I calulate the curvature:

//-----calculate curvature-----
   
   let curvatures = [];

   //leave out the beginning and end points as they only have one neighbour
   for (let i = 1; i < riverPoints.length - 2; i++) {

   	//calulate the curvature at point i using the method given here: 
   	//https://uk.mathworks.com/matlabcentral/answers/284245-matlab-code-for-computing-curvature-equation#answer_222173

   	let a = p5.Vector.dist(riverPoints[i - 1], riverPoints[i]);
   	let b = p5.Vector.dist(riverPoints[i], riverPoints[i + 1]);
   	let c = p5.Vector.dist(riverPoints[i - 1], riverPoints[i + 1]);

   	//area of the triangle that connects the three points
   	s = (a+b+c)/2;
  	A = sqrt(s*(s-a)*(s-b)*(s-c)); //Area of triangle
   	
   	//calulate whether the curvature is positive or negative
   	let gradientA = (riverPoints[i].y -riverPoints[i-1].y)/(riverPoints[i].x -riverPoints[i-1].x);
   	let angleA = atan(gradientA) + ((riverPoints[i].x -riverPoints[i-1].x) < 0 ? -PI : 0);
   	let gradientB = (riverPoints[i+1].y -riverPoints[i].y)/(riverPoints[i+1].x -riverPoints[i].x);
   	let angleB = atan(gradientB) + ((riverPoints[i+1].x -riverPoints[i].x) < 0 ? -PI : 0);
   	let angleChange = angleB - angleA;
   	if(angleChange > TAU/2) angleChange -= TAU;
   	if(angleChange < -TAU/2) angleChange += TAU;
   	if(angleChange < 0) A*=-1;

   	//curvature of the circumscribing circle (K)
   	curvatures[i] = 4 * A / (a * b * c);
   }

I then add on the curvature to the bitangent:
let modifiedBitangent = p5.Vector.mult(bitangent, curvature);.

This seems to work.
This is what I get: https://imgur.com/a/v0pk02Y and this is what Robert shows in his write-up: http://roberthodgin.com/v2/wp-content/uploads/2020/05/process_2_bitangent_mod.png.

  1. Add the tangent to the modified bitangent to get the velocity (or the blended directional vector as Hodgin calls it). I am not adding the tangent to the velocity here (the velocity is just the modified bitangent) as I am pretty sure the problem lies in the modified bitangent not the tangent part.

  2. Move the points by adding the velocity to each of them.

So the problem is this:

Once I start trying to move the river (by adding the velocity to each point), after a few iterations I get these crazy zig-zag patterns: https://imgur.com/a/4zyOUFs and after a while it becomes this: https://imgur.com/a/PkL5Gqx

Why is it doing this? Before I start moving the points everything seems to be calculated properly but once it starts I get these wacky patterns.

Full code from the Open Processing sketch linked above 
const speed = 1;
const numberOfPoints = 400;

let riverPoints = [];

function setup() {
	createCanvas(windowWidth, windowHeight);

	//intital set of points
	//used perlin noise so that the line we start off with is a bit wiggly
	for (let x = 0; x <= width; x += width/(numberOfPoints-1)) {
		let y = height / 2 + map(noise(x * 0.002), 0, 1, -400, 400);
		riverPoints.push(createVector(x + map(noise(x * 0.002, 345), 0, 1, -100, 100), y));
	}

	strokeWeight(2);
	stroke("#f5ea80");
	
}

function draw() {

	background("#3A2D3D");

	//-----calculate curvature-----
	
	let curvatures = [];

	//leave out the beginning and end points as they only have one neighbour
	for (let i = 1; i < riverPoints.length - 2; i++) {

		//calulate the curvature at point i using the method given here: 
		//https://uk.mathworks.com/matlabcentral/answers/284245-matlab-code-for-computing-curvature-equation#answer_222173

		let a = p5.Vector.dist(riverPoints[i - 1], riverPoints[i]);
		let b = p5.Vector.dist(riverPoints[i], riverPoints[i + 1]);
		let c = p5.Vector.dist(riverPoints[i - 1], riverPoints[i + 1]);

		//area of the triangle that connects the three points
		s = (a+b+c)/2;
   	A = sqrt(s*(s-a)*(s-b)*(s-c)); //Area of triangle
		
		//calulate whether the curvature is positive or negative
		let gradientA = (riverPoints[i].y -riverPoints[i-1].y)/(riverPoints[i].x -riverPoints[i-1].x);
		let angleA = atan(gradientA) + ((riverPoints[i].x -riverPoints[i-1].x) < 0 ? -PI : 0);
		let gradientB = (riverPoints[i+1].y -riverPoints[i].y)/(riverPoints[i+1].x -riverPoints[i].x);
		let angleB = atan(gradientB) + ((riverPoints[i+1].x -riverPoints[i].x) < 0 ? -PI : 0);
		let angleChange = angleB - angleA;
		if(angleChange > TAU/2) angleChange -= TAU;
		if(angleChange < -TAU/2) angleChange += TAU;
		if(angleChange < 0) A*=-1;

		//curvature of the circumscribing circle (K)
		curvatures[i] = 4 * A / (a * b * c);
	}
	
	//-----calculate tangent, bitangent and modified bitangent and draw lines-----
	
	let previousPoint = createVector(riverPoints[0].x - width/numberOfPoints, riverPoints[0].y + 0.001);

	for (let i = 0; i < riverPoints.length - 1; i++) {

		let currentPoint = riverPoints[i];

		//draw line
		line(currentPoint.x, currentPoint.y, previousPoint.x, previousPoint.y);

		//m=Δy/Δx
		let tangentGradient = (currentPoint.y - previousPoint.y) / (currentPoint.x - previousPoint.x);

		//convert from gradient to angle
		let tangentAngle = atan(tangentGradient) + ((currentPoint.x - previousPoint.x) < 0 ? -PI : 0);
		// if (tangentAngle > TAU / 4 && tangentAngle < TAU * 3 / 4) tangentAngle -= PI;

		//rotate 90 degrees
		let bitangentAngle = tangentAngle - HALF_PI;

		let curvature = curvatures[i] ?? 0;
		//curvature = curvature > 0 ? 1 : -1;

		let bitangent = createVector(1, 0).rotate(bitangentAngle);

		let modifiedBitangent = p5.Vector.mult(bitangent, curvature);
		//modifiedBitangent.mult(10);

		let tangent = createVector(1, 0).rotate(tangentAngle);
		//tangent.setMag(1);

		let velocity = p5.Vector.add(modifiedBitangent, tangent);
		//velocity.setMag(1);

		//draw the modified bitangent line
		stroke("#D3414D"); //#1DC1A1
		line(currentPoint.x, currentPoint.y, currentPoint.x + modifiedBitangent.x * 700, currentPoint.y + modifiedBitangent.y * 700);
		stroke("#f5ea80");

		previousPoint = currentPoint.copy();

		velocity = modifiedBitangent;
		
		velocity.mult(speed);
		currentPoint.add(velocity);

	}
}

function keyPressed() {
	if (key === "s") save();
}

Sorry for the external links to images, I couldn’t get inline images to work.

1 Like