1 /+ 2 Vasaro Copyright © 2018 Andrea Fontana 3 This file is part of Vasaro. 4 5 Vasaro is free software: you can redistribute it and/or modify 6 it under the terms of the GNU General Public License as published by 7 the Free Software Foundation, either version 3 of the License, or 8 (at your option) any later version. 9 10 Vasaro is distributed in the hope that it will be useful, 11 but WITHOUT ANY WARRANTY; without even the implied warranty of 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 GNU General Public License for more details. 14 15 You should have received a copy of the GNU General Public License 16 along with Vasaro. If not, see <http://www.gnu.org/licenses/>. 17 +/ 18 module generator; 19 20 import std.concurrency; 21 import std.conv : to; 22 23 import opensimplexnoise; 24 import vec3; 25 import viewer; 26 27 public: 28 29 __gshared Tid generatorTid; 30 __gshared size_t currentModel = 0; // 0 or 1 31 __gshared float minDiameter = 60; 32 __gshared float maxDiameter = 100; 33 __gshared int resolution = 300; 34 __gshared float vaseHeight = 100; 35 __gshared float layerHeight = 0.2; 36 37 __gshared float[10] vaseProfileCheckPoints = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]; 38 39 struct Noise 40 { 41 long seed = 0; 42 float xScale = 3; 43 float yScale = 3; 44 float amplitude = 2; 45 float rotation = 0; 46 float[10] strengthPoints = [1,1,1,1,1,1,1,1,1,1]; 47 bool visible = true; 48 } 49 50 __gshared Noise[] noises; 51 52 void start() 53 { 54 if (gs_actual != 0) 55 return; 56 57 generatorTid = spawn(&initGenerator); 58 generatorTid.setMaxMailboxSize(1, OnCrowding.ignore); // Avoid spam 59 60 debug { 61 import std.stdio : writeln; 62 writeln("Generator started."); 63 } 64 } 65 66 void build() 67 { 68 import core.thread; 69 gs_requested++; 70 send(generatorTid, true); 71 } 72 73 bool isRunning() { return true; /*gs_actual != GeneratorStatus.STOPPED;*/ } 74 bool isCompleted() { return gs_actual == gs_requested; /*gs_actual == GeneratorStatus.COMPLETED;*/ } 75 76 size_t lastGenerated() { return gs_actual; } 77 78 private: 79 80 __gshared size_t gs_requested = 1; 81 __gshared size_t gs_actual = 0; 82 83 import std.datetime.stopwatch : StopWatch; 84 85 void initGenerator() 86 { 87 try 88 { 89 // Sync with other threads 90 while(true) 91 { 92 bool received = receiveOnly!bool; 93 if (gs_actual != gs_requested) createVase(); 94 } 95 } 96 97 // That's ok, app was closed. 98 catch (OwnerTerminated ot) { debug { import std.stdio; writeln("Generator cleared."); } } 99 } 100 101 102 void createVase() 103 { 104 import std.stdio; 105 import std.typecons : Tuple, tuple; 106 import std.array : uninitializedArray; 107 import std.parallelism : parallel; 108 import std.range : iota; 109 import std.math : sin, cos, PI; 110 import std.datetime.stopwatch : StopWatch; 111 112 113 // Calculate spline coeff. from ten points equally-separated 114 float[4][10] naturalSpline(float[10] y) 115 { 116 float[4][10] results; 117 float[9] mu; 118 float[10] z; 119 float g = 0; 120 121 mu[0] = 0; 122 z[0] = 0; 123 124 for (int i = 1; i < 9; i++) { 125 g = 4 - mu[i -1]; 126 mu[i] = 1 / g; 127 z[i] = (3 * (y[i + 1] - y[i] * 2 + y[i - 1]) - z[i - 1]) / g; 128 } 129 130 float[9] b; 131 float[10] c; 132 float[9] d; 133 134 z[9] = 0; 135 c[9] = 0; 136 137 for (int j = 8; j >=0; j--) { 138 c[j] = z[j] - mu[j] * c[j + 1]; 139 b[j] = (y[j + 1] - y[j]) - 1 * (c[j + 1] + 2 * c[j]) / 3; 140 d[j] = (c[j + 1] - c[j]) / 3; 141 } 142 143 for (int i = 0; i < 9; i++) { 144 results[i][0] = y[i]; 145 results[i][1] = b[i]; 146 results[i][2] = c[i]; 147 results[i][3] = d[i]; 148 } 149 150 return results; 151 } 152 153 // Using coeff. calculated above we can approximate y value for x. 154 // Used to interpolate points set by user 155 float interpolateSpline(float[4][10] coeff, float x) 156 { 157 size_t xx = x.to!int; 158 float val = x-xx; 159 auto cur = coeff[xx]; 160 return cur[0] + cur[1]*val + cur[2]*val*val + cur[3]*val*val*val; 161 } 162 163 // Some benchmark 164 StopWatch swGen; 165 swGen.start(); 166 167 auto noisesCopy = noises; 168 169 size_t current = gs_requested; // Asked revision 170 size_t candidateModel = (currentModel+1)%2; // Buffer we're working on. No errors? It will become the rendered model. 171 172 // If another request is queued, we can interrupt this one. 173 bool continueProcessing() { return (gs_requested == current); } 174 175 float min_diameter = minDiameter; 176 float max_diameter = maxDiameter; 177 int res = resolution; 178 179 float height = vaseHeight; 180 float layer = layerHeight; 181 182 auto layersCnt = (height/layer).to!size_t; 183 184 Vec3[][] sideMeshVertex = uninitializedArray!(Vec3[][])(res,layersCnt); 185 Vec3[][] sideMeshVertexNormals = new Vec3[][](res,layersCnt); 186 187 alias Coords = Tuple!(size_t, size_t); 188 Coords[] sideMeshVertexNormalsMap; 189 190 float xC = 0.0f; 191 192 193 // 194 // RADIUS VARIANCE 195 // 196 197 float[] layersRadiusFactor = uninitializedArray!(float[])(layersCnt); 198 199 { 200 // I want to limit function inside [-1;1] interval 201 float maxDiameterFactor = 1; 202 float minDiameterFactor = -1; 203 204 205 auto layerSplineCoeff = naturalSpline(vaseProfileCheckPoints); 206 foreach(i, ref lf; layersRadiusFactor) 207 { 208 lf = interpolateSpline(layerSplineCoeff, 9.0f*i*1.0f/layersCnt)*2-1; 209 if (lf < minDiameterFactor) minDiameterFactor = lf; 210 else if (lf > maxDiameterFactor) maxDiameterFactor = lf; 211 } 212 } 213 214 // 215 // NOISE 216 // 217 218 auto noiseMesh = new float[][](res, layersCnt); 219 220 immutable meanDiameter = (min_diameter+max_diameter)/2; 221 immutable diameterDelta = (max_diameter-min_diameter)/2; 222 223 double min = double.max; 224 double max = -double.max; 225 226 227 float[4][10][] noiseCoeff; 228 noiseCoeff.length = noisesCopy.length; 229 230 OpenSimplexNoise!float[] sn; 231 sn.length = noisesCopy.length; 232 233 // Noise strength 234 foreach(i, ref n; noisesCopy) 235 { 236 if (!n.visible) continue; 237 sn[i] = new OpenSimplexNoise!float(n.seed); 238 noiseCoeff[i] = naturalSpline(n.strengthPoints); 239 } 240 241 242 auto circumference = meanDiameter * PI; 243 bool hasNoise = false; 244 245 // Calculate noise for each point of mesh 246 foreach(i, ref n; noisesCopy) 247 { 248 if (!n.visible) continue; 249 250 auto twist = n.rotation/layersCnt; 251 252 immutable xScaleNorm = n.xScale*circumference/200; 253 immutable yScaleNorm = 4*n.yScale*height/200; 254 255 foreach(size_t x; parallel(iota(cast(size_t)0, res))) 256 { 257 auto xx = cast(float)(x)/(res-1); 258 auto noiseMeshPtr = &noiseMesh[x][0]; 259 260 if (!continueProcessing()) continue; 261 262 foreach(size_t y; 0 .. layersCnt) 263 { 264 auto yy = cast(float)(y)/(layersCnt-1); 265 float noise = sn[i].eval( 266 cos(2*PI*(xx+twist*y))*xScaleNorm, 267 sin(2*PI*(xx+twist*y))*xScaleNorm, 268 yy*yScaleNorm 269 ) 270 *n.amplitude 271 *interpolateSpline(noiseCoeff[i], 9.0f*y/layersCnt)*2-1; 272 273 if (!hasNoise) noiseMeshPtr[y] = noise; 274 else noiseMeshPtr[y] += noise; 275 276 if (noise > max) max = noise; 277 else if (noise< min) min = noise; 278 } 279 } 280 281 hasNoise = true; 282 } 283 284 285 // Sum all the things up 286 foreach(size_t x; parallel(iota(cast(size_t)0, res))) 287 { 288 if (!continueProcessing()) continue; 289 auto sideMeshVertexPtr = &sideMeshVertex[x][0]; 290 291 foreach(size_t y; (iota(0,layersCnt))) 292 { 293 immutable curWidth = meanDiameter + diameterDelta * layersRadiusFactor[y]; 294 sideMeshVertexPtr[y] = Vec3((curWidth/2)*cos(2*PI/(res-1)*x),y*layer,(curWidth/2)*sin(2*PI/(res-1)*x)); 295 296 if (hasNoise) 297 { 298 immutable noise = (noiseMesh[x][y] - min); 299 sideMeshVertexPtr[y].x += noise * cos(2*PI/(res-1)*x); 300 sideMeshVertexPtr[y].z += noise * sin(2*PI/(res-1)*x); 301 } 302 303 sideMeshVertexNormals[x][y] = Vec3(0,0,0); 304 } 305 } 306 307 if (!continueProcessing) return; 308 309 // Ok, we can guess how many triangle we need 310 immutable size_t side_vertex_coords_count = 311 (layersCnt-2) // Number of layers 312 * (res-1) // Number of points per layer 313 * 3 // 3 coords for each point 314 * 3 // 3 points for each triangle 315 * 2 // 2 triangle for each point 316 317 + 2 // Last two layers 318 * (res-1) // Number of points per layer 319 * 3 // 3 coords for each point 320 * 3 // 3 points for each triangle 321 * 1; // 1 triangle for each point 0 322 323 immutable size_t base_vertex_coords_count = 324 (res-1) // Triangles 325 * 3 // 3 coords for each point 326 * 3; // 3 points for each triangle 327 328 immutable size_t total_vertex_coords_count = 329 side_vertex_coords_count 330 + 2 // Top and bottom layer 331 * base_vertex_coords_count; 332 333 334 model[candidateModel].vertex.reserve(total_vertex_coords_count); 335 sideMeshVertexNormals.reserve(total_vertex_coords_count); 336 sideMeshVertexNormalsMap.reserve(total_vertex_coords_count); 337 338 model[candidateModel].vertex.length = side_vertex_coords_count; // For each vertex we store x,y,z 339 sideMeshVertexNormalsMap.length = side_vertex_coords_count / 3; // For each vertex we store a normal 340 341 import core.atomic; 342 shared(size_t) globalParallelIdx = 0; 343 344 // Parallelization: to avoid multithread problems first we do all even rows, then all the odd. 345 foreach(idx; 0..2) 346 foreach (size_t x; parallel(iota(size_t(idx),res-1,2))) 347 { 348 if (!continueProcessing()) continue; 349 350 foreach(size_t y; 0..layersCnt) 351 { 352 353 if (y < layersCnt-1) 354 { 355 immutable parallelIdx = globalParallelIdx.atomicOp!"+="(1) - 1; 356 auto vertexSlice = model[candidateModel].vertex[parallelIdx*9+0..parallelIdx*9+9]; 357 auto sideMeshVertexNormalsMapSlice = sideMeshVertexNormalsMap[parallelIdx*3..parallelIdx*3+3]; 358 359 const cur = &sideMeshVertex[x][y]; 360 const right = &sideMeshVertex[x+1][y]; 361 const top = &sideMeshVertex[x][y+1]; 362 363 immutable normal = (*top-*cur).crossProduct(*right-*cur); 364 365 sideMeshVertexNormals[x][y] += normal; 366 sideMeshVertexNormals[x+1][y] += normal; 367 368 if (x+1 == res-1) sideMeshVertexNormals[0][y] += normal; 369 370 sideMeshVertexNormals[x][y+1] += normal; 371 372 vertexSlice[0..3] = cur.data[0..3]; 373 vertexSlice[3..6] = top.data[0..3]; 374 vertexSlice[6..9] = right.data[0..3]; 375 376 sideMeshVertexNormalsMapSlice[0] = tuple(x,y); 377 sideMeshVertexNormalsMapSlice[1] = tuple(x+1,y); 378 sideMeshVertexNormalsMapSlice[2] = tuple(x,y+1); 379 } 380 381 382 if (y > 0) 383 { 384 immutable parallelIdx = globalParallelIdx.atomicOp!"+="(1) - 1; 385 auto vertexSlice = model[candidateModel].vertex[parallelIdx*9+0..parallelIdx*9+9]; 386 auto sideMeshVertexNormalsMapSlice = sideMeshVertexNormalsMap[parallelIdx*3..parallelIdx*3+3]; 387 388 auto cur = &sideMeshVertex[x][y]; 389 auto right = &sideMeshVertex[x+1][y]; 390 auto bottom = &sideMeshVertex[x+1][y-1]; 391 392 immutable normal = (*right-*cur).crossProduct(*bottom-*cur); 393 394 sideMeshVertexNormals[x][y] += normal; 395 sideMeshVertexNormals[x+1][y] += normal; 396 sideMeshVertexNormals[x+1][y-1] += normal; 397 398 if (x+1 == res-1) 399 { 400 sideMeshVertexNormals[0][y] += normal; 401 sideMeshVertexNormals[0][y-1] += normal; 402 } 403 404 vertexSlice[0..3] = cur.data[0..3]; 405 vertexSlice[3..6] = right.data[0..3]; 406 vertexSlice[6..9] = bottom.data[0..3]; 407 408 sideMeshVertexNormalsMapSlice[0] = tuple(x,y); 409 sideMeshVertexNormalsMapSlice[1] = tuple(x+1,y); 410 sideMeshVertexNormalsMapSlice[2] = tuple(x+1,y-1); 411 } 412 413 } 414 } 415 416 if (!continueProcessing) return; 417 418 // Top and bottom base 419 for(size_t x = 1; x < res; ++x) 420 { 421 if (!continueProcessing()) return; 422 423 auto b = &sideMeshVertex[x-1][0]; 424 auto c = &sideMeshVertex[x][0]; 425 426 model[candidateModel].vertex ~= [0, 0, 0]; // All triangles have a vertex in base center 427 model[candidateModel].vertex ~= b.data[0..3]; 428 model[candidateModel].vertex ~= c.data[0..3]; 429 } 430 431 for(size_t x = 1; x < res; ++x) 432 { 433 if (!continueProcessing()) return; 434 435 auto b = &sideMeshVertex[x][layersCnt-1]; 436 auto c = &sideMeshVertex[x-1][layersCnt-1]; 437 438 model[candidateModel].vertex ~= [0, layer*(layersCnt-1), 0]; 439 model[candidateModel].vertex ~= b.data[0..3]; 440 model[candidateModel].vertex ~= c.data[0..3]; 441 } 442 443 // Sum up all the normals 444 model[candidateModel].vertexNormals.length = model[candidateModel].vertex.length; 445 foreach(i, nIdx; parallel(sideMeshVertexNormalsMap)) 446 { 447 if (!continueProcessing()) continue; 448 449 // Normals of side mesh 450 if (i < side_vertex_coords_count / 3) model[candidateModel].vertexNormals[i*3..i*3+3] = sideMeshVertexNormals[nIdx[0]][nIdx[1]].normalized().data[0..3]; 451 } 452 453 // Normals of two bases 454 size_t startingIdx = sideMeshVertexNormalsMap.length; 455 foreach(i; parallel(iota(startingIdx,startingIdx+(res-1)*3))) 456 model[candidateModel].vertexNormals[i*3..i*3+3] = [0, -1, 0]; 457 458 startingIdx += (res-1)*3; 459 foreach(i; parallel(iota(startingIdx,startingIdx+(res-1)*3))) 460 model[candidateModel].vertexNormals[i*3..i*3+3] = [0, 1, 0]; 461 462 swGen.stop(); 463 464 if (!continueProcessing()) return; 465 466 gs_actual = current; 467 currentModel = candidateModel; 468 writeln("Vase regenerated in: ", swGen.peek.total!"msecs", "ms"); 469 }