CG · Maths · Python

Barycentric Coordinates

Barycentric Coordinates is a type of coordinate system to define the coordinates of an entity in “n-dimensional space” (for CG it is 3 dimensions). One of the interesting properties of this coordinate system is that it is affine and depends on the object space.

Barycentric Coordinates were formulated  by Mobius (who is also known for the famous Mobius strip) in 1827.

This coordinate system bears significant importance especially in Computer Graphics but it can be a daunting task to understand fully the concept of Barycentric Coordinates for not so mathematically inclined.

So in this post, I will try to explain in the simplest terms, what Barycentric Coordinates are and how they are used in Computer Graphics.

As is the case with many other DCC Applications, Barycentric Coordinates are exploited to a great extent in Softimage. Although I have no official proof, but I can most safely say that the Location Data we so much love and use in ICE are nothing but an implementation of Barycentric Coordinates.

To begin with, Consider a triangle.

1

If I was asked to find a point on this piece so that it could be balanced perfectly on a point of a pencil, what would I do? A little bit of intuition would suggest taking the central point of the triangle.

2

That is all correct, but let’s now complicate things a little bit. If the shape of the piece is an irregular triangle with all sides having different lengths then the central point is not easier to guess.

3

I would solve this by applying a little geometry and divide each side in half and connect the opposite vertex of the triangle to these points. The lines thus created will yield the central point and I will be able to balance the triangle on that point with the pencil.

4

All good, now let us raise the bar some more and tie some weights to each vertex of the triangle. There are different weights on each vertex. Now the task of finding a point at which the triangle can be balanced is no more pure geometry. I have to consider that each of these weights is pulling the triangle with different forces.

5

This is where the Barycentric Coordinates can be applied to find the answer. In fact, the word barycentre originates from the Greek root “barus” which generally refers to weighty or heavy.

So to solve our problem we should divide each side of the triangle in the same ratio as the weights in the vertex forming that side, as shown below.

After getting this three points, it is just of matter of connecting this points to the opposite vertex and bingo! We have our balancing point:

6

 

What we just did was to find the Barycentric coordinates of the point with given weights at the vertices.

We can also write this as:

mathtex1

 

 

Where P is our barycentric point and A, B and C are the weights at the vertex. 10 comes from adding the three weights 2 + 3 + 5. Also, note that:

mathtex2

 

 

So now that we have the notion of barycenter, let us see how this relates to CG coordinates with vector geometry. We know that any triangle in the 3D Coordinate system can be defined as (x1, y1, z1), (x2, y2, z2), (x3, y3, z3) representing the x, y,z-coordinates of each vertex. Also in vector notation, this can be written as, A, B, C where A, B and C are the vectors defining the point positions (vertices) of the triangle.

Let P be any arbitrary point on the triangle, then P can always be defined with the help of vectors A, B and C as :

mathtex3

Where,

mathtex4

Here ,

mathtex5

are called the Barycentric coordinates (remember these are just three real numbers for a triangle).

So now we totally understand the notion as well as the mathematical representation of Bary Centric Coordinates. How do we use it to make life easier in the CG World?

There are numerous ways to use them, but here I will discuss the locations in Softimage which are available both in ICE and the API methods.

The whole idea of locations magically following the mesh (or sticking to location) happens in the following way:

  1. First a point on the triangle is found which will be used as initial location. This can be done in many ways, raycasting is one of them. We throw a ray from a certain point, in a well defined direction towards the triangle. The point where the ray hits(intersects) with the triangle is determined.
  2. Once the coordinates(vector) of the point on the triangle is known, Barycentric coordinates this point are calculated (this calculation is not very hard, but still requires some mathematical knowledge, so I am skipping them for the purpose of being simple).
  3. Once we know the Barycentric coordinates for this point, we can simply get this point at any point of time just using the three position vectors of the triangle.
  4. In this way, even when the triangle is deformed in any way, we always know our point, in other words, the location of the point sticks to the triangle.
  5. For example, after doing a ray cast we find that barycentric coordinates for the location are (0.35, 0.25, 0.40). These will be fixed throughout the simulation/animation/mesh move. After the triangle has been transformed, to get our location, we have to simply get the vectors for the point position at this point and multiply them by the scalars 0.35, 0.25 and 0.40, add the result. This will give a new vector, which is our location.

This ICE Tree exhibits the barycentric coordinates in action. Even if you move the point of the triangle, the point will always lie in the triangle.

 

