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 import std.experimental.allocator; 112 import std.experimental.allocator.mallocator : Mallocator; 113 114 // Calculate spline coeff. from ten points equally-separated 115 float[4][10] naturalSpline(float[10] y) 116 { 117 float[4][10] results; 118 float[9] mu; 119 float[10] z; 120 float g = 0; 121 122 mu[0] = 0; 123 z[0] = 0; 124 125 for (int i = 1; i < 9; i++) { 126 g = 4 - mu[i -1]; 127 mu[i] = 1 / g; 128 z[i] = (3 * (y[i + 1] - y[i] * 2 + y[i - 1]) - z[i - 1]) / g; 129 } 130 131 float[9] b; 132 float[10] c; 133 float[9] d; 134 135 z[9] = 0; 136 c[9] = 0; 137 138 for (int j = 8; j >=0; j--) { 139 c[j] = z[j] - mu[j] * c[j + 1]; 140 b[j] = (y[j + 1] - y[j]) - 1 * (c[j + 1] + 2 * c[j]) / 3; 141 d[j] = (c[j + 1] - c[j]) / 3; 142 } 143 144 for (int i = 0; i < 9; i++) { 145 results[i][0] = y[i]; 146 results[i][1] = b[i]; 147 results[i][2] = c[i]; 148 results[i][3] = d[i]; 149 } 150 151 return results; 152 } 153 154 // Using coeff. calculated above we can approximate y value for x. 155 // Used to interpolate points set by user 156 float interpolateSpline(float[4][10] coeff, float x) 157 { 158 size_t xx = x.to!int; 159 float val = x-xx; 160 auto cur = coeff[xx]; 161 return cur[0] + cur[1]*val + cur[2]*val*val + cur[3]*val*val*val; 162 } 163 164 // Some benchmark 165 StopWatch swGen; 166 swGen.start(); 167 168 auto noisesCopy = noises; 169 170 size_t current = gs_requested; // Asked revision 171 size_t candidateModel = (currentModel+1)%2; // Buffer we're working on. No errors? It will become the rendered model. 172 173 // If another request is queued, we can interrupt this one. 174 bool continueProcessing() { return (gs_requested == current); } 175 176 float min_diameter = minDiameter; 177 float max_diameter = maxDiameter; 178 int res = resolution; 179 180 float height = vaseHeight; 181 float layer = layerHeight; 182 183 auto layersCnt = (height/layer).to!size_t; 184 185 Vec3[][] sideMeshVertex = Mallocator.instance.makeMultidimensionalArray!Vec3(res, layersCnt); 186 Vec3[][] sideMeshVertexNormals = Mallocator.instance.makeMultidimensionalArray!Vec3(res, layersCnt); 187 188 scope(exit) 189 { 190 Mallocator.instance.disposeMultidimensionalArray(sideMeshVertex); 191 Mallocator.instance.disposeMultidimensionalArray(sideMeshVertexNormals); 192 } 193 194 struct Coords { size_t x; size_t y; } 195 Coords[] sideMeshVertexNormalsMap = null; 196 197 float xC = 0.0f; 198 199 200 // 201 // RADIUS VARIANCE 202 // 203 204 float[] layersRadiusFactor = Mallocator.instance.makeArray!float(layersCnt); 205 scope(exit) Mallocator.instance.dispose(layersRadiusFactor); 206 207 { 208 // I want to limit function inside [-1;1] interval 209 float maxDiameterFactor = 1; 210 float minDiameterFactor = -1; 211 212 213 auto layerSplineCoeff = naturalSpline(vaseProfileCheckPoints); 214 foreach(i, ref lf; layersRadiusFactor) 215 { 216 lf = interpolateSpline(layerSplineCoeff, 9.0f*i*1.0f/layersCnt)*2-1; 217 if (lf < minDiameterFactor) minDiameterFactor = lf; 218 else if (lf > maxDiameterFactor) maxDiameterFactor = lf; 219 } 220 } 221 222 // 223 // NOISE 224 // 225 226 auto noiseMesh = Mallocator.instance.makeMultidimensionalArray!float(res, layersCnt); 227 scope(exit) Mallocator.instance.disposeMultidimensionalArray(noiseMesh); 228 229 immutable meanDiameter = (min_diameter+max_diameter)/2; 230 immutable diameterDelta = (max_diameter-min_diameter)/2; 231 232 double min = double.max; 233 double max = -double.max; 234 235 236 float[4][10][] noiseCoeff; 237 noiseCoeff.length = noisesCopy.length; 238 239 OpenSimplexNoise!float[] sn; 240 sn.length = noisesCopy.length; 241 242 // Noise strength 243 foreach(i, ref n; noisesCopy) 244 { 245 if (!n.visible) continue; 246 sn[i] = new OpenSimplexNoise!float(n.seed); 247 noiseCoeff[i] = naturalSpline(n.strengthPoints); 248 } 249 250 251 auto circumference = meanDiameter * PI; 252 bool hasNoise = false; 253 254 // Calculate noise for each point of mesh 255 foreach(i, ref n; noisesCopy) 256 { 257 if (!n.visible) continue; 258 259 auto twist = n.rotation/layersCnt; 260 261 immutable xScaleNorm = n.xScale*circumference/200; 262 immutable yScaleNorm = 4*n.yScale*height/200; 263 264 foreach(size_t x; parallel(iota(cast(size_t)0, res))) 265 { 266 auto xx = cast(float)(x)/(res-1); 267 auto noiseMeshPtr = &noiseMesh[x][0]; 268 269 if (!continueProcessing()) continue; 270 271 foreach(size_t y; 0 .. layersCnt) 272 { 273 auto yy = cast(float)(y)/(layersCnt-1); 274 float noise = sn[i].eval( 275 cos(2*PI*(xx+twist*y))*xScaleNorm, 276 sin(2*PI*(xx+twist*y))*xScaleNorm, 277 yy*yScaleNorm 278 ) 279 *n.amplitude 280 *interpolateSpline(noiseCoeff[i], 9.0f*y/layersCnt)*2-1; 281 282 if (!hasNoise) noiseMeshPtr[y] = noise; 283 else noiseMeshPtr[y] += noise; 284 285 if (noise > max) max = noise; 286 else if (noise< min) min = noise; 287 } 288 } 289 290 hasNoise = true; 291 } 292 293 // Sum all the things up 294 foreach(size_t x; parallel(iota(cast(size_t)0, res))) 295 { 296 if (!continueProcessing()) continue; 297 auto sideMeshVertexPtr = &sideMeshVertex[x][0]; 298 299 foreach(size_t y; (iota(0,layersCnt))) 300 { 301 immutable curWidth = meanDiameter + diameterDelta * layersRadiusFactor[y]; 302 sideMeshVertexPtr[y] = Vec3((curWidth/2)*cos(2*PI/(res-1)*x),y*layer,(curWidth/2)*sin(2*PI/(res-1)*x)); 303 304 if (hasNoise) 305 { 306 immutable noise = (noiseMesh[x][y] - min); 307 sideMeshVertexPtr[y].x += noise * cos(2*PI/(res-1)*x); 308 sideMeshVertexPtr[y].z += noise * sin(2*PI/(res-1)*x); 309 } 310 311 sideMeshVertexNormals[x][y] = Vec3(0,0,0); 312 } 313 } 314 315 if (!continueProcessing) return; 316 317 // Ok, we can guess how many triangle we need 318 immutable size_t side_vertex_coords_count = 319 (layersCnt-2) // Number of layers 320 * (res-1) // Number of points per layer 321 * 3 // 3 coords for each point 322 * 3 // 3 points for each triangle 323 * 2 // 2 triangle for each point 324 325 + 2 // Last two layers 326 * (res-1) // Number of points per layer 327 * 3 // 3 coords for each point 328 * 3 // 3 points for each triangle 329 * 1; // 1 triangle for each point 0 330 331 immutable size_t base_vertex_coords_count = 332 (res-1) // Triangles 333 * 3 // 3 coords for each point 334 * 3; // 3 points for each triangle 335 336 immutable size_t total_vertex_coords_count = 337 side_vertex_coords_count 338 + 2 // Top and bottom layer 339 * base_vertex_coords_count; 340 341 if (model[candidateModel].vertex !is null) 342 Mallocator.instance.dispose(model[candidateModel].vertex); 343 344 if (model[candidateModel].vertexNormals !is null) 345 Mallocator.instance.dispose(model[candidateModel].vertexNormals); 346 347 model[candidateModel].vertex = Mallocator.instance.makeArray!float(total_vertex_coords_count); 348 model[candidateModel].vertexNormals = Mallocator.instance.makeArray!float(total_vertex_coords_count); 349 350 sideMeshVertexNormalsMap = Mallocator.instance.makeArray!Coords(side_vertex_coords_count/3); 351 scope(exit) Mallocator.instance.dispose(sideMeshVertexNormalsMap); 352 353 import core.atomic; 354 shared(size_t) globalParallelIdx = 0; 355 356 // Parallelization: to avoid multithread problems first we do all even rows, then all the odd. 357 foreach(idx; 0..2) 358 foreach (size_t x; parallel(iota(size_t(idx),res-1,2))) 359 { 360 if (!continueProcessing()) continue; 361 362 foreach(size_t y; 0..layersCnt) 363 { 364 365 if (y < layersCnt-1) 366 { 367 immutable parallelIdx = globalParallelIdx.atomicOp!"+="(1) - 1; 368 auto vertexSlice = model[candidateModel].vertex[parallelIdx*9+0..parallelIdx*9+9]; 369 auto sideMeshVertexNormalsMapSlice = sideMeshVertexNormalsMap[parallelIdx*3..parallelIdx*3+3]; 370 371 const cur = &sideMeshVertex[x][y]; 372 const right = &sideMeshVertex[x+1][y]; 373 const top = &sideMeshVertex[x][y+1]; 374 375 immutable normal = (*top-*cur).crossProduct(*right-*cur); 376 377 sideMeshVertexNormals[x][y] += normal; 378 sideMeshVertexNormals[x+1][y] += normal; 379 380 if (x+1 == res-1) sideMeshVertexNormals[0][y] += normal; 381 382 sideMeshVertexNormals[x][y+1] += normal; 383 384 vertexSlice[0..3] = cur.data[0..3]; 385 vertexSlice[3..6] = top.data[0..3]; 386 vertexSlice[6..9] = right.data[0..3]; 387 388 sideMeshVertexNormalsMapSlice[0] = Coords(x,y); 389 sideMeshVertexNormalsMapSlice[1] = Coords(x+1,y); 390 sideMeshVertexNormalsMapSlice[2] = Coords(x,y+1); 391 } 392 393 394 if (y > 0) 395 { 396 immutable parallelIdx = globalParallelIdx.atomicOp!"+="(1) - 1; 397 auto vertexSlice = model[candidateModel].vertex[parallelIdx*9+0..parallelIdx*9+9]; 398 auto sideMeshVertexNormalsMapSlice = sideMeshVertexNormalsMap[parallelIdx*3..parallelIdx*3+3]; 399 400 auto cur = &sideMeshVertex[x][y]; 401 auto right = &sideMeshVertex[x+1][y]; 402 auto bottom = &sideMeshVertex[x+1][y-1]; 403 404 immutable normal = (*right-*cur).crossProduct(*bottom-*cur); 405 406 sideMeshVertexNormals[x][y] += normal; 407 sideMeshVertexNormals[x+1][y] += normal; 408 sideMeshVertexNormals[x+1][y-1] += normal; 409 410 if (x+1 == res-1) 411 { 412 sideMeshVertexNormals[0][y] += normal; 413 sideMeshVertexNormals[0][y-1] += normal; 414 } 415 416 vertexSlice[0..3] = cur.data[0..3]; 417 vertexSlice[3..6] = right.data[0..3]; 418 vertexSlice[6..9] = bottom.data[0..3]; 419 420 sideMeshVertexNormalsMapSlice[0] = Coords(x,y); 421 sideMeshVertexNormalsMapSlice[1] = Coords(x+1,y); 422 sideMeshVertexNormalsMapSlice[2] = Coords(x+1,y-1); 423 } 424 425 } 426 } 427 428 if (!continueProcessing) return; 429 430 // Top and bottom base 431 { 432 size_t startingIdx = side_vertex_coords_count; 433 for(size_t x = 1; x < res; ++x) 434 { 435 if (!continueProcessing()) return; 436 437 auto b = &sideMeshVertex[x-1][0]; 438 auto c = &sideMeshVertex[x][0]; 439 440 model[candidateModel].vertex[startingIdx+0..startingIdx+3] = [0.0f, 0.0f, 0.0f]; // All triangles have a vertex in base center 441 model[candidateModel].vertex[startingIdx+3..startingIdx+6] = b.data[0..3]; 442 model[candidateModel].vertex[startingIdx+6..startingIdx+9] = c.data[0..3]; 443 startingIdx += 9; 444 } 445 446 for(size_t x = 1; x < res; ++x) 447 { 448 if (!continueProcessing()) return; 449 450 auto b = &sideMeshVertex[x][layersCnt-1]; 451 auto c = &sideMeshVertex[x-1][layersCnt-1]; 452 453 model[candidateModel].vertex[startingIdx+0..startingIdx+3] = [0.0f, layer*(layersCnt-1), 0.0f]; 454 model[candidateModel].vertex[startingIdx+3..startingIdx+6] = b.data[0..3]; 455 model[candidateModel].vertex[startingIdx+6..startingIdx+9] = c.data[0..3]; 456 startingIdx += 9; 457 } 458 } 459 460 // Sum up all the normals 461 foreach(i, nIdx; parallel(sideMeshVertexNormalsMap)) 462 { 463 if (!continueProcessing()) continue; 464 465 // Normals of side mesh 466 if (i < side_vertex_coords_count / 3) model[candidateModel].vertexNormals[i*3..i*3+3] = sideMeshVertexNormals[nIdx.x][nIdx.y].normalized().data[0..3]; 467 } 468 469 // Normals of two bases 470 { 471 size_t startingIdx = sideMeshVertexNormalsMap.length; 472 foreach(i; parallel(iota(startingIdx,startingIdx+(res-1)*3))) 473 model[candidateModel].vertexNormals[i*3..i*3+3] = [0, -1, 0]; 474 475 startingIdx += (res-1)*3; 476 foreach(i; parallel(iota(startingIdx,startingIdx+(res-1)*3))) 477 model[candidateModel].vertexNormals[i*3..i*3+3] = [0, 1, 0]; 478 } 479 480 swGen.stop(); 481 482 if (!continueProcessing()) return; 483 484 gs_actual = current; 485 486 auto oldModel = currentModel; 487 currentModel = candidateModel; 488 489 Mallocator.instance.dispose(model[oldModel].vertex); 490 Mallocator.instance.dispose(model[oldModel].vertexNormals); 491 492 model[oldModel].vertex = null; 493 model[oldModel].vertexNormals = null; 494 495 writeln("Vase regenerated in: ", swGen.peek.total!"msecs", "ms"); 496 }