Граница множества Мандельброта

Я смотрю, это будет первая тема в данном разделе.
Понадобилось мне тут найти приближённо границу множества Мандельброта для [math]z^4+c[/math]. Создал шейдер:

Код (C):

1. float2 Mult(float2 z1,float2 z2)
2. {
3.     return (float2)(z1.x*z2.x-z1.y*z2.y, z1.y*z2.x+z1.x*z2.y);
4. }
5.  
6. float2 DoIter(float2 z,float2 c)
7. {
8.     return Mult(Mult(Mult(Mult(z,z),z),z),z)+c;
9. }
10.  
11. __kernel void CalcMandel(float Border,float Threshold,uint NumOfIter,uint Size,__global uchar *Data/*,__global uint *NumOfPoints*/)
12. {
13.     size_t _X=get_global_id(0);
14.     size_t _Y=get_global_id(1);
15.     if((_X<Size) && (_Y<Size))
16.     {
17.         float2 c=(float2)(Border*2.0f*((float)_X/(float)(Size-1)-0.5f),Border*2.0f*((float)(Size-_Y-1)/(float)(Size-1)-0.5f));
18.         float2 z=(float2)(0,0);
19.         for(uint i=0;(i<NumOfIter) && (length(z)<Threshold);i++)
20.         {
21.             z=DoIter(z,c);
22.         }
23.         uchar Val;
24.         if(length(z)<Threshold)
25.         {
26.             Val=1;
27.         }
28.         else
29.         {
30.             Val=0;
31.         }
32.         Data[_X+Size*_Y]=Val;
33.         barrier(CLK_GLOBAL_MEM_FENCE);
34.         if(Val==1)
35.         {
36.             if((Data[_X+Size*_Y-1]==1) && (Data[_X+Size*_Y+1]==1))
37.             {
38.                 if((Data[_X+Size*(_Y-1)-1]==1) && (Data[_X+Size*(_Y-1)]==1) && (Data[_X+Size*(_Y-1)+1]==1))
39.                 {
40.                     if((Data[_X+Size*(_Y+1)-1]==1) && (Data[_X+Size*(_Y+1)]==1) && (Data[_X+Size*(_Y+1)+1]==1))
41.                     {
42.                         Val=0;
43.                     }
44.                 }
45.             }
46.         }
47.         barrier(CLK_GLOBAL_MEM_FENCE);
48.         Data[_X+Size*_Y]=Val;
49.     }
50.     else
51.     {
52.         barrier(CLK_GLOBAL_MEM_FENCE);
53.         barrier(CLK_GLOBAL_MEM_FENCE);
54.     }
55. }

Передаю в kernel параметры:
Border=1.5
Threshold=3.0
NumOfIter=500
Size=1000
Data - массив байт.
Жду на выходе границу, а получаю:

на этой картинке чёрные точки соответствуют значению 1, а белые - 0.
Что я делаю не так?

 

Добавлю все исходники.

 

Просьба также выкладывать EXE. Спасибо.

 

Вот EXE, сохраняет картинку в C:\temp\Mandel.bmp

 

_qwe8013
Хоть бы написал, что bmp сохраняется в C:\temp. У меня в системе например нет этой папки, пришлось создать. А вообще для сохранения в temp нужно использовать переменную среды %TEMP%.

 

Не совсем то но прикольно

Код (C):

1. float2 Mult(float2 z1,float2 z2)
2. {
3.    return (float2)(z1.x*z2.x-z1.y*z2.y, z1.y*z2.x+z1.x*z2.y);
4. }
5.  
6. float2 DoIter(float2 z,float2 c)
7. {
8.    return Mult(Mult(Mult(Mult(z,z),z),z),z)+c;
9. }
10.  
11. __kernel void CalcMandel(float Border,float Threshold,uint NumOfIter,uint Size,__global uchar *Data/*,__global uint *NumOfPoints*/)
12. {
13.    size_t _X=get_global_id(0);
14.    size_t _Y=get_global_id(1);
15.    if((_X<Size) && (_Y<Size))
16.    {
17.      float2 c=(float2)(Border*2.0f*((float)_X/(float)(Size-1)-0.5f),Border*2.0f*((float)(Size-_Y-1)/(float)(Size-1)-0.5f));
18.      float2 z=(float2)(0,0);
19.      for(uint i=0;(i<NumOfIter) && (length(z)<Threshold);i++)
20.      {
21.        z=DoIter(z,c);
22.      }
23.      uchar Val;
24.      if(round(length(z)*0.2)==round(Threshold*0.2))
25.      {
26.        Val=1;
27.      }
28.      else
29.      {
30.        Val=0;
31.      }
32.      Data[_X+Size*_Y]=Val;
33.      barrier(CLK_GLOBAL_MEM_FENCE);
34.      barrier(CLK_GLOBAL_MEM_FENCE);
35.    }
36.    else
37.    {
38.      barrier(CLK_GLOBAL_MEM_FENCE);
39.      barrier(CLK_GLOBAL_MEM_FENCE);
40.    }
41. }

 

