Generate Polygon from points

[cp51]

Hi all,

This question is regarding creating a custom polygon from a set of semi-ordered points. Let me try and explain this clearly:

The Problem:
I have a set of points (xyz coords) that define the outline of a shape. Lets pretend the shape looks like a rain drop. The Z axis goes through the raindrop from bottom to top(straight through the center of the rain drop when the drop is falling)

Now, if you think about slicing the drop with the x-y plane for various z depths, you end up with a bunch of circles, each with a different radius, large towards the bottom of the drop and getting smaller towards the top. I have the points grouped into these circles by z-depth, and points are ordered so that they come one after another around the circle(they are not just randomly placed around the cirle, they are in order)

Knowing all this, I would like to connect all the points together with triangles to create a polygon.

If that doesnt make any sense:
In simplest terms, I have the hull of an object that I would like to make into a polygon. I dont know if the semi-ordering would be of any use.

Does anyone know how I might go about generating the polygon? I've looked at some kind of convex hull algorithm maybe, but im having trouble getting information about doing that in 3D. Eventually I would like to make these models and draw tem in irrlicht, but I dont know the order to connect all the points together.

If anyone has any tips or could point me in some direction, i would greatly appreciate it.

Thank you.

[chronologicaldot]

I've been interested in doing the same thing for awhile, and I intend to get around to doing it. My method involves using an iterator class that cycles when it goes through a certain number of points, but you could also use an integer.

For example, suppose for your polygon, you have three sets of points. We'll keep the example simple: each set has the same number of points:
1, 2, 3, 4
5, 6, 7, 8
9, 10, 11, 12

You need to make the following sets of three points for your faces:
(1,2,5) ; (2,5,6) ; (2,3,6) ; (3,6,7) ...
or if we label the bottom row "z" and the top row "w":
(z,z+1,w) and (z+1,w,w+1)
Hence, you really only need two points and then you just repeat this for every row in your axis minus 1 (since you won't need to do this for the top row). When moving on to the next row, you set z to cycle over the values w once was at and then w gets to cycle over the next range of values. For an integer, that means going from row 1 to row 2, you add 4 (the size of row 1).

Hopefully you can imagine at this point what the role of the iterator would be and why you would need to cycle it. The cycling is merely so I don't have to worry about resetting the value back to the start.

Once you reach the top, if you want to close your model (so it isn't a cup or hollow barrel), you make/set a center point and repeat the process but only creating triangles of the type (z,z+1,w) since there is only one point at the top.

[cp51]

Thanks for the advice! I was thinking I may have to do something similar to this, but I didnt know if there was a simpler method.

I guess the biggest problem is going to be that not all layers have the same number of points, and the points dont necessarily all start around the same location. So while they all follow the circle around in order, the ordering may start with a point at 3 oclock for the bottom layer, and the next layer up could start at 6 oclock around the circle, though I could probably just do a nearest neighbor search to get a good starting point.

To correct for not having an equal number of points, I figured I could just interpolate points for each layer until I have the same number everywhere.

I guess Ill give it a shot and see what happens.

Thanks again!

[chronologicaldot]

Hm... reducing the number of points makes things a bit more complicated. I can think of ways to do it (none I've worked out on paper yet). You can do that. You could also change the step size for the integer or iterator (whatever you're using). The first important thing, at least in my mind, is to make sure the first point of each layer is directly below the first point of the layer above it - it makes for an easy start at least.

I'm interested in how you are generating your points. Are you creating them from a picture? From your own built-in line/curve system? Or just automatic circles of varying sizes?

[cp51]

Basically, I have a bunch of images. Each image is a picture of an object that has been sliced at various depths so you just end up with a bunch of cross section images.

I do some image processing to find the object in each image, then use some active contours to pick out outline of the shape.

Depending on the size of the object in the image the contour can have a varying number of points(more points when the cross section is largest)

After all that I end up with points defining the outline of the object, and each image is a slice at a different z-depth.

To end up with the same number of points in each layer I wrote a little python script that takes the contour and interpolates some set number of points:
I just wrote this quickly yesterday so there are probably much more efficient ways of coding this up, but I just wanted to make sure it would work firrst. The process I did was:
1) find center of mass of contour and shift everything so it is centered at (0,0)
2) calculate the angle to each point from the origin
3) Sort the set of points based on their angle(so they are ordered based on the angle from 0-2pi)
4) Interpolate
5) Generate some number of points back
6) Repeat for each contour

Ill point out that this works for convex contours, or atleast contours where the center of mass is inside the shape. if it gets too concave then the interpolation doesnt work because the com is outside of the shape, and that messes things up.
I just discard those layers, I have enough without them so its not a big deal to me.

