How do I fix this heart?

Question

Seeing as Valentine's Day is fast approaching, I decided to create a heart. So I found this heart from mathematica.se:

I played around in Mathematica (solved for z, switching some variables around) to get this equation for the z-value of the heart, given the x and y values (click for full-size):

I faithfully ported this equation to Java, dealing with a couple out-of-bounds cases:

import static java.lang.Math.cbrt;
import static java.lang.Math.pow;
import static java.lang.Math.sqrt;

...

public static double heart(double xi, double yi) {
    double x = xi;
    double y = -yi;
    double temp = 5739562800L * pow(y, 3) + 109051693200L * pow(x, 2) * pow(y, 3)
                  - 5739562800L * pow(y, 5);
    double temp1 = -244019119519584000L * pow(y, 9) + pow(temp, 2);
    //
    if (temp1 < 0) {
        return -1; // this is one possible out of bounds location
                   // this spot is the location of the problem
    }
    //
    double temp2 = sqrt(temp1);
    double temp3 = cbrt(temp + temp2);
    if (temp3 != 0) {
        double part1 = (36 * cbrt(2) * pow(y, 3)) / temp3;
        double part2 = 1 / (10935 * cbrt(2)) * temp3;
        double looseparts = 4.0 / 9 - 4.0 / 9 * pow(x, 2) - 4.0 / 9 * pow(y, 2);
        double sqrt_body = looseparts + part1 + part2;
        if (sqrt_body >= 0) {
            return sqrt(sqrt_body);
        } else {
            return -1; // this works; returns -1 if we are outside the heart
        }
    } else {
        // through trial and error, I discovered that this should
        // be an ellipse (or that it is close enough)
        return Math.sqrt(Math.pow(2.0 / 3, 2) * (1 - Math.pow(x, 2)));
    }
}

The only problem is that when temp1 < 0, I cannot simply return -1, like I do:

if (temp1 < 0) {
    return -1; // this is one possible out of bounds location
               // this spot is the location of the problem
}

That's not the behavior of the heart at that point. As it is, when I try to make my image:

import java.awt.image.BufferedImage;
import java.io.File;
import java.io.IOException;
import javax.imageio.ImageIO;

import static java.lang.Math.cbrt;
import static java.lang.Math.pow;
import static java.lang.Math.sqrt;

public class Heart {
    public static double scale(int x, int range, double l, double r) {
        double width = r - l;
        return (double) x / (range - 1) * width + l;
    }

    public static void main(String[] args) throws IOException {
        BufferedImage img = new BufferedImage(1000, 1000, BufferedImage.TYPE_INT_RGB);
        // this is actually larger than the max heart value
        final double max_heart = 0.679;
        double max = 0.0;
        for (int x = 0; x < img.getWidth(); x++) {
            for (int y = 0; y < img.getHeight(); y++) {
                double xv = scale(x, img.getWidth(), -1.2, 1.2);
                double yv = scale(y, img.getHeight(), -1.3, 1);
                double heart = heart(xv, yv); //this isn't an accident
                // yes I don't check for the return of -1, but still
                // the -1 values return a nice shade of pink: 0xFFADAD
                // None of the other values should be negative, as I did
                // step through from -1000 to 1000 in python, and there 
                // were no negatives that were not -1
                int r = 0xFF;
                int gb = (int) (0xFF * (max_heart - heart));
                int rgb = (r << 16) | (gb << 8) | gb;
                img.setRGB(x, y, rgb);
            }
        }
        ImageIO.write(img, "png", new File("location"));
    }
    // heart function clipped; it belongs here
}

I get this:

Look at that dip at the top! I tried changing that problematic -1 to a .5, resulting in this:

Now the heart has horns. But it becomes clear where that problematic if's condition is met.

How can I fix this problem? I don't want a hole in my heart at the top, and I don't want a horned heart. If I could clip the horns to the shape of a heart, and color the rest appropriately, that would be perfectly fine. Ideally, the two sides of the heart would come together as a point (hearts have a little point at the join), but if they curve together like shown in the horns, that would be fine too. How can I achieve this?

Answer 1

The problem is simple. If we look at that horseshoe region, we get imaginary numbers. For part of it, it should belong to our heart. In that region, if we were to evaluate our function (by math, not by programming), the imaginary parts of the function cancel. So it should look like this (generated in Mathematica):