Download Scene File here. In the softimage scene, move the vertices of the triangle and you will see that the location represented by the red sphere will remain on the triangle:

8

 

I hope that now you have a better idea of what barycentric coordinates are.

Advertisement
CG · Maths · Python

Convex Hull – Maths can be fun

I always strive to break the mantle of sophistication that is so attached with mathematics. In fact in all of the sciences, if the fundamentals are understood clearly, in simple terms, it is easy for anybody to understand the advanced concepts. But we have to have a solid fundamental base to start with. I always maintain that if we can learn the 26 alphabets of English, we can learn almost anything !!

Recently, I am working on creating a new algorithm for finding the convex hull from a given set of 2d points in a plane. In this post, we will discover the  concept of the convex hull. It will be a fun ride. In fact, the method that I devised for finding the convex hull is so easy that even a six-year-old can do it given a piece of paper and pencil.

So first thing first, what is a Convex Hull ?

Convex Hull

If we keep ourselves restricted to the computational geometry and 2 dimensions:

Convex Hull can be defined as the boundary polygon for a given set of points.

It is clearer from the illustration below.

convexHull

Recently, I am working on creating a new algorithm for finding the convex hull from a given set of 2d points in a plane. In this post, we will discover the  concept of the convex hull. It will be a fun ride, do believe me. In fact, the method that I devised for finding the convex hull is so easy that even a six-year-old can do it given a piece of paper and pencil.

 

So first thing first, what is a Convex Hull ?

 

Convex Hull

If we keep ourselves restricted to the computational geometry and 2 dimensions:

 

Convex Hull can be defined as the boundary polygon for a given set of points.

 

It is clearer from the illustration below.

convexHull

You can also imagine the convex hull by the following analogy. Consider this points as nails on a piece of wood. If you were to wrap a rubber band around these nails, it will snap to the boundary nails. The shape of the rubber band represents the convex hull for the given nail positions.

The convex hull has numerous application in CG. Namely, better bounding boxes for meshes, BSP calculations for render times, finding projection areas to drive asset resolutions and numerous other. Creativity has no bounds, you can use convex hulls to do many things, use your imagination.

Let’s  start,

Suppose we have a set of points in a given 2D plane. We want to find out a way to find all the points that define the boundary of all the given points.

As I said, it is going to be fun. Let us start with a few random points in a 2d plane. We assign a name to each of this points.

convexHullNamed1

Looking at the figure above, we prepare two lists of these points. First, going to from left to right, we take the left most point, and then the next left most point and so on, then we make another list going from top to bottom for the same points.

These are the two list that we will have:

 

Left to Right      : [f, j, d, a, b, h, l, m, c, g, i, e, k]

Top to Bottom  : [h, m, e, j, a, c, b, k, f, i, g, d, l]

The cool thing now is we have all the information we need to process the convex hull. Just grab a paper and pencil and play with these lists according to the rules of the game I define below. In fact, what we will be doing from now is the gist of my algorithm. Remember, we are working with alphabets, so no calculations are needed, not even an addition. It is so simple. Even the alphabets do not represent any order. They are just a bunch of symbols to distinguish between the points,  you can use smileys if you want instead !! and it will still work. We also don’t need to look at the figure even. Just the lists.

 

OK, so let’s start the game, here’s how we play:

 

We start from the Left-Right list. We take the first element of this list. In this ‘f’. This definitely is a border point. So our first point of the convex hull is ‘f’.

 

Convex Hull = [f, ]

 

Next, we look at the point next to ‘f’ in the Left-Right list, it is ‘j’. Now we should see whether ‘j’ also comes after ‘f’ in the Top-Bottom list. No. So we move to the next point in our Left-Right list, it is ‘d’. Does ‘d’ comes after ‘f’ in the Top-Bottom list. Yes, it does. So check-in this point next in our convexHull,

 

Convex Hull = [f, d, ]

As soon as we check-in a point in the Convex Hull list. We should remove the previous ‘f’ from the Left-Right and Top-Bottom lists. So these now become.

 

Left to Right      : [ j, d, a, b, h, l, m, c, g, i, e, k]

Top to Bottom  : [h, m, e, j, a, c, b, k, i, g, d, l]

 

So our focus now moves to ‘d’. What is the next point after ‘d’ in the Left-Right list –  ‘a’. Does ‘a’ comes after ‘d in the Top-Bottom list. No. So we move on, till we reach ‘l’. This comes after ‘d’ in both Left-Right and Top-Bottom list and so this is a convex hull point.

