2015-03-08 7 views
1

Я хотел написать свой собственный фрактальный генератор Ньютона. Он использует OpenCL ... но это не проблема. Моя проблема в том, что atm. только veerryy несколько пикселей сходятся.Newton Fractal generation

Так, чтобы объяснить, что я сделал до сих пор:

  • Я выбрал функцию Я хотел использовать: F (г) = г^4-1 (для тестирования)
  • Я вычислил корни этой функции: 1 + 0i, -1 + 0i, 0 + 1i, 0-1î
  • Я написал OpenCL хоста и ядро:
    • ядро ​​использует-структуру с 4 удваивает: re (real), im (мнимый), r (as abs), phi (в качестве аргумента, полярный угол или как вы его называете)
    • вычисляет от разрешения, масштабировании и global_work_id «типа» пикселя и интенсивности - где типа является корнем метод ньютона сходится к/ли это расходящийся

Вот что я получаю оказанный : Fractal1

Вот целое ядро:

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

У меня был некоторые ошибки в сложных числовых вычислениях, но проверить все сегодня снова и снова .. Я думаю, что они должны быть хорошо сейчас .. но что еще может быть причиной того, что есть только эти несколько начальных значений сходящихся? Я сделал что-то не так с методом Ньютона?

Спасибо за помощь!

+0

вы обеспечили, что свойства сходимости выполнены? Метод Ньютона не является гарантией конвергенции, если функция «хорошо себя ведет» рядом с корнем, а начальная начальная точка «близка» к корню. Посмотрите http://en.wikipedia.org/wiki/Newton%27s_method#Proof_of_quadratic_convergence_for_Newton.27s_iterative_method –

+0

Если что-то странное происходит, я всегда конвертирую свой код ядра в собственный метод (например, Java, C++) и выполняю отладку с помощью странные значения. Если это приведет к такому же поведению, вы можете быть абсолютно уверены, что ваш код OpenCL работает нормально. Таким образом, проблема может быть граничными условиями (или пользователем перед ПК). – Christian

+0

эй, спасибо за ваши комментарии до сих пор - я посмотрю статью wiki. и Yep Я конвертировал в C++ Code, и он работает так же, как и код OpenCL - всего за 15 минут вместо 15 секунд :-) – fodinabor

ответ

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

надеюсь, что это может помочь кому-нибудь .. :)