java - 我该如何修复这颗心？

Question

眼看情人节快到了，我决定做一颗心。所以我找到了this heart from mathematica.se :

我在 Mathematica 中尝试(求解 z，切换一些变量)得到心脏 z 值的方程，给定 x 和 y 值(点击查看全尺寸):

我忠实地将这个等式移植到 Java，处理了几个越界的情况:

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)));
    }
}

唯一的问题是当temp1 < 0 ，我不能简单地返回 -1 ，就像我一样:

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

那不是当时心脏的行为。事实上，当我尝试制作我的形象时:

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
}

我明白了:

看看顶部的凹陷!我试着改变那个有问题的 -1到 .5 ，结果是:

现在心脏有角了。但是很明显哪里有问题if的条件得到满足。

我该如何解决这个问题？不想心上有个洞，也不想心有角。如果我可以将牛角剪成心形，并适本地为其余部分涂上颜色，那就太好了。理想情况下，心的两侧会作为一个点连接在一起(心在连接处有一个小点)，但如果它们像角中所示那样弯曲在一起，那也很好。我怎样才能做到这一点？

Answer 1

问题很简单。如果我们观察那个马蹄形区域，我们会得到虚数。对于它的一部分，它应该属于我们的心。在那个区域，如果我们要评估我们的函数(通过数学，而不是编程)，函数的虚部会取消。所以它应该看起来像这样(在 Mathematica 中生成):

基本上，该部分的功能几乎相同；我们只需要用复数而不是实数进行算术运算。下面是一个函数，正是这样做的:

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;
    }
}

要将其加入您的程序，只需更改您的函数以反射(reflect)这一点:

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

运行它，我们得到了想要的结果:

---

很明显，这个新函数实现了完全相同的公式。但是每个部分是如何工作的？

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

cbrt_complex采用 a + b i 形式的复数.这就是为什么第二个参数只是 sqrt(-temp1) (注意 temp1 < 0 ，所以我使用 - 而不是 Math.abs ； Math.abs 可能是一个更好的主意)。 cbrt_complex以极坐标形式返回复数的立方根:r e<sup>iθ</sup> . We can see from wolframalpha积极的 r和 θ ，我们可以写出一个复数的 n 次方根:

这正是 cbrt_complex 的代码和 sqrt_complex工作。请注意，两者都采用直角坐标 (a + b i) 中的复数并返回极坐标中的复数 (r e<sup>iθ</sup>)

double[] part1 = polar_reciprocal(temp3);

取极坐标复数的倒数比取直角复数的倒数容易。如果我们有 r e<sup>iθ</sup> ，它的倒数(幸运的是，这遵循标准的幂规则)就是 1/r e<sup>-iθ</sup> .这实际上就是我们保持极地形式的原因；极坐标形式使乘法运算更容易，而加法运算更难，而矩形形式则相反。

请注意，如果我们有一个极坐标复数 r e<sup>iθ</sup>我们想乘以一个实数 d ，答案很简单 d r e<sup>iθ</sup> .

toRect函数的作用和它看起来的一样:它将极坐标复数转换为直角坐标复数。

您可能已经注意到 if 语句不检查是否有没有虚部，但仅在虚部非常小的情况下才检查。这是因为我们使用的是 float ，所以检查 result[1] == 0很可能会失败。

你来了!请注意，我们实际上可以使用这种复数运算来实现整个心脏功能，但避免这种情况可能会更快。