How to fix this heart?

Question

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

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

I correctly ported this equation to Java, dealing with a few cases outside the box:

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 can't just return -1 , like me:

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

This is not the behavior of the heart at this point. Be that as it may, 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 this drop above! I tried to change this -1 problem to .5 , resulting in the following:

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

How can I fix this problem? I do not need a hole in my heart above, and I do not want a horned heart. If I could fix the horns in the shape of a heart and color everything else correctly, that would be great. Ideally, the two sides of the heart would come together as a dot (the hearts have a small dot in the joint), but if they are bent together, as shown on the horns, that would be good too. How can I achieve this?

Answer 1

The problem is simple. If we look at this area of ​​the horseshoe, we get imaginary numbers. For part of it, this must belong to our heart. In this area, if we were to evaluate our function (in mathematics, but not in programming), the imaginary part of the function was canceled. Therefore, it should look like this (created in Mathematica):

In principle, the function for this part is almost identical; we just have to do arithmetic with complex numbers instead of real numbers. Here's a function that does just 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 program, simply modify it 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 a + bi . So the second argument is just sqrt(-temp1) (note that temp1 < 0 , so I use - instead of Math.abs ; Math.abs is probably the best idea). cbrt_complex returns the cube root of a complex number, in polar form: re iθ . We can see from wolframalpha that with positive r and θ we can write the nth root of complex numbers as follows

And this is exactly how the code for cbrt_complex and sqrt_complex . Note that both have a complex number in rectangular coordinates ( a + bi ) and return a complex number in polar coordinates ( re iθ )

 double[] part1 = polar_reciprocal(temp3); 

It is easier to take the inverse polar complex number than a rectangular complex number. If we have re iθ , then its inverse (this, fortunately, corresponds to the standard power rules) is simply 1/re -iθ . That is why we remain in polar form; the polar form simplifies operations like multiplication, and operations like input are more complicated, and the rectangular form is vice versa.

Note that if we have a polar complex number re iθ and we want to multiply by a real number d , the answer will be as simple as dre iθ .

The toRect function does exactly what it seems: it converts complex complex coordinate numbers into rectangular coordinate complex numbers.

You may have noticed that the if statement does not check for an imaginary part, but only if the imaginary part is really small. This is due to the fact that we use floating point numbers, so checking result[1] == 0 will most likely fail.

And here you are! Please note that we could implement the whole function of the heart with this complex arithmetic of numbers, but it is probably faster to avoid it.