Basically, the function for that part is almost identical; we just have to do arithmetic with complex numbers instead of real numbers. Here's a function that does exactly that:

private static double topOfHeart(double x, double y, double temp, double temp1) {
    //complex arithmetic; each double[] is a single number
    double[] temp3 = cbrt_complex(temp, sqrt(-temp1));
    double[] part1 = polar_reciprocal(temp3);
    part1[0] *= 36 * cbrt(2) * pow(y, 3);
    double[] part2 = temp3;
    part2[0] /= (10935 * cbrt(2));
    toRect(part1, part2);
    double looseparts = 4.0 / 9 - 4.0 / 9 * pow(x, 2) - 4.0 / 9 * pow(y, 2);
    double real_part = looseparts + part1[0] + part2[0];
    double imag_part = part1[1] + part2[1];
    double[] result = sqrt_complex(real_part, imag_part);
    toRect(result);

    // theoretically, result[1] == 0 should work, but floating point says otherwise
    if (Math.abs(result[1]) < 1e-5) {
        return result[0];
    }
    return -1;
}

/**
 * returns a specific cuberoot of this complex number, in polar form
 */
public static double[] cbrt_complex(double a, double b) {
    double r = Math.hypot(a, b);
    double theta = Math.atan2(b, a);
    double cbrt_r = cbrt(r);
    double cbrt_theta = 1.0 / 3 * (2 * PI * Math.floor((PI - theta) / (2 * PI)) + theta);
    return new double[]{cbrt_r, cbrt_theta};
}

/**
 * returns a specific squareroot of this complex number, in polar form
 */
public static double[] sqrt_complex(double a, double b) {
    double r = Math.hypot(a, b);
    double theta = Math.atan2(b, a);
    double sqrt_r = Math.sqrt(r);
    double sqrt_theta = 1.0 / 2 * (2 * PI * Math.floor((PI - theta) / (2 * PI)) + theta);
    return new double[]{sqrt_r, sqrt_theta};
}

public static double[] polar_reciprocal(double[] polar) {
    return new double[]{1 / polar[0], -polar[1]};
}

public static void toRect(double[]... polars) {
    for (double[] polar: polars) {
        double a = Math.cos(polar[1]) * polar[0];
        double b = Math.sin(polar[1]) * polar[0];
        polar[0] = a;
        polar[1] = b;
    }
}

To join this with your program, simply change your function to reflect this:

if (temp1 < 0) {
    return topOfHeart(x, y, temp, temp1);
}

And running it, we get the desired result:

---

It should be pretty clear that this new function implements exactly the same formula. But how does each part work?

double[] temp3 = cbrt_complex(temp, sqrt(-temp1));

cbrt_complex takes a complex number in the form of a + b i. That's why the second argument is simply sqrt(-temp1) (notice that temp1 < 0, so I use - instead of Math.abs; Math.abs is probably a better idea). cbrt_complex returns the cube root of the complex number, in polar form: r eiθ. We can see from wolframalpha that with positive r and θ, we can write an n-th root of a complex numbers as follows:

And that's exactly how the code for the cbrt_complex and sqrt_complex work. Note that both take a complex number in rectangular coordinates (a + b i) and return a complex number in polar coordinates (r eiθ)

double[] part1 = polar_reciprocal(temp3);

It is easier to take the reciprocal of a polar complex number than a rectangular complex number. If we have r eiθ, its reciprocal (this follows standard power rules, luckily) is simply 1/r e-iθ. This is actually why we are staying in polar form; polar form makes multiplication-type operations easier, and addition type operations harder, while rectangular form does the opposite.

Notice that if we have a polar complex number r eiθ and we want to multiply by a real number d, the answer is as simple as d r eiθ.

The toRect function does exactly what it seems like it does: it converts polar coordinate complex numbers to rectangular coordinate complex numbers.

You may have noticed that the if statement doesn't check that there is no imaginary part, but only if the imaginary part is really small. This is because we are using floating point numbers, so checking result[1] == 0 will likely fail.

And there you are! Notice that we could actually implement the entire heart function with this complex number arithmetic, but it's probably faster to avoid this.