Мне кажется, что это условие нужно запускать после того, как массив Data будет заполнен и результат выводить в другой массив (например Data2)

Код (C):

1.  
2. barrier(CLK_GLOBAL_MEM_FENCE);
3.         if(Val==1)
4.         {
5.             if((Data[_X+Size*_Y-1]==1) && (Data[_X+Size*_Y+1]==1))
6.             {
7.                 if((Data[_X+Size*(_Y-1)-1]==1) && (Data[_X+Size*(_Y-1)]==1) && (Data[_X+Size*(_Y-1)+1]==1))
8.                 {
9.                     if((Data[_X+Size*(_Y+1)-1]==1) && (Data[_X+Size*(_Y+1)]==1) && (Data[_X+Size*(_Y+1)+1]==1))
10.                     {
11.                         Val=0;
12.                     }
13.                 }
14.             }
15.         }
16.         barrier(CLK_GLOBAL_MEM_FENCE);
17.         Data2[_X+Size*_Y]=Val;

Не?

 

В #4 написано.

Зачем Data2, если Data весь заполнен (для этого я и поставил barrier(CLK_GLOBAL_MEM_FENCE)).

 

Извиняюсь - прочитал немного про OpenCL теперь стало понятнее.

Если втупую расчитывать каждую соседнюю точку - результат нормальный

Код (C):

1. float2 Mult(float2 z1,float2 z2)
2. {
3.    return (float2)(z1.x*z2.x-z1.y*z2.y, z1.y*z2.x+z1.x*z2.y);
4. }
5.  
6. constant float Threshold;
7. constant uint  Size;
8. constant float Border;
9. constant uint  NumOfIter;
10.  
11. uchar Mandel(size_t _X, size_t _Y){
12.    float2 c=(float2)(Border*2.0f*((float)_X/(float)(Size-1)-0.5f),Border*2.0f*((float)(Size-_Y-1)/(float)(Size-1)-0.5f));
13.    float2 z=(float2)(0,0);
14.    for(uint i=0;(i<NumOfIter) && (length(z)<Threshold);i++)
15.      z=Mult(Mult(Mult(Mult(z,z),z),z),z)+c;
16.    if(length(z)<Threshold)
17.      return 1;
18.    else
19.      return 0;
20. }
21.  
22.  
23. __kernel void CalcMandel(float _Border,float _Threshold,uint _NumOfIter,uint _Size,__global uchar *Data/*,__global uint *NumOfPoints*/)
24. {
25.    Threshold=_Threshold;
26.    Size  =_Size;
27.    Border  =_Border;
28.   NumOfIter=_NumOfIter;
29.    size_t _X=get_global_id(0);
30.    size_t _Y=get_global_id(1);
31.    Data[_X+_Size*_Y]=Mandel(_X,_Y);
32.    if((Mandel(_X-1,_Y)==1)&&(Mandel(_X+1,_Y)==1))
33.      if((Mandel(_X-1,_Y-1)==1)&&(Mandel(_X,_Y-1)==1)&&(Mandel(_X+1,_Y-1)==1))
34.        if((Mandel(_X-1,_Y+1)==1)&&(Mandel(_X,_Y+1)==1)&&(Mandel(_X+1,_Y+1)==1))
35.          Data[_X+_Size*_Y]=0;
36. }

Я вот думаю val - это ведь private memory. И при синхронизации потоков в GPU просто не хватает памяти, чтобы сохранить значение переменной для каждого потока. Если сохранять промежуточные вычисления в __global массиве, то всё должно работать.

 

Так не должно быть, это - глюк.
В вашем варианте производится больше вычислений.