Also, I should also say that each image has multiple objects, so on top of this I have to find out which contours belong to the same object in each slice. I do a dbscan clustering of the center of mass of each object to group them all together. thats what the label is for in the code below.(not really relevant to the process, but just incase it was confusing)

Code: Select all

 
def do_interp(cx, cy,z,l):                               #cx,cy are my contour point array, z is my depth, l is a label(can ignor z and l for this)
        #Generate angles to each point
        comx = mean(cx)
        comy = mean(cy)
        #Calc center of mass
        cx = array(cx) - comx
        cy = array(cy) - comy
        angle = arctan2(cy,cx) + pi
 
        #Sort points by angle
        comb = zip(angle,cx)
        comb.sort()
        angle_sorted, cx_sorted = zip(*comb)
 
 
        comb = zip(angle,cy)
        comb.sort()
        angle_sorted, cy_sorted = zip(*comb)
 
        #For interpolate function, angle must be strictly increasing. This causes a problem if for some reason the angle is the same for 2 points, so I just delete if not strictly increasing
        angle = []
        cx = []
        cy = []
        for i in range(1,len(angle_sorted)):
                if angle_sorted[i-1] < angle_sorted[i]:
                        angle.append(angle_sorted[i-1])
                        cx.append(cx_sorted[i-1])
                        cy.append(cy_sorted[i-1])
        if angle_sorted[-1] > angle_sorted[-2]:
                angle.append(angle_sorted[-1])
                cx.append(cx_sorted[-1])
                cy.append(cy_sorted[-1])
 
        #Probably not needed, loop back around to fix end points not matching up
        for i in range(len(angle)):
                angle.append(angle[i]+2*pi)
                cx.append(cx[i])
                cy.append(cy[i])
        #Scipy interpolate function, nice and easy
        s0x = UnivariateSpline(angle,cx,s=0)
        s0y = UnivariateSpline(angle,cy,s=0)
 
        #Angle to recalculate points at
        angle = arange(0,2*pi,0.5)
        #images are 600x600, make sure we didnt accidentally interpolate points way outside the image
        #If contour is concave and the center of mass is outside of the shape, then we have problems. Just disscard those, I dont need them
        s0xs = s0x(angle) + comx
       for i in range(len(s0xs)):
                if s0xs[i] > 600 or s0xs[i] < 0:
                        return 0,0
                     
 
        for i in range(len(s0ys)):
                if s0ys[i] > 600 or s0ys[i] < 0:
                        return 0,0
                       
        #Print out the values
        for x,y in zip(s0xs, s0ys):
                print x,y,z,l
        print 'x'
        return s0xs, s0ys
 
 

So this should give me the same number of points for each contour, all starting at the same angle with respect to the x axis, so they should all be ontop of eachother.

I just generated all 68 objects in the image stack into polygons (basically just the 3 corners of each triangle) and im gonna try drawing them now, but I suspect ill have a bit of tweaking to do.

[chronologicaldot]

Looks good for the most part. Personally, I wouldn't try to put multiple object on the same image because it makes it more difficult to tell which is which when you are combining layers. What you'll find is that the different objects start connecting with each other with seemingly random triangles.

But dang - you're work is really getting me motivated to work on this. XD

In order to avoid the center of gravity issue and worrying about concavity and such, I plan on walking-around-the-edge so-to-speak when generating my points, and then rearranging my points such that the first one I grab is from the same general location in the image and almost directly above the starting point on the previous layer (the distance being determined by pixel location).

Let us know how your images turn out.

[cp51]

So... I always appreciate it when someone gives a example of wat they are talking about, so here we go

I am going to generate a "raindrop" looking shape.

The radius as a fn of Z will vary like this:


The center of mass of the drop will vary along the z axis like this:


So you could imagine, a raindrop, fat at the bottom, skinny on the top, and kind of bending over. Almost looks like a golf club when its all done...

To Generate the points that make up the drop, I generate a circle of radius r for that Z layer based on my radius fn plotted above.
Then I shift the points by the center of mass(along x axis only) based on my COM function above:

Code: Select all

#Generate Raindrop given R and COM
def raindrop(radius,com):
    layer = []
    for r,c in zip(radius,com):
        #Number of points generated will change based
        #on the circumference of this layer
        angle = linspace(0,2*pi,2*pi*r)
 
        cx = r*cos(angle) + c #Shift everything by center of mass
        cy = r*sin(angle)
 
        #Save Our points for each layer:
        #layer[z][x/y]
        layer.append([cx,cy])
 
    return layer

Rendering all these points looks like this:


