Я хотел написать свой собственный фрактальный генератор Ньютона. Он использует 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 «типа» пикселя и интенсивности - где типа является корнем метод ньютона сходится к/ли это расходящийся
Вот что я получаю оказанный :
Вот целое ядро:
#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);
}
}
}
У меня был некоторые ошибки в сложных числовых вычислениях, но проверить все сегодня снова и снова .. Я думаю, что они должны быть хорошо сейчас .. но что еще может быть причиной того, что есть только эти несколько начальных значений сходящихся? Я сделал что-то не так с методом Ньютона?
Спасибо за помощь!
вы обеспечили, что свойства сходимости выполнены? Метод Ньютона не является гарантией конвергенции, если функция «хорошо себя ведет» рядом с корнем, а начальная начальная точка «близка» к корню. Посмотрите http://en.wikipedia.org/wiki/Newton%27s_method#Proof_of_quadratic_convergence_for_Newton.27s_iterative_method –
Если что-то странное происходит, я всегда конвертирую свой код ядра в собственный метод (например, Java, C++) и выполняю отладку с помощью странные значения. Если это приведет к такому же поведению, вы можете быть абсолютно уверены, что ваш код OpenCL работает нормально. Таким образом, проблема может быть граничными условиями (или пользователем перед ПК). – Christian
эй, спасибо за ваши комментарии до сих пор - я посмотрю статью wiki. и Yep Я конвертировал в C++ Code, и он работает так же, как и код OpenCL - всего за 15 минут вместо 15 секунд :-) – fodinabor