Four points are chosen at random inside a square. Each point is chosen by choosing a random x-coordinate and a random y-coordinate.

A convex quadrilateral is drawn with the the four random points as the vertices.

Determine the probability that the center of the square is inside this quadrilateral.

def angl(dy,dx):
    if dx == 0:
        if dy ==0:
            return 0
        elif dy > 0:
            return 90
        elif dy < 0:
            return 270
        
    a = math.atan(dy/dx) * 180 / math.pi
    if dy >= 0 and dx >= 0:  # quadrant I
        return  a
    elif dy >= 0 and dx < 0:  # quadrant II
        return  a + 180
    elif dy < 0 and dx < 0:  # quadrant III
        return a + 180
    elif dy < 0 and dx >= 0:  # quadrant IV
        return a + 360
    return

def isconvex(coords):
    """ bool  can these 4 points form a convex quadrilateral?  
    coords is [x1,y1,x2,y2,x3,y3,x4,y4]    
    True if no point is inside the triangle formed by the other 3 points """
    inside = -1
    for i in range(4):
        startX = coords[2*i]
        startY = coords[2*i+1]
        jjj = [k for k in range(4) if k != i]
        angles = []
        for j in jjj:
            vector = angl(coords[2*j+1] - startY ,  coords[2*j] - startX)
            angles.append(vector)
        a = (angles[0] + 180)%360
        b = (angles[1] + 180)%360
        arc = abs(a-b)
        if arc <= 180:
            lo = min(a,b)
            hi = max(a,b)
        elif arc > 180:
            lo = max(a,b)
            hi = 360 + min(a,b)
        if angles[2] > lo and angles[2] < hi:
            inside += 1
    return inside == -1

def isinside(coords):
    """ bool  is the center inside any of the 4 triangles you can form from these 4 points?  
    coords is [x1,y1,x2,y2 etc ]    """
    inside = -1
    for i in range(4):
        startX = 0
        startY = 0
        jjj = [k for k in range(4) if k != i]
        angles = []
        for j in jjj:
            vector = angl(coords[2*j+1] - startY ,  coords[2*j] - startX)
            angles.append(vector)
        a = (angles[0] + 180)%360
        b = (angles[1] + 180)%360
        arc = abs(a-b)
        if arc <= 180:
            lo = min(a,b)
            hi = max(a,b)
        elif arc > 180:
            lo = max(a,b)
            hi = 360 + min(a,b)
        if angles[2] > lo and angles[2] < hi:
            inside += 1
    return inside > -1

# RANDOM METHOD
size = 1000
reps = 10000000
insideConvexCount = 0
outsideConvexCount = 0
insideHullCount = 0
outsideHullCount = 0

from random import randint
for i in range(reps):
    coordinates = []
    angles = []
    for k in range(4):
        coordinates.append(randint(-size,size))
        coordinates.append(randint(-size,size))
    
    if isconvex(coordinates):
        if isinside(coordinates):
            insideConvexCount += 1
        else:
            outsideConvexCount += 1
    else:
        if isinside(coordinates):
            insideHullCount += 1
        else:
            outsideHullCount += 1

print(insideConvexCount, outsideConvexCount, insideHullCount, outsideHullCount )
print(insideConvexCount /  (insideConvexCount + outsideConvexCount),
      insideHullCount / (insideHullCount+outsideHullCount) )
print()

# GRID METHOD
size = 4
for a in range(-size,size+1):
    for b in range(-size,size+1):
        if a==0 and b==0:
            continue
        for c in range(-size,size+1):
            for d in range(-size,size+1):
                if c==0 and d==0:
                    continue
                for e in range(-size,size+1):
                    for f in range(-size,size+1):
                        if e==0 and f==0:
                            continue
                        for g in range(-size,size+1):
                            for h in range(-size,size+1):
                                if g==0 and h==0:
                                    continue
                                coordinates = [a,b,c,d,e,f,g,h]
                                if isconvex(coordinates):
                                    if isinside(coordinates):
                                        insideConvexCount += 1
                                    else:
                                        outsideConvexCount += 1
                                else:
                                    if isinside(coordinates):
                                        insideHullCount += 1
                                    else:
                                        outsideHullCount += 1
print(insideConvexCount, outsideConvexCount, insideHullCount, outsideHullCount )
print(insideConvexCount /  (insideConvexCount + outsideConvexCount),
      insideHullCount / (insideHullCount+outsideHullCount) )

Posted by Larry on 2023-01-12 10:45:24