牛顿分形生成

Question

我想编写自己的牛顿分形生成器。它使用的是 OpenCL...但这不是问题。我的问题是 atm。只有很少的像素会聚。

所以解释一下我到目前为止所做的事情：

我选择了一个我想使用的函数：f(z)=z^4-1（用于测试目的）
我计算了这个函数的根：1+0，-1+0，0+1，0-1
我写了一个 OpenCL 主机和内核：
- 内核使用具有 4 个双精度数的结构：re（实数）、im（虚数）、r（作为 abs）、phi（作为参数、极角或您如何称呼它）
- 根据分辨率、缩放和 global_work_id 计算像素的 "type" 和强度 - 其中类型是牛顿法收敛到的根/是否发散

这是我得到的渲染：

这是整个内核：

#pragma OPENCL EXTENSION cl_khr_fp64 : enable
#define pi 3.14159265359

struct complex {
    double im;
    double re;
    double r;
    double phi;
};

struct complex createComplexFromPolar(double _r, double _phi){
    struct complex t;
    t.r = _r;
    t.phi = _phi;

    t.re = cos(t.phi)*t.r;
    t.im = sin(t.phi)*t.r;

    return t;
}

struct complex createComplexFromKarthes(double real, double imag){
    struct complex t;
    t.re = real;
    t.im = imag;

    t.phi = atan(imag / real);
    t.r = sqrt(t.re*t.re + t.im*t.im);

    return t;
}

struct complex recreateComplexFromKarthes(struct complex t){
    return t = createComplexFromKarthes(t.re, t.im);
}

struct complex recreateComplexFromPolar(struct complex t){
    return t = createComplexFromPolar(t.r, t.phi);
}

struct complex addComplex(const struct complex z, const struct complex c){
    struct complex t;
    t.re = c.re + z.re;
    t.im = c.im + z.im;
    return recreateComplexFromKarthes(t);
}

struct complex subComplex(const struct complex z, const struct complex c){
    struct complex t;
    t.re = z.re - c.re;
    t.im = z.im - c.im;
    return recreateComplexFromKarthes(t);
}

struct complex addComplexScalar(const struct complex z, const double n){
    struct complex t;
    t.re = z.re + n;
    return recreateComplexFromKarthes(t);
}

struct complex subComplexScalar(const struct complex z, const double n){
    struct complex t;
    t.re = z.re - n;
    return recreateComplexFromKarthes(t);
}

struct complex multComplexScalar(const struct complex z, const double n) {
    struct complex t;
    t.re = z.re * n;
    t.im = z.im * n;
    return recreateComplexFromKarthes(t);
}

struct complex multComplex(const struct complex z, const struct complex c) {
    return createComplexFromPolar(z.r*c.r, z.phi + c.phi);
}

struct complex powComplex(const struct complex z, int i){
    struct complex t = z;
    for (int j = 0; j < i; j++){
        t = multComplex(t, z);
    }
    return t;
}

struct complex divComplex(const struct complex z, const struct complex c) {
    return createComplexFromPolar(z.r / c.r, z.phi - c.phi);
}

bool compComplex(const struct complex z, const struct complex c, float comp){
    struct complex t = subComplex(z, c);
    if (fabs(t.re) <= comp && fabs(t.im) <= comp)
        return true;
    return false;
}

__kernel void newtonFraktal(__global const int* res, __global const int* zoom, __global int* offset, __global const double* param, __global int* result, __global int* resType){
    const int x = get_global_id(0) + offset[0];
    const int y = get_global_id(1) + offset[1];

    const int xRes = res[0];
    const int yRes = res[1];

    const double a = (x - (xRes / 2)) == 0 ? 0 : (double)(zoom[0] / (x - (double)(xRes / 2)));
    const double b = (y - (yRes / 2)) == 0 ? 0 : (double)(zoom[1] / (y - (double)(yRes / 2)));

    struct complex z = createComplexFromKarthes(a, b);
    struct complex zo = z;

    struct complex c = createComplexFromKarthes(param[0], param[1]);

    struct complex x1 = createComplexFromKarthes(1,0);
    struct complex x2 = createComplexFromKarthes(-1, 0);
    struct complex x3 = createComplexFromKarthes(0, 1);
    struct complex x4 = createComplexFromKarthes(0, -1);

    resType[x + xRes * y] = 3;

    int i = 0;
    while (i < 30000 && fabs(z.r) < 10000){
        z = subComplex(z, divComplex(subComplexScalar(powComplex(z, 4), 1), multComplexScalar(powComplex(z, 3), 4)));

        i++;
        if (compComplex(z, x1, 0.05)){
            resType[x + xRes * y] = 0;
            break;
        } else if (compComplex(z, x2, 0.05)){
            resType[x + xRes * y] = 1;
            break;
        } else if (compComplex(z, x3, 0.05)){
            resType[x + xRes * y] = 2;
            break;
        }
    }
    if (fabs(z.r) >= 10000){
        resType[x + xRes * y] = 4;
    }
    result[x + xRes * y] = i;
}

这是图像的着色：

