1

With reference to this answer to a post concerning finding which points in a set are inside or outside a convex hull: https://stackoverflow.com/a/4903615/6376297

The suggested technique relies on testing, for a point P, whether all the vectors P-H_i, where H_i is the i-th point on the hull, are on the 'same side' of a (hyper)plane.
The post however does not explain how to test that in practice. One of the comments says:

To check for "on less than 1/2 of a hypersphere" you check the dot product between the two directions and see if any 2 are negative (i.e. in opposite direction)?

And the answer by another user to that was:

It is not sufficient to do that, Evans. But you could do this: sum all vectors from the point to the other points. Then, use that summed vector as axis to project again all the vectors from the point in question. If all projections lay on the same side of the axis, the point is outside the hull.

It seemed reasonable, so I tried it, and unfortunately it failed in some cases, even in 2D.
See the example below (R):

# make a set of 15 points in 2D
pss <- structure(c(-1.00537976836315, -0.87500222326089, 1.24730199591969, 0.10896916887902, -0.019985323116284, 1.07410848070838, -0.587685189152222, 0.511733218973978, 0.288174663946251, 0.647930319959558, 1.85518190486922, -0.0773021851847363, -0.597293019621372, -1.15643543271745, -1.41431661184, -0.665028124047381, 1.54922010009077, -0.324241735276436, -0.0788617244295097, 1.26736978063882, 0.914849025731909, 0.491408952001827, -1.03783441172396, -0.0468133979737429, -1.06852848099562, -1.22208926690273, -2.25715256174601, 1.61782957811296, -0.152370162480533, 1.01224242899964), .Dim = c(15L, 2L), "`scaled:center`" = c(0.1286127953793, 0.286128946770637))

# represent the points in a 2D plot
plot(pss)
text(pss,label=1:N,pos=4)

# collect the indices of all points
p <- 1:N
# define the indices of the points that are on the convex hull, and remove them from p
b <- c(15,5,6,3,12)
# remove the points in b from p
p <- p[!(p %in% b)]
NB = length(b)

# colour in black the points that are on the hull and connect them by solid lines
points(pss[b,],col=1,pch=16)
lines(rbind(pss[c(b,b[1]),,drop=F]))

# store the length (or norm) of each vector
pss_lengths <- apply(pss,1,function(v) sqrt(sum(v^2)))

# test for each point if it is inside or outside the hull, based on: https://stackoverflow.com/a/4903615/6376297

pb0 <- apply(pss[b,,drop=F],2,sum)
ip = 1

while (ip <= length(p)) {

  flag = 0
  ip1 = p[ip]
  p1 = pss[ip1,,drop=F]
  p1n = p1 / pss_lengths[ip1]
  v1i_sum = pb0 - NB*p1
  v1i_sum_n = v1i_sum / sqrt(sum(v1i_sum^2))
  lines(rbind(p1,p1+v1i_sum_n/5),lty=1,col=2)

  # test all the points on the hull (b)
  for (ib in 1:NB) {
    ibi = b[ib]
    pbi = pss[ibi,,drop=F]
    v1i = pbi - p1
    v1in = v1i / sqrt(sum(v1i^2))
    lines(rbind(p1,p1+v1in/5),lty=2,col=4)
    prpr = sum(v1in*v1i_sum_n)
    if (ib == 1) {sign_prpr = as.numeric(prpr >= 0)} else {
      new_sign_prpr = as.numeric(prpr >= 0)
      if (new_sign_prpr != sign_prpr) {
        cat("\nPoint ",ip1," is inside the convex hull (angle with point ",ibi,").")
        # colour in red the point that was found to be inside the hull
        points(p1,pch=16,col=2)
        flag = 1
        break
      }
    }
  }

  ip = ip + 1; next  

}

This should result in this plot:

enter image description here

As you can see, the method correctly detects that points 2, 11 and 13 are outside the defined convex hull, and that points 4, 7, 8, 9 and 10 are inside it (coloured in red). However, it wrongly concludes that points 1 and 14 are inside the hull (because the angle between at least one P-H_i and the sum of all P-H_i's is obtuse).
If I understand correctly, this is due to the fact the the sum of P-H_i's is often biased due to where the points of the hull are located, so it doesn't offer a good reference for testing whether all vectors 'point in the same half-plane'.

Still, from the picture it is clear that all vectors departing from any P outside the hull are indeed contained in a specific half-plane.

Question: can you suggest any way to make this work?

The only way I can think of is to take as the reference the vector bisecting the largest angle between any pair of P-H_i vectors; that vector is guaranteed to have no angles > 90 degrees when P is outside the hull. But I would not know how to find that vector without looping over all pairs of points on the hull.

Thanks!

user6376297
  • 575
  • 2
  • 15

1 Answers1

0

You are on the right track to think about using angles. But you can do it in a simpler way without looping over the vectors an extra time. Stick to a convention for measuring angles between two vectors - say, angles measured counterclockwise are positive and those measured clockwise are negative. If you measure the angle between two vectors based on dotproduct, it will always be between 0 and pi. You can then enforce the convention by attaching a +ve or -ve sign to this number based on the crossproduct of the two vectors, i.e. +ve if the crossproduct is positive and -ve if the cross product is negative.

(This is pseudocode since I am not familar with R)

#loop over the hull points
for each input_point ip{
    left = first_hull_point - ip
    right = first_hull_point - ip
    for each hull_point h{
        #skip the first hull point
        vi = h - ip # vector from ip to h
        angle = angle between vi and left
        if (angle < 0){left = vi}
        angle = angle between right and vi
        if (angle < 0){right = vi}
    }
    angle = absolute(angle between left and right) # absolute value
    if (angle < pi) { "The point is outside" }
    else if (angle == pi) { "The point is on the edge" }
    else { "Point is inside" }
}

Note: Because of our angle measuring convention that uses cross product, the order of vectors is important. Flipping the order changes the sign of the angle. Keep this in mind when measuring the angles (between vi and left) and (between right and vi).

Ranjeeth Mahankali
  • 391
  • 2
  • 10