Convex Hull = [f, d, l, ]

 

again we remove previous point ‘d’ from both the lists,

 

Left to Right      : [ j, a, b, h, l, m, c, g, i, e, k]

Top to Bottom  : [h, m, e, j, a, c, b, k, i, g, l]

 

Now as we have reached the last point of the Top-Bottom list. We can no longer check whether there is any point after it. Whenever we reach the last point in any of the list. We should switch the list and change direction. Since we have reached ‘l’, which is the last point of the Top-Bottom list, we should now start looking for points in the Top-Bottom list but in the reverse direction. So what point comes before ‘l’ in the Top-Bottom list. ‘g’, does ‘g’ also comes after ‘l’ in the Left-Right list, yes. So we check in ‘g’ to the Convex Hull.

 

Convex Hull = [f, d, l, g]

 

removing ‘l’ we get,

Left to Right      : [ j, a, b, h, m, c, g, i, e, k]

Top to Bottom  : [h, m, e, j, a, c, b, k, i, g]

 

You get the idea, again, continuing in the same fashion we will have the convex hull as :

Convex Hull = [f, d, l, g, i, k]

 

and the lists as :

Left to Right      : [ j, a, b, h, m, c, e, k]

Top to Bottom  : [h, m, e, j, a, c, b, k]

when we reach ‘k’, this is the last element in the Left-Right list. So according to our rules, we switch the lists and start to look in the reverse in the Left-Right list.  At this moment, we are looking in both the lists in the reverse directions remember. So the point before ‘k’ in Left-Right list is ‘e’, does it also come before ‘k’ in the Top-Bottom list, the answer is yes and so we check in ‘e’ in the convex hull and remove the previous point ‘k’ from the list,

Convex Hull = [f, d, l, g, i, k, e]

Left to Right      : [ j, a, b, h, m, c, e]

Top to Bottom  : [h, m, e, j, a, c, b]

 

moving in the same fashion in both the lists we have,

 

Convex Hull = [f, d, l, g, i, k, e, m, h]

Left to Right      : [ j, a, b, h, c]

Top to Bottom  : [h,  j, a, c, b]

 

‘h’ being the first point of the Top-Bottom list, we again switch the lists and look in the Top-Bottom list in the forward direction, while reverse in the Left-Right list. And we keep on going to find the same patterns of occurrence, we are almost there,

 

in the Top-Bottom list, the element next to ‘h’ is ‘j’. Does ‘j’ comes before ‘h’ in the Left-Right list. Yes, so we check in ‘j’ and remove ‘h’ from the lists. At this point, we stop. Why ? because we have reached where we started, the first element of the Left-Right list. So our convex hull at this point look like.

 

Convex Hull = [f, d, l, g, i, k, e, m, h, j]

 

Let’s draw our hull and check what we have,

convexHullRough

Not bad, without any calculation. But we have a few concave points still inside. These are ‘g’ , ‘i’ and ‘m’. These points might or might not  appear in the hull depending upon the uniqueness of our data set. But we should always as a final step remove these out of our way. To do this we take three points of the hull at a time and check for convexity. Please refer to my post for more details finding convex points.

So we take three point sets from the hull in the order they appear, like ‘fdl’, ‘dlg’, ‘lgi’ and so on and so forth and remove all the concave points, we call this a single iteration. We do a few iterations of the same logic till we have no more convex points left.

This is how these points are removed in the iteration :

First Iteration : ‘i’ is removed from ‘gik’, ‘m’ is removed from ‘emh’.

convexHullFirstIteration

Second Iteration : ‘g’ is removed from ‘lgk’.

convexHullSecondIteration2

And lo! we finally have our convex hull as:

 

Final Convex Hull = [  f, d, l, k, e,  h, j  ]

convexHullFinal

Let’s use our algorithm inside Softimage to find the convex hull for meshes. To use the algorithm, we use a python implementation of the algorithm and feed to it the position array of a mesh. We create a curve based on the convex hull. Just select any mesh and run the below code. It will create a curve representing the convex hull.

 

Note that as we are on a 2D plane, we get the hull as we see the points in front view. Even if you rotate, translate or scale the mesh, the hull will still be created at the current position array. To create the hull for a different view of the mesh, select the points of the mesh and then rotate, scale or transform them. After applying your transformations to the mesh points, freeze the mesh.

 

Here are a few screenshots.