const int xRes = core->getXRes();
const int yRes = core->getYRes();
for (int y = 0; y < fraktal->getHeight(); y++){
    for (int x = 0; x < fraktal->getWidth(); x++){
        int conDiv = genCL->result[x + y * xRes];
        int type = genCL->typeRes[x + y * xRes];
        if (type == 0){
            //converging to x1
            fraktal->setPixel(x, y, 1*conDiv, 1*conDiv, 0, 1);
        } else if (type == 1){
            //converging to x2
            fraktal->setPixel(x, y, 0, 0, 1*conDiv, 1);
        } else if (type == 2){
            //converging to x3
            fraktal->setPixel(x, y, 0, 1*conDiv, 0, 1);
        } else if (type == 3){
            //diverging and interrupted by loop end
            fraktal->setPixel(x, y, 1*conDiv, 0, 0, 1);
        } else {
            //diverging and interrupted by z.r > 10000
            fraktal->setPixel(x, y, 1, 1, 1, 0.1*conDiv);
        }
    }
}

我在复数计算中犯了一些错误，但我今天一遍又一遍地检查所有内容。我认为它们现在应该没问题了。但是还有什么可能是这几个起始值收敛的原因？我是不是用牛顿法做错了什么？

感谢您的帮助！！

Answer 1

在某种程度上，它确实有助于运行将代码作为普通 C 代码。因为这使调试更容易：所以问题是一些我现在能够解决的编码问题。例如我的pow 函数已损坏，当我添加或减去时，我忘记将虚部设置为临时复数。所以这是我最终的 OpenCL 内核：

#pragma OPENCL EXTENSION cl_khr_fp64 : enable
#define pi 3.14159265359

struct complex {
    double im;
    double re;
    double r;
    double phi;
};

struct complex createComplexFromPolar(double _r, double _phi){
    struct complex t;
    t.r = _r;
    t.phi = _phi;

    t.re = _r*cos(_phi);
    t.im = _r*sin(_phi);

    return t;
}

struct complex createComplexFromKarthes(double real, double imag){
    struct complex t;
    t.re = real;
    t.im = imag;

    t.phi = atan2(imag, real);
    t.r = sqrt(t.re*t.re + t.im*t.im);

    return t;
}

struct complex recreateComplexFromKarthes(struct complex t){
    return createComplexFromKarthes(t.re, t.im);
}

struct complex recreateComplexFromPolar(struct complex t){
    return createComplexFromPolar(t.r, t.phi);
}

struct complex addComplex(const struct complex z, const struct complex c){
    return createComplexFromKarthes(c.re + z.re, c.im + z.im);
}

struct complex subComplex(const struct complex z, const struct complex c){
    return createComplexFromKarthes(z.re - c.re, z.im - c.im);
}

struct complex addComplexScalar(const struct complex z, const double n){
    return createComplexFromKarthes(z.re + n,z.im);
}

struct complex subComplexScalar(const struct complex z, const double n){
    return createComplexFromKarthes(z.re - n, z.im);
}

struct complex multComplexScalar(const struct complex z, const double n){
    return createComplexFromKarthes(z.re * n,z.im * n);
}

struct complex multComplex(const struct complex z, const struct complex c) {
    return createComplexFromKarthes(z.re*c.re-z.im*c.im, z.re*c.im+z.im*c.re);
    //return createComplexFromPolar(z.r*c.r, z.phi + c.phi);
}

struct complex powComplex(const struct complex z, int i){
    struct complex t = z;
    for (int j = 0; j < i-1; j++){
        t = multComplex(t, z);
    }
    return t;
}

struct complex divComplex(const struct complex z, const struct complex c) {
        return createComplexFromPolar(z.r / c.r, z.phi-c.phi);
}

bool compComplex(const struct complex z, const struct complex c, float comp){
    if (fabs(z.re - c.re) <= comp && fabs(z.im - c.im) <= comp)
        return true;
    return false;
}

__kernel void newtonFraktal(__global const int* res, __global const int* zoom, __global int* offset, __global const double* param, __global int* result, __global int* resType){
    const int x = get_global_id(0) + offset[0];
    const int y = get_global_id(1) + offset[1];

    const int xRes = res[0];
    const int yRes = res[1];

    const double a = (x - (xRes / 2)) == 0 ? 0 : (double)((x - (double)(xRes / 2)) / zoom[0]);
    const double b = (y - (yRes / 2)) == 0 ? 0 : (double)((y - (double)(yRes / 2)) / zoom[1]);

    struct complex z = createComplexFromKarthes(a, b);

    //struct complex c = createComplexFromKarthes(param[0], param[1]);

    struct complex x1 = createComplexFromKarthes(0.7071068, 0.7071068);
    struct complex x2 = createComplexFromKarthes(0.7071068, -0.7071068);
    struct complex x3 = createComplexFromKarthes(-0.7071068, 0.7071068);
    struct complex x4 = createComplexFromKarthes(-0.7071068, -0.7071068);

    struct complex f, d;

    resType[x + xRes * y] = 11;

    int i = 0;
    while (i < 6000 && fabs(z.r) < 10000){
        f = addComplexScalar(powComplex(z, 4), 1);
        d = multComplexScalar(powComplex(z, 3), 3);

        z = subComplex(z, divComplex(f, d));

        i++;
        if (compComplex(z, x1, 0.0000001)){
            resType[x + xRes * y] = 0;
            break;
        } else if (compComplex(z, x2, 0.0000001)){
            resType[x + xRes * y] = 1;
            break;
        } else if (compComplex(z, x3, 0.0000001)){
            resType[x + xRes * y] = 2;
            break;
        } else if (compComplex(z, x4, 0.0000001)){
            resType[x + xRes * y] = 3;
            break;
        }
    }
    if (fabs(z.r) >= 1000){
        resType[x + xRes * y] = 10;
    }
    result[x + xRes * y] = i;
}

希望有一天它能对某人有所帮助.. :)

牛顿分形生成

Newton Fractal generation

c++

c

opencl

fractals

newtons-method