Add a primitive surface area function ... for polyhedron with 4 to 8 sides (ARB8)BRL-CAD
Status: ClosedTime to complete: 72 hrs Mentors: Sean, Matt S.Tags: C, C++, math, geometry, surface area

BRL-CAD provides more than two dozen types of geometry "primitives" such as ellipsoids, boxes, and cones. Every primitive is described by a collection of callback functions, for example rt_ell_bbox() returns the bounding box dimensions for an ellipsoid. Wikipedia, Wolfram Mathworld, and various other math sites (and research papers) around the web include the equations for most of our basic primitives while others are a little more tricky to compute.

This task involves writing a new callback function that takes an rt_db_internal object and calculates the surface area (units are mm^2). There are numerous examples in our code where we compute surface area for other primitives. The primitives that do not already have a centroid callback are itemized in following.

References:

Code:

  • src/librt/primitives/arb8/arb8.c
Uploaded Work
File name/URLFile sizeDate submitted
area.patch3.6 KBNovember 30 2012 21:14 UTC
areasimple.patch2.0 KBDecember 02 2012 16:01 UTC
areamatrix.patch3.7 KBDecember 02 2012 18:06 UTC
Comments
Silvrouson November 28 2012 04:46 UTCTask Claimed

I would like to work on this task.

Harmanpreet Singh on November 28 2012 04:47 UTCTask Assigned

This task has been assigned to Silvrous. You have 72 hours to complete this task, good luck!

Silvrouson November 28 2012 16:57 UTCextracting points

I'm trying to see how the points are accessed. In rt_arb_volume, for instance, the rt_db_internal  *ip-idb_ptr is passed to a rt_arb_internal, which stores the points. But the idb_ptr is a pointer. To what does it point?

Silvrouson November 28 2012 17:07 UTCnr points

Also, a related question, how are ARB's 4-7 handled, since the general definition is 8-vertexed? The primitives docs only state that they are handled by arb8

Sean on November 28 2012 21:48 UTCpoints are replicated

For performance and simplicity, the points are just duplicated in a specific order.  There are functions that will tell you whether a given arb8 is a 4-7 and from there you can just read the first N points iirc.


Understanding idb_ptr is getting very low-level but basically idb_ptr IS an rt_arb_internal.  It's not passed, it's cast.  The generic idb_ptr pointer is converted into a struct rt_arb_internal pointer. 

Silvrouson November 30 2012 13:41 UTCone more thing

 I think the arb8 volume function might be wrong. For calculating the base area, it uses the cross product of (OB-OA)X(OC-OA) = ABXAC, but that gives twice the actual area of the triangle ABC, iirc. Should I halve the result of that calculation when calculating the area?

Matt S. on November 30 2012 21:00 UTCElaboration

What cross product, and I guess more specifically, what function are you referring to?

Silvrouson November 30 2012 21:14 UTCReady for review

The work on this task is ready to be reviewed.

Silvrouson November 30 2012 21:18 UTCrt_arb_volume

In the RT_ARB_VOLUME function, on line 2309, the function computes the area of the base of the arb4 to get the volume, by using vector cross products, but the formula is 1/2*(cross_product).

Matt S. on November 30 2012 21:52 UTCNo, it's correct

Note that the area of a tetrahedron is:


V = \frac{1}{3} A_0 h


which follows from the volume formula for a general polyhedra.  Notice that in rt_arb_volume we have:


        /* calculate area of arb4 base */


        VSUB2(b_a, aip-pt[b], aip-pt[a]);


        VSUB2(c_a, aip-pt[c], aip-pt[a]);


        VCROSS(area, b_a, c_a);


 


        *vol += MAGNITUDE(area) * arb4_height;


    }


    *vol /= 6.0;


so the factor of 1/2 has just been factored out.

Silvrouson November 30 2012 22:02 UTCOh, right

My bad, I should have seen that. Especially since I did something similar in my own function.


 

Matt S. on November 30 2012 22:19 UTCOK, but...

Since you have the vertes coordinates, could you not simply compute the area via Green's Theorem?

Silvrouson November 30 2012 22:35 UTCI suppose

I know the 2D version of green's theorem, but I thought this way was easier to implement

Matt S. on November 30 2012 22:47 UTCBut it is 2D!

Each face is planar, and thus the so-called 2D version of Green's Theorem applies directly. 


As an algorithm, one could get the vertex coordinates of a single face (which is always a triangle in this case) and compute an affine transformation that would translate the face in a manner such  that one vertex is at the origin, and rotate it so that it is on the x-y plane.  This would allow the aree of each face to then be computed rather trivially via Green's Theorem.