meshHulls

Please note that for the implementation of the algorithm in python, I am not using any optimizing techniques. This code creates a convex hull for 15000 mesh points in about 6 seconds. To deploy this as a production ready solution, I would use numpy, cython or build my own extensions written in c/c++ for performance enhancement and better memory management to handle a few million points. Since that is beyond the scope of this post, so let’s not discuss it.

This  is the  implementation:

from win32com.client import constants as c
import random

def printOnOff(inStr):
    printIt = False
    if printIt:
        print inStr
    else:
        pass

    return

def makeNulls(xCoords, yCoords, sortedPoints,curve):

    nbPoints = len(xCoords)
    mdl = Application.ActiveSceneRoot.AddModel(None,'randomPoints')
    mdl.Properties('Visibility').Parameters('viewvis').Value = False

    points = mdl.AddNull('points')
    points.Properties('Visibility').Parameters('viewvis').Value = False

    if curve:
        points.AddChild(curve)    

    for i in range(nbPoints):
        if i < 10:
            zeros = '000'
        elif i < 100:
            zeros = '00'
        elif i < 1000: zeros = '0' else: zeros = '' null = points.AddNull('point_%s%s'%(zeros, i)) null.Kinematics.Global.Parameters('PosX').Value = xCoords[i] null.Kinematics.Global.Parameters('PosY').Value = yCoords[i] null.Kinematics.Global.Parameters('PosZ').Value = 0 null.Parameters('primary_icon').Value = 0 null.Parameters('shadow_icon').Value = 5 null.Parameters('shadow_colour_custom').Value = True if i in sortedPoints: r = 0.0 else: r = 1.0 null.Parameters('R').Value = r null.Parameters('G').Value = 1.0 null.Parameters('B').Value = 0.0 return def createRandomPoints(nbPoints=20): xCoords = [] yCoords = [] for i in range(nbPoints): xRand = nbPoints * random.random() #* random.randint(1,3) yRand = nbPoints * random.random() #* random.randint(1,3) x = random.uniform(-1 * xRand, xRand) y = random.uniform(-1 * yRand, yRand) xCoords.append(x) yCoords.append(y) return xCoords, yCoords def makeOrderedTuples(xCoords, yCoords): data = None nbPoints = len(xCoords) xCoordsTuples = [] yCoordsTuples = [] for i in range(nbPoints): xRand = nbPoints * random.random() #* random.randint(1,3) yRand = nbPoints * random.random() #* random.randint(1,3) x = xCoords[i] y = yCoords[i] xCoordsTuples.append((x, i)) yCoordsTuples.append((y, i)) sortedXIndex = [] sortedXIndex = [ xTuple[1] for xTuple in sorted(xCoordsTuples, key=lambda xCoord: xCoord[0])] sortedYIndex = [] sortedYIndex = list(reversed([ yTuple[1] for yTuple in sorted(yCoordsTuples, key=lambda yCoord: yCoord[0])])) data = xCoordsTuples, yCoordsTuples, sortedXIndex, sortedYIndex return data def makeCoordString(coords): xCoords, yCoords,zCoords = coords pStr = '' nbPoints = len(xCoords) for i in range(nbPoints): pStr += '(' pStr += str(xCoords[i]) pStr += ',' pStr += str(yCoords[i]) pStr += ',' pStr += str(zCoords[i]) pStr += ')' if i != nbPoints - 1: pStr += ',' return pStr def makeCurve(pStr): curve = Application.CreateCurve( 1, 0, pStr)('Value') curve.Name = 'convexHull' Application.ApplyTopoOp('CrvOpenClose', curve, 3, c.siPersistentOperation) Application.FreezeObj(curve) #curve.Properties('Visibility').Parameters('selectability').Value = False return curve def isConvex(pointA, pointP, pointB): x1, y1 = pointA x2, y2 = pointB x3, y3 = pointP quadB = 0 if x2 > x1 and y2 > y1:
        quadB = 0
    elif x2 < x1 and y2 > y1:
        quadB = 1
    elif x2 < x1 and y2 < y1: quadB = 2 elif x2 > x1 and y2 < y1: quadB = 3 y = y1 + ((x3 - x1) * (y2 - y1) / (x2 - x1)) if quadB == 0 or quadB == 3: if y > y3:
            return 1
        else:
            return 0
    elif quadB == 1 or quadB == 2:
        if y > y3:
            return 0
        else:
            return 1

