Why this Graham Scan algorithm of convex hull works wrong when number of points more than 10k

[obry]

Here is the code. It works fine with 10k and less.
NUM variable it is number of points

import java.util.Arrays;
import java.util.Comparator;

ArrayList<PVector> init, sorted;
int NUM = 10000;
int border = 100;
void setup() {
	size(600, 600);
	background(205);
	
	//noLoop();
	// exit();
}
void newConvex(){
	init = new ArrayList<PVector>();
	// strokeWeight(2);
	int index = -1;
	float maxY = -1;
	float prevMax = maxY;
	for (int i = 0; i < NUM; i++) {
		PVector p = new PVector(border+random(width-2*border),border+random(height-2*border));
		init.add(p);
		// point(p.x, p.y);
		maxY=max(maxY, init.get(i).y);
		if (maxY!=prevMax) {
			index = i;
			prevMax = maxY;
		}
	}
	sorted = new ArrayList<PVector>();
	sorted.add(init.get(index));
	float[][] angles = new float[init.size()-1][2];
  PVector p0 = sorted.get(0);
	int j=0;
	for (int i = 0; i < init.size(); ++i) {
		PVector p = init.get(i);
		if (i!=index) {
			PVector a = PVector.sub(p,p0);
			angles[j][0] = a.heading();
			angles[j][1] = i;
			j++;
		}
	}
	Arrays.sort(angles, new Comparator<float[]>(){
		int compare(float[] a, float[] b) {
			return Float.compare(a[0],b[0]);
		}
	});
	sorted.add(init.get(int(angles[0][1])));
	for (int i = 0; i < init.size()-1; ++i) {
  if(i!=int(angles[0][1])) {
		PVector pp = sorted.get(sorted.size()-2);
		PVector p = sorted.get(sorted.size()-1);
		PVector p1 = init.get(int(angles[i][1]));
		PVector a = PVector.sub(p1,p);
		PVector b = PVector.sub(p,pp);
		PVector c = b.cross(a);
			while (c.z<0) {
        sorted.remove(sorted.size()-1);
				pp = sorted.get(sorted.size()-2);
				p = sorted.get(sorted.size()-1);
				a = PVector.sub(p1,p);
				b = PVector.sub(p,pp);
				c = b.cross(a);
			}
      sorted.add(init.get(int(angles[i][1])));
		}
    
		// text(i,init.get(int(angles[i][1])).x, init.get(int(angles[i][1])).y);
	}

	strokeWeight(1);
  stroke(0);
	noFill();
	beginShape();
	for (PVector p : sorted) {
		vertex(p.x, p.y);
	}
	endShape(CLOSE);
}

void drawPoints(){
  strokeWeight(1);
  stroke(0);
  noFill();
  for (PVector p : init) {
    point(p.x, p.y);
  }
}

void draw() {
  background(205);
	newConvex();
  drawPoints();
	
}

[josephh]

Hi @obry,

Welcome to the forum!

When I run the code it does that:

It looks like it’s working (convex hull of a dense point cloud approximate a square), could you describe precisely the behavior you are seeing?

[obry]

Yep. Everything works fine almost always. But as the number of points grows more than 10000 the more bugs like this come up

In this example I don’t draw points, only convex hull. The thing is that it is not convex as you see.
I wander if it is bug in algorithm. Maybe the reason is that some points are so close that this happens, which I doubt. Probably, the real reason is wrong algorithm

[josephh]

I tried with different values ranging from 10_000 to 100_000 and it seems to appear randomly. Maybe it’s a rounding error caused by the proximity of the points as you said so when you compare the angles, it breaks…

[quark]

Like @josephh I tried your code and as you increase the number of points the more frequent the artifacts. I even tried with 100,000 points and although the artifact appeared more often it didn’t happen all the time.

This tells me that the problem almost certainly not the algorithm, it is something to do with the data set.

I suspect that some of the intermediate calculations are producing NaN(s) I had this problem using floats and had to uses doubles instead.

[quark]

In your method newConvex you have this conditional statement

@obry

if (maxY!=prevMax) {

You should test 2 floating point number for equality (==) or inequality (!=) because it can lead to logical errors. What would you expect from this sketch?

float a, b, c;
a = 3.1415927;
b = 3.1415928;
c = 3.1415926;
println(a == b);
println(a == c);

You get

true
false

which is obviously wrong.

You can improve the first part of the newConvex method to avoid this issue and also improve efficiency by avoiding reading from the list like this -

void newConvex(){
  init = new ArrayList<PVector>();
  // strokeWeight(2);
  int index = -1;
  float maxY = -1;
  for (int i = 0; i < NUM; i++) {
    PVector p = new PVector(border+random(width-2*border),border+random(height-2*border));
    init.add(p);
    if (p.y > maxY) {
      index = i;
      maxY = p.y;
    }
  }

BTW this does not get rid of the visual artifact unfortunately

[obry]

Thank you, guys!
Seems like this is the case.
I appreciate your help.