c++ - 牛顿分形生成

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

因此，请解释一下我到目前为止所做的事情:

我选择了一个我想使用的函数:f(z)= z ^ 4-1(出于测试目的)

我已经计算出此函数的根:1 +0î，-1 +0î，0 +1î，0-1î

我已经编写了OpenCL主机和内核:

内核使用具有4个 double 的结构: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);
        }
    }
}

我在复数计算中有一些错误，但是我一次又一次检查所有内容..我认为它们现在应该可以了..但是还有什么可能是因为这几个初始值收敛的原因呢？我用牛顿的方法做错了吗？

感谢你的帮助!!

最佳答案

好吧，这确实有助于将代码作为普通的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;
}

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

关于c++ - 牛顿分形生成，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/28930151/