def getConvexHull(xCoords, yCoords):
    xCoords, yCoords, sortedX, sortedY = makeOrderedTuples(xCoords, yCoords)
    startPoint = sortedX[0]
    firstPoint = startPoint
    sortedPoints = [firstPoint]
    inflection = 0
    normalInflection = False
    inflectionPoints = []
    nextPoint = 0
    nbPoints = len(sortedX)
    for i in range(nbPoints):
        yIndexOfFirstPoint = sortedY.index(firstPoint)
        xIndexOfFirstPoint = sortedX.index(firstPoint)
        for xIndex in range(xIndexOfFirstPoint + 1, nbPoints):
            pointToCheck = sortedX[xIndex]
            yIndexOfPointToCheck = sortedY.index(pointToCheck)
            if yIndexOfPointToCheck > yIndexOfFirstPoint:
                nextPoint = pointToCheck
                sortedPoints.append(nextPoint)
                if yIndexOfPointToCheck + 1 == nbPoints:
                    printOnOff('First Inflection at %s !!'%nextPoint)
                    inflectionPoints.append(nextPoint)
                    inflection += 1
                    normalInflection = True
                break
        if inflection == 1:
            break

        firstPoint = nextPoint

    if not normalInflection:
        printOnOff('First Inflection at %s - Abnormal!!'%firstPoint)
        inflection += 1

    normalInflection = False
    firstPoint = nextPoint
    for i in range(nbPoints):
        xIndexOfFirstPoint = sortedX.index(firstPoint)
        yIndexOfFirstPoint = sortedY.index(firstPoint)
        for yIndex in range(yIndexOfFirstPoint -1, -1,-1):
            pointToCheck = sortedY[yIndex]
            xIndexOfPointToCheck = sortedX.index(pointToCheck)
            if xIndexOfPointToCheck > xIndexOfFirstPoint:
                nextPoint = pointToCheck
                sortedPoints.append(nextPoint)
                if xIndexOfPointToCheck + 1 == nbPoints:
                    printOnOff('Second Inflection at %s !!'%nextPoint)
                    inflectionPoints.append(nextPoint)
                    inflection += 1
                    normalInflection = True
                break
        if inflection == 2:
            break

        firstPoint = nextPoint

    if not normalInflection:
        printOnOff('Second Inflection at %s - Abnormal!!'%firstPoint)
        inflection += 1

    normalInflection = False
    firstPoint = nextPoint
    for i in range(nbPoints):
        yIndexOfFirstPoint = sortedY.index(firstPoint)
        xIndexOfFirstPoint = sortedX.index(firstPoint)
        for xIndex in range(xIndexOfFirstPoint - 1, -1,-1):
            pointToCheck = sortedX[xIndex]
            yIndexOfPointToCheck = sortedY.index(pointToCheck)
            if yIndexOfPointToCheck < yIndexOfFirstPoint:
                nextPoint = pointToCheck
                sortedPoints.append(nextPoint)
                if yIndexOfPointToCheck == 0:
                    printOnOff('Third Inflection at %s !!'%nextPoint)
                    inflectionPoints.append(nextPoint)
                    inflection += 1
                    normalInflection = True
                break
        if inflection == 3:
            break

        firstPoint = nextPoint

    if not normalInflection:
        printOnOff('Third Inflection at %s - Abnormal!!'%nextPoint)
        inflection += 1

    printOnOff('\n\n\n')
    normalInflection = False
    firstPoint = nextPoint
    for i in range(nbPoints):
        xIndexOfFirstPoint = sortedX.index(firstPoint)
        yIndexOfFirstPoint = sortedY.index(firstPoint)
        for yIndex in range(yIndexOfFirstPoint + 1, nbPoints):
            pointToCheck = sortedY[yIndex]
            xIndexOfPointToCheck = sortedX.index(pointToCheck)

            if xIndexOfPointToCheck < xIndexOfFirstPoint:
                if pointToCheck in sortedPoints:
                    printOnOff('Final Point %s added !!'%firstPoint)
                    inflection += 1
                    normalInflection = True
                    break
                #printOnOff('breaking out . . . ')
                nextPoint = pointToCheck
                sortedPoints.append(nextPoint)
                break
        if inflection == 4:
            break

        firstPoint = nextPoint

    if not normalInflection:
        printOnOff('Fourth Inflection at %s - Abnormal!!'%nextPoint)
        inflection += 1

    printOnOff('\n\n\n')
    printOnOff('No. of Inflections : %s'%inflection)

    sortedXCoords = [ xTuple[0] for xTuple in sorted(xCoords,
                                                  key=lambda xCoord: xCoord[1])]

    sortedYCoords = [ yTuple[0] for yTuple in sorted(yCoords,
                                                  key=lambda yCoord: yCoord[1])]

    printOnOff('Inflection Points : %s'%inflectionPoints)

    filterConcave = True
    if filterConcave:
        noConcaveLeft = False
        filteredPoints = sortedPoints
        while not noConcaveLeft:
            nbFilteredPoints = len(filteredPoints)
            concavePoints = []
            for i in range(nbFilteredPoints):
                if i == nbFilteredPoints - 1:
                    break

                j = i + 1
                if i == nbFilteredPoints - 2:
                    k = 0
                else:
                    k = i + 2

                indexStartPoint = filteredPoints[i]
                indexMiddlePoint = filteredPoints[j]
                indexEndPoint = filteredPoints[k]

                startPoint = (sortedXCoords[indexStartPoint], sortedYCoords[indexStartPoint])
                middlePoint = (sortedXCoords[indexMiddlePoint],
                                       sortedYCoords[indexMiddlePoint])
                endPoint = (sortedXCoords[indexEndPoint], sortedYCoords[indexEndPoint])

                if not isConvex(startPoint, middlePoint, endPoint):
                    concavePoints.append(indexMiddlePoint)

            if not len(concavePoints):
                noConcaveLeft = True

            for point in filteredPoints:
                if point in concavePoints:
                    filteredPoints.remove(point) 

        sortedPoints = filteredPoints

    makeConvexHullCurve = True
    curve = None
    if makeConvexHullCurve:
        nbBorderPoints = len(sortedPoints)

        xSortedPoints = [sortedXCoords[i] for i in sortedPoints]
        ySortedPoints = [sortedYCoords[i] for i in sortedPoints]
        zSortedPoints = [0.0 for i in range(nbBorderPoints)]

        coordsForString = (xSortedPoints, ySortedPoints, zSortedPoints)

        pStr = makeCoordString(coordsForString)
        curve = makeCurve(pStr)

    return sortedPoints, curve

