Search code examples
c++copenclfractalsnewtons-method

Newton Fractal generation

I wanted to write my own newton fractal generator.. It's using OpenCL... but that's not the problem.. my problem is that atm. only veerryy few pixels are converging.

So to explain what I've done so far:

  • I've selected a function I wanted to use: f(z)=z^4-1 (for testing purposes)
  • I've calculated the roots of this function: 1+0î, -1+0î, 0+1î, 0-1î
  • I've written a OpenCL Host and Kernel:
    • the kernel uses a struct with 4 doubles: re (real), im (imaginary), r (as abs), phi (as argument, polar angle or how you call it)
    • computes from resolution, zoom and global_work_id the "type" of the pixel and the intensity - where type is the root the newton method is converging to / whether it's diverging

Here's what I get rendered: Fractal1

Here is the whole kernel:

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

And here is the coloration of the image:

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

I had some mistakes in the complex number computations but I check everything today again and again.. I think they should be okay now.. but what else could be the reason that there are just this few start values converging? Did I do something wrong with newton's method?

Thanks for all your help!!

Solution

  • Well somewhat it really helped to run the code as normal C code.. as this makes Debugging easier: so the issue were some coding issues which I have been able to solve now.. for example my pow function was corrupted and when I added or subtracted I forgot to set the imaginary part to the temp complex number .. so here's my final OpenCL kernel:

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

    hope it might help someone someday.. :)