The (possible) advantage of this approach is that:


A; it could be extended to compute the area of a BOT primitave fairly easily, although counting vertices / faces becomes a non-trivial issue


B: It can be extended VERY easily to the arbn primitave as well, assuming such a function does not exist (I have not looked to be honest).


Just my thoughts of course, feel free to ignore my ramblings :-)

Silvrouson November 30 2012 23:16 UTCWell,

That was actually my first idea, but I didn't know how to do the transformation.

Matt S. on November 30 2012 23:33 UTCThat I can help with

Give me a few minutes, I have enough material kicking around here somewhere on just this topic...

Matt S. on December 1 2012 01:11 UTCMethod Attached

Please refer to the attached .pdf file, I've outlined the math required to implement the method I outlined above.


-Matt

Matt S. on December 1 2012 01:14 UTCAck!

Can I not attach documents?!?

Sean on December 1 2012 01:40 UTCyou should

You should be able to.  Does it not give you that option when it's in a needing review state?  Someone else reported that they still could.

Sean on December 1 2012 01:43 UTCOH

My misunderstood who was askin.  Mentors can't upload docs and can't even modify the task once it's claimed to make edits.  Best to just post a URL here in the comments.  If you send me the PDF and I can it up on the server or you can add it to the wiki (accepts pdf files).


 

Sean on December 1 2012 01:45 UTCDeadline extended

The deadline of the task has been extended with 2 days and 0 hours.

Matt S. on December 1 2012 01:47 UTCSent

OK, I've sent it to you via the dev mailing list.


Thanks!

Matt S. on December 1 2012 02:15 UTCBackground Math

For real this time!


http://brlcad.org/wiki/Image:Affine_transformations.pdf

Silvrouson December 1 2012 09:42 UTCthanks

OK, so I take each face, translate and rotate its points to a plane and calculate the surface using the derived green's formula, abs(1/2 sum(xi yi+1 - xi+1 yi))


I have 2 questions:


-When doing the transformations, wouldn't it be easier to just subtract the coordinates of x0 from the other points on the face to translate the face to the origin, instead of doing matrix transformations?


-What does the exponent notation do on a point? For instance, like it's done in the b vector in 1.1, with [x^1, x^2 x^3]

Silvrouson December 1 2012 09:47 UTCmatrix

Also, where are the matrix operations defined?

Matt S. on December 1 2012 12:48 UTCTwo answers

Not sure where you get abs() from, I was thinking using the Green's Theorem in a manner like this guy has laid out:


http://shuisman.com/?p=359


where he has said that


 A = \frac 12 \sum \limits_{n=1}^{N} det \begin{pmatrix} x_n & y_n \\ x_{n+1} & y_{n+1} \end{pmatrix}


is this what you meant?  Anyway, as for your questions...


First, what you have said is exactly what I have suggested, so we're on the same page here.  Building the augmented matrix A_T that I have suggested simply keeps things compact and easily managed--the algorithm can be implemented in a single line of code in numpy for example.


Second, what you're calling an exponent is simply a superscript.  This is pretty common notaiton within the realm of differential geometry, and you can read up on if you like--it's really useful stuff!  For now though, all it means is that instead of having a frame in R^3 represented as a triple of numbers in the x-y-z directions, think of them as a triple of numbers in the 1-2-3 directions.  In this way, the letter x in x^1, x^2, x^3 represents a coordinate system x, and the numbers represent the coordinate values.


This is more flexible than x-y-z notation because we can have as many dimensions as we like (countable, of course) without the problem of running out of letters that we do with x-y-z.  As well, the letter being used can also represent the coordinate system employed.  That is, we could say that


(y^1, y^2, y^3) = T(x^1, x^2 x^3)


and not have to worry about running out of letters or numbers, but saying the same thing in x-y-z terms we would have to do something like


(u, v, w) = T(x, y, z)


which can get surprisingly confusing quite quickly.


That make sense?

Matt S. on December 1 2012 12:52 UTCMatrix operations?

What do you mean where are they defined?  You simply construct them as I've described, thus defining them.  Or am I not understanding your question?

Silvrouson December 1 2012 13:44 UTCmatrix

About the formula, yes, that's what I meant, I used the abs so it wouldn't matter in which order the points were traversed.


I was referring to the functions/classes in the brl-cad libraries I could use to compute matrices, something similar to vect_t for vectors, and VCROSS for cross products. I don't know where to find them, I didn't see any in vmath.h.


About the superscripts, I get it, you're naming the x y z axes to be 1 2 3, but then what is the coordinate system x used in b? Also, why is b defined as minus the vector, and then used as -b again in the matrix A_T, which factors out the minus?


 