def createHullFromSelectedMesh():
    mesh = Application.Selection(0)
    if not mesh or mesh.Type != 'polymsh':
        Application.LogMessage('Invalid Selection !!!\nPlease select a polymesh.', 2)
        return
    geom = mesh.ActivePrimitive.Geometry
    points = geom.Points
    posArr = points.PositionArray

    xCoords = list(posArr[0])
    yCoords = list(posArr[1])
    sortedPoints, curve = getConvexHull(xCoords, yCoords)
    curve.Name = 'convexHull%s'%mesh.FullName
    mesh.AddChild(curve)

    printOnOff('Convex Hull Created Successfully')

    return

def createHullFromRandomPoints(nbPoints = 100):
    xCoords, yCoords = createRandomPoints(nbPoints)
    sortedPoints, curve = getConvexHull(xCoords, yCoords)
    makeNulls(xCoords, yCoords, sortedPoints, curve)
    printOnOff('Convex Hull Created Successfully')

    return

#--------------------------------------------------------------------------
# If you wish test the algorithm on a random points then set the flag
#
# hullFromMesh = False
#
# else the hull will be created from the selected mesh
#
# This will create a bunch of random point and claculate the hull on these
# points. You can set the total number of points with the variable
# nbPoints.
#
# If you care for data logged out of algorithm then set the flag
# printIt = True in the top most function printOnOff. It is off by default.
#-------------------------------------------------------------------------

hullFromMesh = True

if hullFromMesh:
    createHullFromSelectedMesh()
else:
    createHullFromRandomPoints(nbPoints=100)

Enjoy!

CG · Maths · Python

Testing a point inside triangle – Math at work in cg

Here we discuss ways to derive and implement the logic for testing a point inside a triangle. In other words, we define a system for which given three points as the vertex of a triangle, we find whether a point lies inside the triangle or not. For the purpose of simplicity, we will confine these test in a 2d coordinate system.

Back to High School

Before approaching any logic/methodology we first examine the fundamentals of mathematics involved. We revisit some high school math.

Every point in a two-dimensional space can be represented as a pair of coordinates, usually represented by x and y.

point p = (x, y)

 

coordinateSystem

 

I would like to touch python implementation in parallel here as well.

