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 }