How can I generate an object from a mathematical function?

Question

I have a formula that mathematically describes the heart surface in Cartesian coordinates:

$$\left(x^2+\frac{9y^2}{4}+z^2-1\right)^3-x^2z^3-\frac{9y^2z^3}{80}=0$$

Unfortunately, I have no idea how to implement it in Blender.

Explicitly, I want to use for the Mesh Math Function. I already activated the Addon ADD Mesh: Extra Objects.

I searched for a tutorial, also for a syntax, I have found nothing, something near it.

I created (guess) from other functions this syntax.

(X**2+ ((9*y**2)/4)+z**2-1)**3-((x**2)*(z**3)-((9*y**2*z**3)/80)=0

I found another formula:

5 + (-sqrt(1 - x**2 - (y-abs(x))**2)) * cos(30*((1 - x**2 - (y-abs(x))**2)))

x is from -1 to 1, y is from -1 to 1.5, z is from 1 to 8

Answer 1

Have used the equations found here http://www.econym.demon.co.uk/isotut/real.htm#heart1

The first can be crunched into XYZ surface. Notice the use of fabs(...) in the A helper function. The abs(...) method is not in the list of allowable methods and creates an error.

Here is some code for the second. Which required a mapping from spherical coordinates. Run the code to produce a heart mesh. Change the line h = heart(1.0, 1.0, 0.2) to produce other hearts.

import bpy
import bmesh
from math import degrees, radians, sin, cos, tan
from mathutils import Vector
from bpy import context

class SphericalCoordsPoint:
    @property
    def xyz(self):
        theta = self.theta
        zi = self.zi
        x = cos(theta) * sin(zi)
        y = sin(theta) * sin(zi)
        z = cos(zi)
        R = self.R
        return R * Vector((x,y,z))

    def __init__(self, R, theta, zi):
        self.R = R
        self.theta = theta
        self.zi = zi
        #self.xyz = self.point(theta, zi)

    def __repr__(self):
        return "SphericalCoord(%.4f, %.4f)" % (degrees(self.theta),
                                               degrees(self.zi))
# define the heart method.
def heart(a, b, c):
    def heart(v):
        x = a * v.x
        y = b * (v.y -pow(abs(v.x),0.8))
        z = c * (v.z)
        return Vector((x, y, z))
    return heart

# make the heart bmesh
bm = bmesh.new()

# TODO come up with a nicer way to do this.
sphere_rings = [[SphericalCoordsPoint(1, radians(theta), radians(zi)) 
                 for theta in range (0, 360, 10)]
                 for zi in range(0, 181, 10)]

h = heart(1.0, 1.0, 0.2)

verts0 = [bm.verts.new(h(p.xyz)) for p in sphere_rings[0]]
verts0.append(verts0[0])
for ring in range(1, len(sphere_rings)):

    verts1 = [bm.verts.new(h(p.xyz)) for p in sphere_rings[ring]]
    verts1.append(verts1[0])

    faces = [
        bm.faces.new((
            verts0[i], verts1[i],
            verts1[i+1], verts0[i+1]
        ))
        for i in range(len(verts0) - 1)
    ]
    verts0 = verts1

# create mesh link it to scene etc
mesh = bpy.data.meshes.new("heart")
bm.to_mesh(mesh)
obj = bpy.data.objects.new("heart", mesh)
scene = context.scene

context.collection.objects.link(obj)
context.view_layer.objects.active = obj
obj.select_set(True)
obj.location = scene.cursor.location

which produced this result filled like the UV sphere (the node is in the middle of the heart

Can change the shape by tweaking the scalars in the heart equation. For example

def heart(v):
    x = v.x
    y = v.y -pow(abs(v.x),0.8)*1
    z = v.z * 0.5
    return Vector((x, y, z))

EDIT: this can now be done with h = heart(1.0, 1.0, 0.5)

produces

Samples of both using Spherical and Toroidal spaces

Another similarish object is the cardioid

Answer 2

Computing a real mathematical set given implicitly by function(s) f1(x,y,z,...)=0, f2(x,y,z,...)=0, ... is really hard. Here's how I do what you want:

• The implicit function theorem tell us that, off a set of measure 0, an analytic set can be locally parameterized by a number of parameters equal to the dimension of the set.

• Your heart is a two-dimensional algebraic variety. Thus, we can compute a set of local patches, the union of which when taken together with the problematic set of measure 0, will describe the heart.

I wrote a computer program which does just this over the last several years as a postdoc. It's called Bertini_real, and is part of the Bertini suite of software products for computing solutions to algebraic equations. The field of mathematics this software comes from is called Numerical Algebraic Geometry.

First, we need to get some complex points on the surface. For this, we write an input file for Bertini, and call bertini, using tracktype 1 to compute the witness set. In this case, we expect six points, because the component is a hypersurface, which always has degree equal to the degree of the polynomial.

%file `input` for Bertini
%the real part describes a mathematical heart.
CONFIG
tracktype: 1; %compute the positive dimensional components
sharpendigits: 60; %sharpen to ludicrous accuracy
condnumthreshold: 1e300; %it's hard to be singular unless multiplicity>1
securitylevel: 1; %do not stop tracking paths going to infinity
odepredictor: 2; %use RK45
endgamenum: 2; %use the Cauchy endgame for computing singular points
numsamplepoints: 5; %use a total of 6 points in the endgame
END;

INPUT
variable_group x, y, z; %these are our variables
function f;  %we'll have one function.
f = (x^2+9/4*y^2+z^2-1)^3 - x^2*z^3-9/80*y^2*z^3; %here's the function.
% we assume all functions equal 0
END;

Great, now Bertini has confirmed that the surface is indeed of degree 6. Nice. Next is to get just the real part of this complex object. This is where Bertini_real comes in. The software unfortunately depends on Matlab right now, but we are working to replace it with Python so the software is completely free. Almost there as of December 2016.

Next, invoke bertini_real on the same input file. A whole bunch of stuff happens, including computing that set of measure 0 which cannot be well-parameterized by the two parameters. By the way, the two parameters we use for the computation of your surface are random orthogonal linear projections of the coordinates of the surface. The set of measure 0 is called the critical curve, and in this case consists of a curve which depends on the random projection, and a singular curve about the waist of the heart, as show below:

It doesn't look much like the heart you are after, and that's because this is the raw decomposition, and it needs to be refined. Invoking the sampler from Bertini_real will refine it, using your settings. Here's an older refinement I produced a few months ago:

Looks better, right? I used my plotting software in Matlab for these images. Now, we export the surface to .stl format. Perhaps there are better formats. I use .stl. Importing into Blender reveals some faces are mis-oriented. Bertini_real does not compute normals, because the normal is ill-defined in higher ambient dimensions than 3. So, I let Blender do it for me.

Using the recalculate normals function in Edit mode in Blender fixes this issue. At this point, we now have a well-oriented mathematical heart in Blender.

I take it one step further, and 3D print it. Here's a result:

We can do this process for any algebraic surface. You do not need to fuss around with trying to compute z = f(x,y), because that happens implicitly as part of Bertini_real's decomposition. The set on which this parameterization fails is called the critical curve, and is itself decomposed with the curve decomposition routine.

If you have Matlab, you can do this yourself. If you know how to program in Python, perhaps you'd be interested in contributing to my software, in the form of visualization and .stl export directly from Python, so I can eliminate Matlab?

The more complicated the surface, the more challenging the decomposition. Surfaces such as Barth's Dodecahedric are really hard to get to compute 100% correctly. But nice ones such as this heart, aren't that bad.

Check out Herwig Hauser's gallery of algebraic surfaces, and my reproduction in plastic.

If you want to learn more about the software I use, or the theory behind it, I recommend the book that goes with the software Bertini, called Numerically solving polynomial systems with Bertini. Not an intro book, but it tells you how to do things, and has worked examples.

Also, I made a 37 minute real-time video of me going from equation to printable object, using the above process including Blender for a significant portion of post-processing, available here.

I finally conclude that this answer does use software that is not part of the Blender world, but due to the genuine challenge of rendering implicit mathematical objects, a generic solution to this question is probably not going to be present in Blender any time soon.

Thanks for reading!

edit: since I now have sufficient rep, I added more links, and several concluding paragraphs

edit 2021-09-21: updated links

Answer 3

If you require that specific formula, you're in for a rough ride, as it would require solving a sixth degree equation to isolate z. There are no general methods for algebraicly solving equations of degrees higher than 4. Special cases exist, where you can solve higher degree equations by simplifications, substitutions or factorisations. If you do need this formula, you should probably ask over on Math.SE if anyone can help you to either isolate z or convert the equation into a set of parametric equations for x, y and z.

I use hearts to some extent in my scenes, and I'll show you how I make them.
First, you need the Add Mesh: Extra Objects addon. To enable it, go to user preferences, CtrlAltU or File -> User Preferences..., go to the Add-ons tab and click Add Mesh, then click the check box next to Add Mesh: Extra Objects to enable it. To make the changes stick, click Save User Settings, then close the user preferences.

Now click Add -> Mesh -> Math Function -> XYZ Math Surface.

In the options that appear at the botton of the tool shelf, enter the following (if the tool shelf isn't there, press T to show it.

X Equation
sin(u)**3

Y Equation
0

Z Equation
(13*cos(u)-5*cos(2*u)-2*cos(3*u)-cos(4*u))/16

U min 0

U max
2*pi

U step sets how many vertices you want along the curve.

Enable U wrap to close the curve.

Set both V min and V max to 0 and V step to 1. Disable V wrap and Close V. The remaining options can be left untouched.

Tab into edit mode and make sure everything is selected (if it isn't, press A as many times as needed to select everything). Press CtrlV followed by D to remove double vertices. Press E to extrude then Esc to cancel moving the extrusion. Press S followed by some small number, for example .025, to scale the extruded vertices down. Then press Enter to confirm the scaling, then F to make a face. Usually I'd scale to 0 or collapse the vertices, but in this case I don't recomment it. This would create a bunch of triangles, and I normally tend to prefer them over n-gons, but in this case, the smallest angle of the tris around the cusps of the heart would be too small, and cause more distorsions when subdividing than an n-gon.
Press A twice to select everything, then press CtrlN to make the normals consistent. Press CtrlR then hover the mouse over one of the radial edges and scroll the mouse wheel to make some loop cuts, then press Enter to confirm the cuts followed by Esc to finish without sliding the cuts.

Tab out of edit mode and, if you so require or desire, add Solidify and/or Subsurf modifier(s).

The amount of bevelling around the edges depends on the distance between the outher edge and the closest loop cut, so in essence it can be controled by decreasing or increasing the number of loop cuts made in the steps between the third and fourth image above, or by scaling those loop cuts down or up.