This point can be represented as python tuple of two floats and can be declared as :

p = (floatA, floatB)

To define a triangle, we need three points. Let’s call them a, b, c and define them as:

a = (x1, y1)

b = (x2, y2)

c = (x3, y3)

Now we define our test point p as:

p = (x, y)

triangle

 

If we look at the point definitions above, these are nothing but 8 numbers(floats). We will do some calculations on them. Namely, addition, subtraction, multiplication and division and arrive at another number as a result of these operations. And we have our test result. It is as simple as that. Nothing more, nothing less. So now let us begin applying our logic.

Convexity

Our test requires us to examine the concept of concavity/convexity.

 

Convexity is defined as :

‘Having a surface or boundary that curves or bulges outward’

 

The word ‘concave’ has the word ‘cave’ inside it. So it has something to do with a cave. This is how I remembered it from high school.

Which of these two figures below you think represents a cave.

caves

You are right, it’s fig b.

 

This is what a concave point is :

 

In the figure below point ‘p’ is concave to points ‘a’ and ‘b’, if we move from ‘a’ to ‘b’ looking towards the left of ‘a’. Note that we define left on the side of our choice to be considered ‘inside’.

concavity01

Convex is nothing but the opposite of concave.

convexity1

Now that we have the idea what convexity is, we will exploit this concept in our logic. For our given triangle, we will check whether the point ‘p’ is concave to the three lines ‘ab’, ‘bc’ and ‘ca’.

If it is concave to all three of them then it is inside the triangle, else not.

We devised a way to check if a given point lies to the left or right of a line. Here, the right/left side of a line is determined by going from one point to another. If we go from a —-> b, then imagine a person standing on point ‘a’ and facing point ‘b’. All the point lying on the right side of this person are right and all points on the left side of this person are left.

leftRight

Two Point Form Equation of  a Line

To do this, we make use of two-point-form of equation of a line:

y = y1 + [(x – x1) (y2 – y1) / (x2 – x1)]

 

Again going to python implementation, this could be easily written as

y = y1 + (   (x – x1) * (y2 – y1) / (x2 – x1) )

In the above the line is defined between the point (x1, y1) and point (x2, y2). The general point (x, y) represents any point on this line.

line

Now that we have a solid formula to use, let’s apply this to our case. Given the three sides of the triangle ‘ab’, ‘bc’ and ‘ca’, we test whether our test-point point lies to the left or right to each of these lines.

To find whether the point lies to left or right, we take the x coordinate of point ‘p’ and put in it the two-point-form equation of the line.

This will give us the y coordinate of the point ‘p’ if it were on the line.

Next, we check whether this y coordinate that we calculated is greater or smaller than the actual y coordinate of ‘p’. If it is greater than ‘p’ is on the left else it is on the right.

twoPointLineUse

Let the points forming the line be ‘a’ and ‘b’ defined as:

a = (x1, y1)

b= (x2, y2)

and point ‘p’ as:

p = (x3,  y3)

 

Next we find the y coordinate if were to lie on ‘ab’,

this will given by:

y = y1 + [(x3 – x1) (y2 – y1) / (x2 – x1)]

Let’s take a simple view of the above equation. In the above equation, we have five number on the right. x1, x2, x3, y1 and y2. We are doing simple operations of  +, – , x and / on them and we get a sixth number y. Now we just want to calculate whether y is greater than y3 or not. Simple, and we get the answer if lies on the left or right.

Also, we want to be sure that we take all the cases in going from ‘a’ to ‘b’. Point ‘b’ could be in any quadrant relative to ‘a’. There are exactly four possibilities for ‘b’ lying in each quadrant.

quadrants

If ‘b’ lies in quadrant ’0′ or ’3′ then we note that:

y > y3, ‘p’ is on the right of ‘ab’

y < y3, ‘p’ is on  the left of ‘ab’

 

if ‘b’ is in quadrant ’1′ or ’2′ then the conditions are reversed.

y > y3, ‘p’ is on the left of ‘ab’

y < y3, ‘p’ is on right of ‘ab’

 

So we have to check for ‘b’ position with respect to ‘a’. And then use that information along the comparisons to finally get the result.

Let us put the python implementation for this. This code is more verbose than what it should be. It could definitely be more concise. Parts of it are also redundant. But I wanted this to be more readable and hence the verbosity.

 