Silvrouson December 1 2012 14:23 UTCaah

Is it the set of unit vectors (100) (010) (001) ?

Matt S. on December 1 2012 14:35 UTCIt's a typo

It should be b = [x_0^1, x_0^2, x_0^3] and then have a "-b" term in the matrix A_T.


As for the coordinate system x, I'm simply taking that as the original coordinate system.  That is, if you have a point x_0=[5, 6, -8]  -- as well as bunch of other points x_i -- then you translate them all to a new set of values y_i  via:


y1iy2iy3i1=1000010000105681x1ix2ix3i1


As for defining a matrix, isn't that just via mat_t ?

Matt S. on December 1 2012 14:37 UTCJumbled mess!

Sorry about that, it got magled... What was there was simply:


\begin{equation*}


\left[\begin{array}{c}


    y_i^1\\


    y_i^2\\


    y_i^3\\


    1\end{array}\right] =


\left[\begin{array}{cccc}


    1  0  0  -5\\


    0  1  0  -6\\


    0  0  1  8\\


    0  0  0  1\end{array}\right]


\left[\begin{array}{c}


    x_i^1\\


    x_i^2\\


    x_i^3\\


    1\end{array}\right]


\end{equation*} 

Silvrouson December 1 2012 14:43 UTCIt's still jumbled

And the mat_t is defined as a fastf_t 4x4, I was thinking maybe there were some quicker, numpy-like operations I could perform...

Matt S. on December 1 2012 15:29 UTCNo, it's just in TeX format now...

As for defining the matrices, I see what you're asking now.  Unfortunately, I don't have a quick answer for you...  Could you not simply define the transformation(s) yourself though?


 

Silvrouson December 1 2012 15:52 UTCDon't know

I guess I could, but it could be quite involved...


Maybe you could review the function as it is now, and I'll generalize the code for the other primitives after the competition, when I've got more experience?

Sean on December 2 2012 03:00 UTClots of math routines

There are lots of matrix macros in include/vmath.h and other matrix functions throughout the libbn library (see include/bn.h or src/libbn/*.c).


As for the patch itself, be sure to read it with a text editor before submitting.  You'll notice that it includes a unintentional change to the INSTALL documentation file.


 

Sean on December 2 2012 03:00 UTCTask Needs More Work

One of the mentors has sent this task back for more work. Talk to the mentor(s) assigned to this task to satisfy the requirements needed to complete this task, submit your work again and mark the task as complete once you re-submit your work.

Silvrouson December 2 2012 13:08 UTCgot it

Those were just what I was looking for, I think I've got it now.


As for the INSTALL, I have no idea why the change appeared, all I did was compile and then delete the .build folder.

Silvrouson December 2 2012 14:00 UTCu^1^2

In the u miu X u miu matrix, what does (u^1)^2 and (u^2)^2 mean?

Silvrouson December 2 2012 14:27 UTCTensor product!

Disregard that, I've just figured it out.

Silvrouson December 2 2012 14:29 UTCtypo?

In the tensor product, should that be u or umiu?

Silvrouson December 2 2012 18:05 UTCfinally

I assumed it was umiu, the other possibility didn't make sense since umiu wouldn't have been used anywhere.


areasimple.patch is my old function,now fixed, and areamatrix.patch is the new one, using affine transformations.

Silvrouson December 2 2012 18:07 UTCReady for review

The work on this task is ready to be reviewed.

Matt S. on December 3 2012 00:55 UTCThis looks good!

And yes, another typo on my part--it should be u_{\mu} in the computation of the matrix R as you have figured out.  Thanks for spotting that!


As for the code, it looks good.  I have not had time to go over it in detail, nor will I today I don't think.  I should be able to get to it this evening.  In the meantime, should I extend the deadline again to avoid having the clock run out, or is that an issue now that you have submitted some work?  I'm asking because I don't actually know how this system works...


Cheers!


(Oh, and I'll also need to figure how to actually test this--how have you managed to run your code thus far?)

Sean on December 3 2012 02:18 UTCmore tasks for testing

We'll add some more tasks to update the 'analyze' command so they call these new volume, surface area, and centroid functions.


We can extend the deadline arbitrarily on our end, so just ask if you need it but I think you've put enough work in to at least deserve credit for the task.  We can create another to "fix it" if there's a problem.


  

Sean on December 3 2012 02:18 UTCTask Closed

Congratulations, this task has been completed successfully.

Sean on December 18 2012 05:29 UTCfollow-on task

A follow-on task has been posted:


http://www.google-melange.com/gci/task/view/google/gci2012/8098203