Implementing a sextic equation for the heart
I intend to implement the equation of the heart in C which is talked about on Valentine's day. The sextic ( of power 6 ) equation associated with the Three dimensional heart is
( x^2 + 9/4*y^2 + z^2 -1 )^3 - x^2*y^3 - 9/80*y^3z^3 = 0.
Please any ideas would be appreciated.
I am working on a BRLCAD project - a computer simulation software which helps visualise products during design in military and manufacturing. I want to build a heart primitive for this equation and need ideas on how to do this or references on algorithms and papers .
What do you mean by "implement"?
My homepage
Advice: Take only as directed - If symptoms persist, please see your debugger
Linus Torvalds: "But it clearly is the only right way. The fact that everybody else does it some other way only means that they are wrong"
Unfortunately in CAD systems things aren't as simple as making a solid from an equation. I found this out myself when I wanted to make the outline for a Wankel engine in AutoCAD. You'll have to create a mesh or surface composed of a bunch of points from that equation. This is a little difficult because you want mesh points to be evenly spaced so that the created mesh appears smooth. You won't have just a step size or something like that.
You'll need to look for papers on "meshing" 3D solids, most likely.
Epy, no. BRL-CAD uses constructive solid geometry, using primitive solids and their unions, intersections, and exclusions to define the object. It has very limited support for triangular meshes.
Nyah Check, if you only wanted to construct a heart, you could easily do so using ellipsoids and cuboids (cuboids intersecting with the ellipsoids to make it easier to join the parts seamlessly, and build the heart from easier building blocks) -- perhaps just two ellipsoids, one for each half, with a cuboid intersecting with that half, and joined at x=0 plane. (That's what I'd do. I don't like the nipple at the bottom of the heart at all.)
However, I am assuming you want to see how you could add the heart as a primitive solid to the program.
The main work is pretty much math. I didn't bother to check exactly how BRL-CAD does it, so I'm offering the points below as they relate to constructive solid geometry in general.
The equation
(x2 + 9/4 y2 + z2 - 1)3 - x2 y3 - 9/80 y3 z3 = 0defines the surface of the heart. It happens to be symmetric along the x axis, so you can rewrite it into x=±(function(y,z)), will probably turn out to be very useful. (Solve the equation for x2, and you get a cubic equation you can solve algebraically.)
For lighting (ray tracing), you'll probably need to parametrize the surface equation, converting it into form
x = function1(u, v)so that you can compute the normal vector using the cross product of the partial derivatives of the coordinate functions. That is,
y = function2(u, v)
z = function3(u, v)
xu(u,v) = ∂function1(u, v) / ∂uwhere vectors
xv(u, v) = ∂function1(u, v) / ∂v
yu(u, v) = ∂function2(u, v) / ∂u
yv(u, v) = ∂function2(u, v) / ∂v
zu(u, v) = ∂function3(u, v) / ∂u
zv(u, v) = ∂function3(u, v) / ∂vare tangential to the surface; therefore
tu(u, v) = <xu(u,v) yu(u,v) zu(u,v)>
tv(u, v) = <xv(u,v) yv(u,v) zv(u,v)>
p(u, v) = tu(u,v) × tv(u,v)gives a vector perpendicular to the tangent, i.e. normal to the surface.
normal(u, v) = p(u, v) / || p(u, v) ||
The most important formula is the inequality that defines which part of space is inside the heart. This is given by the inequality
(x2 + 9/4 y2 + z2 - 1)3 - x2 y3 - 9/80 y3 z3 < 0
For visualization (rendering), you will absolutely need a line-heart intersection test. You can optimize the test by realizing that any line (ray) that comes closer than 0.61 of the center of the heart will intersect the heart, and only those lines that come closer than 1.1 of the center of the heart may intersect it. (The point of intersection(s) are those that satisfy the line equation and the equation for the surface of the heart at the same coordinates x, y, z.)
For lighting, you'll also need the surface normal at those intersection points; I already explained the common/typical way to find the normal, above.
Depending on the BRL-CAD innards (the kind of processing it supports on analyzing the composite objects modelled), you may also need more complex intersections, probably between the heart and each of the primitive solids supported by the program. These are solved the same way you solve the heart-line intersection, but the solution is a surface in 3D (usually parametrized using u,v).
Some links:
- Constructive solid geometry at Wikipedia
- Line-sphere intersection at Wikipedia
- Surface normal at Wikipedia
- Normal vector at Wolfram MathWorld
- 3D point-line distance at Wolfram MathWorld (for line-sphere intersections!)
- My user page at Wikipedia, describing line-cylinder intersection, with the line passing through origin. (I did that page in 2010, when I needed efficient solutions for line-sphere and line-cylinder intersections for ray-tracing ball-and-stick molecular models with millions of atoms. The trick with rays passing through origin (at optical focal point) turns out to make everything much simpler and faster to calculate.)
When you have the above math roughly solved, and written them as sufficiently efficient C code (to see there are no precision/numeric issues), you can start adding the primitives to BRL-CAD. Note that unless you can do all the math that BRL-CAD needs for the heart, none of it is useful; it's an all-or-nothing deal. Depending on the codebase, you might not need the algebraic solutions for all, numeric ones may suffice. If you look at BRL-CAD source (I have not!), you should find the code that handles the primitives. If it is well commented, you should see which equations you need to implement for a new primitive. And, you'll also need to add the support for the heart parameters (orientation, scaling) to the CAD language.
All in all, I suspect this is much more work than you originally envisioned. That is because the primitives in solid constructive geometry are very fundamental objects, and adding a new one -- especially a complicated one -- is not going to be easy.
I hope this helps you. Good luck with your efforts!
P.S. Did you know that ♥ originally depicted two hearts joined together?
O_oAll in all, I suspect this is much more work than you originally envisioned.
Well, sure, when you describe it with all that... maths.
Well, that is a myth, but at least it is a cute one. ^_^Did you know that ♥ originally depicted two hearts joined together?
Soma
It's awesome to see you back, Nominal Animal.
Fact - Beethoven wrote his first symphony in C
Nominal Animal, thanks for your assistance on the heart primitive. It has been very helpful and I see that you are very versed with geometry of solids. Thanks once more. Check
As someone who is not good at math & trying to learn programming, all of this is very scary
As a general rule for anyone interested in graphics, scientific computation or data processing, or just programming for fun, I really recommend learning the basics of vector analysis and combinatorics.
For example, the math for constructive solid geometry for a sphere, right cuboid (like a cube, except sides are rectangles, not squares), or even a cylinder, can be completely done using basic vector math. (There is no need for any differential math, because the surface normal at any point on the surface is already known -- and can be inferred from the definition; no need for differential algebra.)
In particular, concepts like dot product, cross product, norm (length), and the cosine of the angle between two intersecting vectors, are extremely useful. As a real-world example, consider the Trilateration Wikipedia article (the Preliminary computations section) in 2010. It required Euler angles, finding the tangent of an angle, and a boatload of sine and cosine evaluations. Those are very, very slow operations on small microcontrollers. Compare that to the current version, after the simplifications and explanations I recommended and RHB100 implemented. Now, you only need basic (vector) calculus; no sines, cosines, or arcus tangents to compute at all. It is relatively simple an algorithm to implement even on a microcontroller, although you might want to use fixed point types for the vector components instead of real floating point types, for maximum speed.
As to combinatorics, there are many subfields. The subfield you might wish to concentrate on depends on your interests, and the domain of problems you wish to work on. Since it is such a large field, and so very "math-heavy", I recommend approaching it via puzzles. The Traveling salesman problem is probably the most well known combinatorial optimization problem. Looking at good descriptions on different approaches to solving the problem gives you a very good basic working knowledge -- it's surprising how often these solutions crop up in real life, if you write enough code. There is often no need to crack open a math combinatorics book; good explanations and explorations on combinatorial problems are much better, in my opinion.
Personally, I'm not really a math person at all, myself. To me, math itself is not interesting, it is just a tool I use to solve problems. I find the descriptions of logic and algorithms involved in programming much more intuitive and closer to my own thought flow. Even now, when looking up details on some mathematical thing I know the rough outlines for, but need the details to implement it, I still sometimes struggle. It's not the idea or the concepts! It's the language used to describe them. I honestly think most excellent mathematicians are a bit autistic. (That is a real-world observation, definitely not intended in any way as an insult.) There have been a number of physicists that have had the ability to explain very weird physical phenomena to laypeople -- like Richard Feynman, for example --, but I don't know of any really good mathematics popularizers. (If you do, let me know; I'm often wrong, and appreciate being corrected.)
It is sometimes hard to believe, but software is mathematics. The way we read and write software and the various fields of mathematics are just very different; clearly, at least different languages. The approaches are different, the way the concepts and problems and solutions are described are different.. but the fundamental reality, the stuff lurking underneath the concepts that is being described, those are the one and the same.
In practice, this means that if you struggle with something in math or programming, especially when struggling with one to complete a solution in the other, never get frustrated. It is -- and you must realize it is -- just a communications problem. You then need to find another way to discuss the problem, or another way to approach the problem, to proceed. Also, you very soon learn to be as explicit as you can in describing your programming or math problems, and that will alwayshelp you get better, more on-track advice.
In this particular case, the new primitive shape is a very complex one, much more complex than any of the primitives in any solid constructive geometry program I've ever seen. It's not impossible to implement it, I think; it's just very hard.
If I wanted to add a heart primitive to a solid constructive geometry program, I'd definitely do it using ellipsoids instead, with cuboids used to cut the overlapping parts out. I'd first start with two ellipsoids (and two cuboids) arranged symmetically, tilted a bit from the vertical axis, so they form two lobes above, and only one below.
You can use POV-Ray for visual experimentation; it supports basic solid constructive geometry operations.
Since the two ellipsoids are symmetric, they will join at x=0, just not smoothly: the heart would have a small ridge of surface normal discontinuity. Not visible in silhouette, only visible when lighted.
Anyway, I would still use vector math "on paper" -- in other words, Maxima or Maple to find algebraic solutions, or quick C or Python programs to find numerical solutions -- to find out how to position the ellipsoids so that their intersection is where I want; and that the tangents are similar enough at the intersections for the "seam" to be unobtrusive enough.
(I think it might be possible to arrange three intersecting ellipsoids in an Y shape, their intersections forming an upside-down Y, with C1 continuity: not only the surface being continuous, but the change in the surface (curvature) being continuous. If that is correct, it would yield a perfect, seamless heart, suitable for use as a SCG primitive.)
Sorry for the long, semi-off-topic post.
I can't understand people who don't like math equations.
I've been seeing a particularly attractive quadratic equation I met on line for a couple
of weeks now. No nipples, of course, but very curvy. Complex roots too, so interesting
and a tad exotic.Yup, over the last fortnight we must have gone through all the bases.
The only downside is that I appear to be the one handling all the problems in the
relationship. Only yesterday I said "Isn't it time you solved your own problems,huh?",
and what did I get back?
Nada, nothing!
Nice! But.. we're talking about high-degree multivariates, here. Six degrees and at least three independent variables; that's what causes the problems. They just don't get along.Originally Posted by gemera
Just consider for a moment what the partial derivatives look like. And, being quintics, they're not that much easier than the sextics either. Their cross products (if parametrized for u and v) are likely to be tenth degree! (Isn't that further than Kevin Bacon?)
There's nothing wrong in a relatively straightforward equation; it's the nastiness when things get more complicated. Although, I must admit I like my curves and ups and downs, so cubics are my favourite. Wait, we're still talking about maths, right?
Thanks very much Nominal Animal for your depth in this domain of research. I am actually the one who asked the question on Implementing a heart primitive through my friend Check Nyah's profile. we are about to code during the summer for the BRLCAD and have submitted project proposals for this.Are you a researcher in this domain or is it your application domain for programming. Please, I have submitted a proposal to BRLCAD to participate in a project. 1.First off , I would appreciate if you read and scrutinize my proposal on this link User:Izak - BRL-CAD. on Implementing a heart primitive and propose ideas on Sean's review comments on this link Implementation of a heart primitive . 2.Secondly, I have been asked to submit a convincing code patch for this work and I really don't know what to tackle. Could you give me some ideas on how to go about this. 3.Lastly , Can you propose to me some parametric equations for x,z,y to solve the surface equation.
Thanks
I've worked on 3D visualization for atomic/molecular models, including raycasting and raytracing techniques, but I've only toyed with solid constructive geometry.
I'll try to find time, but I'm traveling until the day after tomorrow, so it'll take a few days at least.
Sure, but only after your first and third points are solved first.
I was afraid you'd ask..
I did try my usual tricks to find expressions
where the limits for u and v depend on the mapping used (rectangular/toroidal, triangular, or circular disc).Code:x(u, v) = f1(u, v) y(u, v) = f2(u, v) z(u, v) = f3(u, v)
No luck. It seems that if it is at all possible, the expressions are piecewise and high-degree (not slow to compute, but very hard to find).
If I were you, I'd look for an easier equation for the heart. (I'm assuming a pseudoprimitive -- a primitive created using other primitives but otherwise usable as a CSG primitive; a "macro" -- is out of the question. Otherwise it'd be my approach.)
Usually constructive solid geometry does not need the general equation for the solid (comparable to the equation you started this discussion with), so you don't need to worry about that.
How about creating a parametric surface first? Say, start from a volume of revolution defined by Euler angles φ (rotation around Z axis) and θ (angle to the Z axis), and distance r:
At φ=0° and φ=180° the cross section is an unit circle.Code:0° ≤ φ < 360° 0° ≤ θ < 180° r(φ, θ) = 1 + 6 sin(φ)2 (1 - θ/π) (θ/π)2
At φ=90° and φ=270° the cross section is 1+6(1-θ/π)(θ/π)2, which describes a lobe towards θ=135° or so. The surface equations can be rewritten as
and can be rewritten using texture/patch coordinatesCode:x(φ, θ) = cos(φ) sin(θ) r(φ, θ) = cos(φ) sin(θ) + 6 cos(φ) sin(θ) sin(φ)2 (1 - θ/π) (θ/π)2 y(φ, θ) = sin(φ) sin(θ) r(φ, θ) = sin(φ) sin(θ) + 6 sin(φ) sin(θ) sin(φ)2 (1 - θ/π) (θ/π)2 z(φ, θ) = cos(θ) r(φ, θ) = cos(θ) + 6 cos(θ) sin(φ)2 (1 - θ/π) (θ/π)2
which can be simplified, and even approximated using series expansions of sine and cosine functions. These are obviously differentiable; for example, the X component of the u and v tangents areCode:0 ≤ u,v < 1 x(u, v) = cos(2πu) sin(πv) + 6 cos(2πu) sin(πv) sin(2πu)2 (1 - v) v2 y(u, v) = sin(2πu) sin(πv) + 6 sin(2πu) sin(πv) sin(2πu)2 (1 - v) v2 z(u, v) = cos(πv) + 6 cos(πv) sin(2πu)2 (1 - v) v2 Note: r(u,v) = 1 + 6 sin(2πu)2 (1 - v) v2
If you calculate all six, the cross products of the ∂u and ∂v vectors give the surface normal vectors, which are needed for lighting and reflections.Code:∂ux(u,v) = ∂x(u,v)/∂u = -2π sin(2πu) sin(πv) (1 + 6 v2 - 18 v2 cos(2πu)2 - 6 v3 + 18 v3 cos(2πu)2) ∂vx(u,v) = ∂x(u,v)/∂v = π cos(2πu) cos(πv) + 6π cos(2πu) cos(πv) sin(2πu)2 (1-v) v2 - 6 cos(2π*u) sin(πv) sin(2πu)2 v2 + 12 cos(2πu) sin(πv) sin(2πu)2 (1-v) v
That's just the starting point (although the rest seems very doable). The other math stuff you need to solve include
- Is point P inside the heart?
(Heart defined using center point and X, Y, and Z unit/scaling vectors)- At which point(s) vector V intersects the heart?
- What is the volume of intersection between the heart and (the other primitives in BRL-CAD)?
- What is the surface of intersection between the heart and (the other primitives in BRL-CAD)?
The best way to solve the above depend heavily on how the other primitives in BRL-CAD are computed. (I'm sure you can imagine the importance of having a consistent logical approach in the software overall, without abrupt changes in approach for something as fundamental as the primitive solids!)
I'll take a look at BRL-CAD in a few days, and return to this topic then. No promises, though.
The new primitive is first added to the librt library in BRL-CAD. Each primitive is described by a structure that defines various callback functions. The structures are listed in src/librt/primitives/table.c, with the constants defined in include/rtgeom.h and include/raytrace.h. The src/librt/primitives/ directory contains the sources for the other primitives, so you can look at what needs to be implemented.
After the library implements the new primitive, support for the new primitive must be added to the other parts of BRL-CAD. You can grep for ID_SPH (sphere primitive) or ID_TOR (toroid primitive) et cetera in the source tree to see where the new primitive needs to be added for it to become useable in BRL-CAD.
The librt already provides a metaball primitive (using 1/(x2+y2+z2+), I believe). You could use it in BRL-CAD to create a heart pseudo-primitive, by parametrizing a heart as a metaball, mapping the heart properties (size? lobe size? aspect ratio?) as metaball properties. Of course, you'd first need to experiment with metaballs, and find out a way to parametrize different types of hearts as metaballs.
However you decide to approach adding the heart primitive, I suggest you look at the metaball implementation in BRL-CAD (ID_METABALL), and use that as a guide on how to add a new primitive. Compare implementations between other primitives, to see what callback functions are pertinent for the heart primitive, and what different approaches there are available.
Good luck!