def isConvex(pointA, pointP, pointB):
    x1, y1 = pointA

    x2, y2 = pointB

    x3, y3 = pointP

    if x2 > x1 and y2 > y1:
        quadB = 0
    elif x2 < x1 and y2 > y1:
        quadB = 1
    elif x2 < x1 and y2 < y1: quadB = 2 elif x2 > x1 and y2 < y1: quadB = 3 y = y1 + ((x3 - x1) * (y2 - y1) / (x2 - x1)) if quadB == 0 or quadB == 3: if y > y3:
            return True
        else:
            return False
    elif quadB == 1 or quadB == 2:
        if y > y3:
            return False
        else:
            return True

Following is the example of the above function at work.

>>a = (1.0, 6.0)
>>b = (3.0, 8.0)
>>p = (2.0, 12.0)
>>isConvex(a, p, b)
>>False

Now that we have a working function to test for convexity(‘right’). It is very easy to get all the three points of a triangle – ‘a’, ‘b’, ‘c’ and test a point – ‘p’ is inside or outside.

Let write the code for that:

def insideTriangle(a, b, c, p):
    if (not isConvex(a, p, b)) and (not isConvex(b, p, c)) and (not isConvex(c, p ,a)):
        return True
    else:
        return False

Please note that we are using the function isConvex() which we defined earlier.

And that’s it, we have our function to check whether the point lies in the triangle or not.

 

Now that  we know the concept of our test inside-out, let’s put it to a simple test. Please download the attached scene file and the script. There is a null called ‘p’ in the scene, that we test inside the triangle defined by three other nulls ‘a’, ‘b’ and ‘c’.

 

Here’s the final script to run the test inside the sample scene file

def getObjectByUserKeyword(inKeyword):
    oColl = Application.ActiveSceneRoot.FindChildren('')

    for o in oColl:
        keywordColl = Application.GetUserKeyword(o)
        for thisKeyword in keywordColl:
            if thisKeyword == inKeyword:
                return o

    return None

def getScenePoints():
    a = getObjectByUserKeyword('a')
    b = getObjectByUserKeyword('b')
    c = getObjectByUserKeyword('c')
    p = getObjectByUserKeyword('p')

    return a, b, c, p

def getPointCoordinates(inPoint):
    x = inPoint.Kinematics.Global.Parameters('PosX').Value
    y = inPoint.Kinematics.Global.Parameters('PosY').Value

    return x, y

def getPointCoordinateList(inList):
    coordinateList = []
    for point in inList:
        coordinateList.append(getPointCoordinates(point))

    return coordinateList

def setVisibility(on=True):
    inside = getObjectByUserKeyword('inside')
    outside = getObjectByUserKeyword('outside')
    inside.Properties('Visibility').Parameters('viewvis').Value = on
    outside.Properties('Visibility').Parameters('viewvis').Value = not on

    return

def isConvex(pointA, pointP, pointB):
    x1, y1 = pointA

    x2, y2 = pointB

    x3, y3 = pointP

    if x2 > x1 and y2 > y1:
        quadB = 0
    elif x2 < x1 and y2 > y1:
        quadB = 1
    elif x2 < x1 and y2 < y1: quadB = 2 elif x2 > x1 and y2 < y1: quadB = 3 y = y1 + ((x3 - x1) * (y2 - y1) / (x2 - x1)) if quadB == 0 or quadB == 3: if y > y3:
            return True
        else:
            return False
    elif quadB == 1 or quadB == 2:
        if y > y3:
            return False
        else:
            return True
def insideTriangle(a, b, c, p):
    if (not isConvex(a, p, b)) and (not isConvex(b, p, c)) and (not isConvex(c, p ,a)):
        return True
    else:
        return False

scenePoints = getScenePoints()

a, b, c, p = getPointCoordinateList(scenePoints)

testValue = insideTriangle(a, b, c, p)

setVisibility(testValue)

To run the test:

1. Move the null ‘p’ , or any of the nulls ‘a’, ‘b’ and/or ‘c’  to make ‘p’ lie inside or outside the triangle in the front view. If you move the points ‘a’ , ‘b’ or ‘c’, please make sure that cyclic order ‘abc’ is maintained in the anticlockwise direction. This is important so that we always test for points on the left of the line.

2. Run the script.

3. The text-curves ‘inside’ and ‘outside’ will be visible according to the result of the test.

 

Note that, as we are testing in a 2D plane , we use the only the front view as our coordinate plane and the points have their position z = 0.

Enjoy !!

 

Download insideTriangleTest.scn