1 /**
2 * Mandelbulber v2, a 3D fractal generator ,=#MKNmMMKmmßMNWy,
3 * ,B" ]L,,p%%%,,,§;, "K
4 * Copyright (C) 2017-21 Mandelbulber Team §R-==%w["'~5]m%=L.=~5N
5 * ,=mm=§M ]=4 yJKA"/-Nsaj "Bw,==,,
6 * This file is part of Mandelbulber. §R.r= jw",M Km .mM FW ",§=ß., ,TN
7 * ,4R =%["w[N=7]J '"5=],""]]M,w,-; T=]M
8 * Mandelbulber is free software: §R.ß~-Q/M=,=5"v"]=Qf,'§"M= =,M.§ Rz]M"Kw
9 * you can redistribute it and/or §w "xDY.J ' -"m=====WeC=\ ""%""y=%"]"" §
10 * modify it under the terms of the "§M=M =D=4"N #"%==A%p M§ M6 R' #"=~.4M
11 * GNU General Public License as §W =, ][T"]C § § '§ e===~ U !§[Z ]N
12 * published by the 4M",,Jm=,"=e~ § § j]]""N BmM"py=ßM
13 * Free Software Foundation, ]§ T,M=& 'YmMMpM9MMM%=w=,,=MT]M m§;'§,
14 * either version 3 of the License, TWw [.j"5=~N[=§%=%W,T ]R,"=="Y[LFT ]N
15 * or (at your option) TW=,-#"%=;[ =Q:["V"" ],,M.m == ]N
16 * any later version. J§"mr"] ,=,," =="""J]= M"M"]==ß"
17 * §= "=C=4 §"eM "=B:m|4"]#F,§~
18 * Mandelbulber is distributed in "9w=,,]w em%wJ '"~" ,=,,ß"
19 * the hope that it will be useful, . "K= ,=RMMMßM"""
20 * but WITHOUT ANY WARRANTY; .'''
21 * without even the implied warranty
22 * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
23 *
24 * See the GNU General Public License for more details.
25 * You should have received a copy of the GNU General Public License
26 * along with Mandelbulber. If not, see <http://www.gnu.org/licenses/>.
27 *
28 * ###########################################################################
29 *
30 * Authors: Krzysztof Marczak (buddhi1980@gmail.com), Sebastian Jennen (jenzebas@gmail.com)
31 *
32 * mathematical algebra code for opencl
33 */
34
35 #ifndef MANDELBULBER2_OPENCL_ALGEBRA_HPP_
36 #define MANDELBULBER2_OPENCL_ALGEBRA_HPP_
37
38 #ifndef OPENCL_KERNEL_CODE
39 #include <QString>
40
41 #include "src/algebra.hpp"
42 #include "src/color_structures.hpp"
43 #endif
44
45 typedef struct
46 {
47 cl_float3 m1;
48 cl_float3 m2;
49 cl_float3 m3;
50 } matrix33;
51
52 #ifndef OPENCL_KERNEL_CODE
toClMatrix33(CRotationMatrix source)53 inline matrix33 toClMatrix33(CRotationMatrix source)
54 {
55 CMatrix33 matrix = source.GetMatrix();
56 matrix33 m;
57 m.m1 = {{cl_float(matrix.m11), cl_float(matrix.m12), cl_float(matrix.m13), cl_float(0.0)}};
58 m.m2 = {{cl_float(matrix.m21), cl_float(matrix.m22), cl_float(matrix.m23), cl_float(0.0)}};
59 m.m3 = {{cl_float(matrix.m31), cl_float(matrix.m32), cl_float(matrix.m33), cl_float(0.0)}};
60 return m;
61 }
toClFloat2(CVector2<double> v)62 inline cl_float2 toClFloat2(CVector2<double> v)
63 {
64 cl_float2 retVal = {{cl_float(v.x), cl_float(v.y)}};
65 return retVal;
66 }
toClFloat3(CVector3 v)67 inline cl_float3 toClFloat3(CVector3 v)
68 {
69 cl_float3 retVal = {{cl_float(v.x), cl_float(v.y), cl_float(v.z), cl_float(0.0)}};
70 return retVal;
71 }
toClInt3(int x,int y,int z)72 inline cl_int3 toClInt3(int x, int y, int z)
73 {
74 cl_int3 retVal = {{cl_int(x), cl_int(y), cl_int(z), cl_int(0)}};
75 return retVal;
76 }
toClInt3(sRGB c)77 inline cl_int3 toClInt3(sRGB c)
78 {
79 cl_int3 retVal = {{cl_int(c.R), cl_int(c.G), cl_int(c.B), cl_int(1)}};
80 return retVal;
81 }
toClFloat3(sRGBFloat c)82 inline cl_float3 toClFloat3(sRGBFloat c)
83 {
84 cl_float3 retVal = {{cl_float(c.R), cl_float(c.G), cl_float(c.B), cl_float(1.0f)}};
85 return retVal;
86 }
toClFloat4(CVector4 v)87 inline cl_float4 toClFloat4(CVector4 v)
88 {
89 cl_float4 retVal = {{cl_float(v.x), cl_float(v.y), cl_float(v.z), cl_float(v.w)}};
90 return retVal;
91 }
toClFloat4(sRGBAfloat v)92 inline cl_float4 toClFloat4(sRGBAfloat v)
93 {
94 cl_float4 retVal = {{cl_float(v.R), cl_float(v.G), cl_float(v.B), cl_float(v.A)}};
95 return retVal;
96 }
97
98 #endif
99
100 #ifdef OPENCL_KERNEL_CODE
101
Matrix33MulFloat4(matrix33 matrix,float4 vect)102 inline float4 Matrix33MulFloat4(matrix33 matrix, float4 vect)
103 {
104 float4 out;
105 out.x = dot(vect.xyz, matrix.m1);
106 out.y = dot(vect.xyz, matrix.m2);
107 out.z = dot(vect.xyz, matrix.m3);
108 out.w = vect.w;
109 return out;
110 }
111
Matrix33MulFloat3(matrix33 matrix,float3 vect)112 inline float3 Matrix33MulFloat3(matrix33 matrix, float3 vect)
113 {
114 float3 out;
115 out.x = dot(vect.xyz, matrix.m1);
116 out.y = dot(vect.xyz, matrix.m2);
117 out.z = dot(vect.xyz, matrix.m3);
118 return out;
119 }
120
Matrix33MulMatrix33(matrix33 m1,matrix33 m2)121 matrix33 Matrix33MulMatrix33(matrix33 m1, matrix33 m2)
122 {
123 matrix33 out;
124 out.m1.x = m1.m1.x * m2.m1.x + m1.m1.y * m2.m2.x + m1.m1.z * m2.m3.x;
125 out.m1.y = m1.m1.x * m2.m1.y + m1.m1.y * m2.m2.y + m1.m1.z * m2.m3.y;
126 out.m1.z = m1.m1.x * m2.m1.z + m1.m1.y * m2.m2.z + m1.m1.z * m2.m3.z;
127 out.m2.x = m1.m2.x * m2.m1.x + m1.m2.y * m2.m2.x + m1.m2.z * m2.m3.x;
128 out.m2.y = m1.m2.x * m2.m1.y + m1.m2.y * m2.m2.y + m1.m2.z * m2.m3.y;
129 out.m2.z = m1.m2.x * m2.m1.z + m1.m2.y * m2.m2.z + m1.m2.z * m2.m3.z;
130 out.m3.x = m1.m3.x * m2.m1.x + m1.m3.y * m2.m2.x + m1.m3.z * m2.m3.x;
131 out.m3.y = m1.m3.x * m2.m1.y + m1.m3.y * m2.m2.y + m1.m3.z * m2.m3.y;
132 out.m3.z = m1.m3.x * m2.m1.z + m1.m3.y * m2.m2.z + m1.m3.z * m2.m3.z;
133 return out;
134 }
135
RotateX(matrix33 m,float angle)136 matrix33 RotateX(matrix33 m, float angle)
137 {
138 matrix33 out, rot;
139 float s = sin(angle);
140 float c = cos(angle);
141 rot.m1 = (float3){1.0f, 0.0f, 0.0f};
142 rot.m2 = (float3){0.0f, c, -s};
143 rot.m3 = (float3){0.0f, s, c};
144 out = Matrix33MulMatrix33(m, rot);
145 return out;
146 }
147
RotateY(matrix33 m,float angle)148 matrix33 RotateY(matrix33 m, float angle)
149 {
150 matrix33 out, rot;
151 float s = sin(angle);
152 float c = cos(angle);
153 rot.m1 = (float3){c, 0.0f, s};
154 rot.m2 = (float3){0.0f, 1.0f, 0.0f};
155 rot.m3 = (float3){-s, 0.0f, c};
156 out = Matrix33MulMatrix33(m, rot);
157 return out;
158 }
159
RotateZ(matrix33 m,float angle)160 matrix33 RotateZ(matrix33 m, float angle)
161 {
162 matrix33 out, rot;
163 float s = sin(angle);
164 float c = cos(angle);
165 rot.m1 = (float3){c, -s, 0.0f};
166 rot.m2 = (float3){s, c, 0.0f};
167 rot.m3 = (float3){0.0f, 0.0f, 1.0f};
168 out = Matrix33MulMatrix33(m, rot);
169 return out;
170 }
171
RotateAroundVectorByAngle(float3 origin,float3 axis,float angle)172 float3 RotateAroundVectorByAngle(float3 origin, float3 axis, float angle)
173 {
174 float3 vector = origin * cos(angle);
175 vector += cross(axis, origin) * sin(angle);
176 vector += axis * dot(axis, origin) * (1.0f - cos(angle));
177 return vector;
178 }
179
RotateAroundVectorByAngle4(float4 origin4d,float3 axis,float angle)180 float4 RotateAroundVectorByAngle4(float4 origin4d, float3 axis, float angle)
181 {
182 float3 origin = origin4d.xyz;
183 float3 vector = origin * cos(angle);
184 vector += cross(axis, origin) * sin(angle);
185 vector += axis * dot(axis, origin) * (1.0f - cos(angle));
186 return (float4){vector.x, vector.y, vector.z, origin4d.w};
187 }
188
TransposeMatrix(matrix33 m)189 inline matrix33 TransposeMatrix(matrix33 m)
190 {
191 matrix33 out;
192 out.m1 = (float3){m.m1.x, m.m2.x, m.m3.x};
193 out.m2 = (float3){m.m1.y, m.m2.y, m.m3.y};
194 out.m3 = (float3){m.m1.z, m.m2.z, m.m3.z};
195 return out;
196 }
197
SmoothConditionAGreaterB(float a,float b,float sharpness)198 float SmoothConditionAGreaterB(float a, float b, float sharpness)
199 {
200 return native_recip(1.0f + native_exp(sharpness * (b - a)));
201 }
202
SmoothConditionALessB(float a,float b,float sharpness)203 float SmoothConditionALessB(float a, float b, float sharpness)
204 {
205 return native_recip(1.0f + native_exp(sharpness * (a - b)));
206 }
207
wrap(float3 x,float3 a,float3 s)208 inline float3 wrap(float3 x, float3 a, float3 s)
209 {
210 x -= s;
211 float3 out;
212 out.x = x.x - a.x * floor(x.x / a.x) + s.x;
213 out.y = x.y - a.y * floor(x.y / a.y) + s.y;
214 out.z = x.z - a.z * floor(x.z / a.z) + s.z;
215 return out;
216 }
217
218 //********** Random ******************************
RandomInt(int * randomSeed)219 int RandomInt(int *randomSeed)
220 {
221 int s = *randomSeed;
222 int i = 0;
223 int const a = 15484817;
224 int const m = 6571759;
225 s = ((long)(s * a)) % m;
226 if (s < 0) s = -s;
227 return *randomSeed = s;
228 }
229
Random(int max,int * randomSeed)230 int Random(int max, int *randomSeed)
231 {
232 return RandomInt(randomSeed) % (max + 1);
233 }
234
vectorMod(float3 vector1,float3 vector2)235 inline float3 vectorMod(float3 vector1, float3 vector2)
236 {
237 return (float3){(vector2.x > 0.0f ? fmod(vector1.x, vector2.x) : vector1.x),
238 (vector2.y > 0.0f ? fmod(vector1.y, vector2.y) : vector1.y),
239 (vector2.z > 0.0f ? fmod(vector1.z, vector2.z) : vector1.z)};
240 }
241
modRepeat(float3 vector1,float3 repeat)242 inline float3 modRepeat(float3 vector1, float3 repeat)
243 {
244 if (length(repeat) == 0.0f) return vector1;
245 return vectorMod((vectorMod((vector1 - repeat * 0.5f), repeat) + repeat), repeat) - repeat * 0.5f;
246 }
247
248 #ifdef ITERATION_WEIGHT
SmoothCVector(const float4 v1,const float4 v2,float k)249 float4 SmoothCVector(const float4 v1, const float4 v2, float k)
250 {
251 float4 result;
252 float nk = 1.0f - k;
253
254 if (k <= 0.0f)
255 {
256 result = v1;
257 }
258 else if (k >= 1.0f)
259 {
260 result = v2;
261 }
262 else
263 {
264 float length1 = length(v1);
265 float length2 = length(v2);
266 float lenInterp = length1 * nk + length2 * k;
267 float4 vTemp = v1 * nk + v2 * k;
268 float lengthTemp = length(vTemp);
269 if (lengthTemp > 0.0f)
270 {
271 result = (vTemp / lengthTemp) * lenInterp;
272 }
273 else
274 {
275 result = v1;
276 }
277 }
278 return result;
279 }
280 #endif
281
LengthPow(float2 vect,float p)282 inline float LengthPow(float2 vect, float p)
283 {
284 return pow(pow(vect.x, p) + pow(vect.y, p), 1.0f / p);
285 }
286
GetAlpha(float3 vect)287 inline float GetAlpha(float3 vect)
288 {
289 return atan2(vect.y, vect.x);
290 }
291
GetBeta(float3 vect)292 inline float GetBeta(float3 vect)
293 {
294 return atan2(vect.z, length(vect.xy));
295 }
296
297 #endif
298
299 #endif // MANDELBULBER2_OPENCL_ALGEBRA_HPP_
300