Notice that each ring doesnt have the same amount of points, so we wont be able to connect them all together by triangles easily.
That was all just generating data, assuming you obtained your points elsewhere(curves or something) you would start here.

So now, I take each ring and interpolate the same amount of points at each z layer:

Code: Select all

#Given Our coordinates, interpolate
def do_interp(cx,cy):
    #Center this ring
    comx = mean(cx)
    comy = mean(cy)
    
    cx = array(cx) - comx
    cy = array(cy) - comy
 
    #Get the angle to each point
    angle = arctan2(cy,cx) + pi
 
    #Sort points by angle so angle is always increasing 0-2pi
    angpt = zip(angle,cx)
    angpt.sort()
    angle_sorted, cx_sorted = zip(*angpt)
 
    angpt = zip(angle,cy)
    angpt.sort()
    angle_sorted, cy_sorted = zip(*angpt)
 
    #For some reason, im still getting angles that arent strictly increasing
    #Even though Everything is calculated 0-2pi
    #So here I just check that it is all strictly increaseing
    #And just orget the angles that arent
    angle = []
    cx = []
    cy = []
    for i in range(len(angle_sorted)-1):
        if angle_sorted[i] < angle_sorted[i+1]:
            angle.append(angle_sorted[i])
            cx.append(cx_sorted[i])
            cy.append(cy_sorted[i])
    angle.append(angle_sorted[-1])
    cx.append(cx_sorted[-1])
    cy.append(cy_sorted[-1])
    
    #Interpolate
    s0x = UnivariateSpline(angle,cx,s=0)
    s0y = UnivariateSpline(angle,cy,s=0)
 
    angle = linspace(0,2*pi,10) #10 pts per layer
 
    s0xs = s0x(angle) + comx
    s0ys = s0y(angle) + comy
 
    return s0xs, s0ys

This ends up looking like:


Now that I have the same number of points per layer, I can connect the points together with triangles:

Code: Select all

def makePolys(layer):
    #Loop over each layer in layer array
    #Z location is index of each layer
    triangles = []  #[ [pt1],[pt2],[pt3],   [pt1],[pt2],[pt3], ...]
    for z in range(len(layer) - 1):     #Leave out last Layer for now
        for i in range(len(layer[z][0]) - 1): #length of x coord array(should be 10)
            pt1 = [layer[z][0][i], layer[z][1][i], z]
            pt2 = [layer[z][0][i+1], layer[z][1][i+1], z]
            pt3 = [layer[z+1][0][i], layer[z+1][1][i], z+1]
            triangles.append(pt1)
            triangles.append(pt2)
            triangles.append(pt3)
 
            pt1 = [layer[z][0][i+1], layer[z][1][i+1], z]
            pt2 = [layer[z+1][0][i], layer[z+1][1][i], z+1]
            pt3 = [layer[z+1][0][i+1], layer[z+1][1][i+1], z+1]
            triangles.append(pt1)
            triangles.append(pt2)
            triangles.append(pt3)
 
    return triangles

And finally Render it:

Code: Select all

def render(triangles):
    f = faces(pos=triangles)
    f.make_twosided()
    return f

And we end up with this:



And here is the code to set this all up and call each function:

Code: Select all

r = sin(linspace(0.1,pi-.6,20))
r = list(r[0:-1])
for i in linspace(2,15,20):
    r.append(2*sin(pi-.6)*1/i)
r = 20*array(r)
figure(1)
title('Changing Radius')
ylabel('radius')
xlabel('z')
plot(r)
 
x = linspace(0,15,len(r))
com = (x-5)*(x-5)
figure(2)
title('Center of Mass (x axis)')
ylabel('com_x')
xlabel('z')
plot(com)
 
layer = raindrop(r, com)
 
 
interp_layer = []
for z in range(len(layer)):
    cx,cy = do_interp(layer[z][0], layer[z][1])
    interp_layer.append([cx,cy])
 
for z in range(len(interp_layer)):
    for i in range(len(interp_layer[z][0])):
        sphere(pos = (interp_layer[z][0][i], interp_layer[z][1][i], z))
triangles = makePolys(interp_layer)
render(triangles)

The first part is just generating my r and com functions and plotting them, the real fun starts when I call layer = raindrop(r,com)

All the rendering I did in visual python. It is just faster to set up than doing it all in irrlicht. However, moving the camera is a pain, the lighting doesnt like to cooperate, documentation is horrible and python is slow. So eventually I will render everything using irrlicht. I just like prototyping in python since its fast and I dont waste a ton of time on an idea that doesnt work.

[chronologicaldot]

You're camera angles slightly confused me for a while there, but I think I get it. Looks good to me... at least from what I can tell is going on.
Nice job

(I like Python also )