be43154fd93e4791d1b72d018facfd3bca5e2706
[blender-staging.git] / intern / cycles / render / mesh_volume.cpp
1 /*
2  * Copyright 2011-2016 Blender Foundation
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16
17 #include "render/mesh.h"
18 #include "render/attribute.h"
19 #include "render/scene.h"
20
21 #include "util/util_foreach.h"
22 #include "util/util_logging.h"
23 #include "util/util_progress.h"
24 #include "util/util_types.h"
25
26 CCL_NAMESPACE_BEGIN
27
28 static size_t compute_voxel_index(const int3 &resolution, size_t x, size_t y, size_t z)
29 {
30         if(x == -1 || x >= resolution.x) {
31                 return -1;
32         }
33
34         if(y == -1 || y >= resolution.y) {
35                 return -1;
36         }
37
38         if(z == -1 || z >= resolution.z) {
39                 return -1;
40         }
41
42         return x + y*resolution.x + z*resolution.x*resolution.y;
43 }
44
45 struct QuadData {
46         int v0, v1, v2, v3;
47
48         float3 normal;
49 };
50
51 enum {
52         QUAD_X_MIN = 0,
53         QUAD_X_MAX = 1,
54         QUAD_Y_MIN = 2,
55         QUAD_Y_MAX = 3,
56         QUAD_Z_MIN = 4,
57         QUAD_Z_MAX = 5,
58 };
59
60 const int quads_indices[6][4] = {
61         /* QUAD_X_MIN */
62         { 4, 0, 3, 7 },
63         /* QUAD_X_MAX */
64         { 1, 5, 6, 2 },
65         /* QUAD_Y_MIN */
66         { 4, 5, 1, 0 },
67         /* QUAD_Y_MAX */
68         { 3, 2, 6, 7 },
69         /* QUAD_Z_MIN */
70         { 0, 1, 2, 3 },
71         /* QUAD_Z_MAX */
72         { 5, 4, 7, 6 },
73 };
74
75 const float3 quads_normals[6] = {
76         /* QUAD_X_MIN */
77         make_float3(-1.0f, 0.0f, 0.0f),
78         /* QUAD_X_MAX */
79         make_float3(1.0f, 0.0f, 0.0f),
80         /* QUAD_Y_MIN */
81         make_float3(0.0f, -1.0f, 0.0f),
82         /* QUAD_Y_MAX */
83         make_float3(0.0f, 1.0f, 0.0f),
84         /* QUAD_Z_MIN */
85         make_float3(0.0f, 0.0f, -1.0f),
86         /* QUAD_Z_MAX */
87         make_float3(0.0f, 0.0f, 1.0f),
88 };
89
90 static void create_quad(int3 corners[8], vector<int3> &vertices, vector<QuadData> &quads, int face_index)
91 {
92         size_t vertex_offset = vertices.size();
93
94         QuadData quad;
95         quad.v0 = vertex_offset + 0;
96         quad.v1 = vertex_offset + 1;
97         quad.v2 = vertex_offset + 2;
98         quad.v3 = vertex_offset + 3;
99         quad.normal = quads_normals[face_index];
100
101         quads.push_back(quad);
102
103         vertices.push_back(corners[quads_indices[face_index][0]]);
104         vertices.push_back(corners[quads_indices[face_index][1]]);
105         vertices.push_back(corners[quads_indices[face_index][2]]);
106         vertices.push_back(corners[quads_indices[face_index][3]]);
107 }
108
109 struct VolumeParams {
110         int3 resolution;
111         float3 cell_size;
112         float3 start_point;
113         int pad_size;
114 };
115
116 static const int CUBE_SIZE = 8;
117
118 /* Create a mesh from a volume.
119  *
120  * The way the algorithm works is as follows:
121  *
122  * - the coordinates of active voxels from a dense volume (or 3d image) are
123  * gathered inside an auxialliary volume.
124  * - each set of coordinates of an CUBE_SIZE cube are mapped to the same
125  * coordinate of the auxilliary volume.
126  * - quads are created between active and non-active voxels in the auxialliary
127  * volume to generate a tight mesh around the volume.
128  */
129 class VolumeMeshBuilder {
130         /* Auxilliary volume that is used to check if a node already added. */
131         vector<char> grid;
132
133         /* The resolution of the auxilliary volume, set to be equal to 1/CUBE_SIZE
134          * of the original volume on each axis. */
135         int3 res;
136
137         size_t number_of_nodes;
138
139         /* Offset due to padding in the original grid. Padding will transform the
140          * coordinates of the original grid from 0...res to -padding...res+padding,
141          * so some coordinates are negative, and we need to properly account for
142          * them. */
143         int3 pad_offset;
144
145         VolumeParams *params;
146
147 public:
148         VolumeMeshBuilder(VolumeParams *volume_params);
149
150         void add_node(int x, int y, int z);
151
152         void add_node_with_padding(int x, int y, int z);
153
154         void create_mesh(vector<float3> &vertices,
155                                          vector<int> &indices,
156                                          vector<float3> &face_normals);
157
158 private:
159         void generate_vertices_and_quads(vector<int3> &vertices_is,
160                                                                          vector<QuadData> &quads);
161
162         void deduplicate_vertices(vector<int3> &vertices,
163                                                           vector<QuadData> &quads);
164
165         void convert_object_space(const vector<int3> &vertices,
166                                                           vector<float3> &out_vertices);
167
168         void convert_quads_to_tris(const vector<QuadData> &quads,
169                                                            vector<int> &tris,
170                                                            vector<float3> &face_normals);
171 };
172
173 VolumeMeshBuilder::VolumeMeshBuilder(VolumeParams *volume_params)
174 {
175         params = volume_params;
176         number_of_nodes = 0;
177
178         const size_t x = divide_up(params->resolution.x, CUBE_SIZE);
179         const size_t y = divide_up(params->resolution.y, CUBE_SIZE);
180         const size_t z = divide_up(params->resolution.z, CUBE_SIZE);
181
182         /* Adding 2*pad_size since we pad in both positive and negative directions
183          * along the axis. */
184         const size_t px = divide_up(params->resolution.x + 2*params->pad_size, CUBE_SIZE);
185         const size_t py = divide_up(params->resolution.y + 2*params->pad_size, CUBE_SIZE);
186         const size_t pz = divide_up(params->resolution.z + 2*params->pad_size, CUBE_SIZE);
187
188         res = make_int3(px, py, pz);
189         pad_offset = make_int3(px - x, py - y, pz - z);
190
191         grid.resize(px*py*pz, 0);
192 }
193
194 void VolumeMeshBuilder::add_node(int x, int y, int z)
195 {
196         /* Map coordinates to index space. */
197         const int index_x = (x/CUBE_SIZE) + pad_offset.x;
198         const int index_y = (y/CUBE_SIZE) + pad_offset.y;
199         const int index_z = (z/CUBE_SIZE) + pad_offset.z;
200
201         assert((index_x >= 0) && (index_y >= 0) && (index_z >= 0));
202
203         const size_t index = compute_voxel_index(res, index_x, index_y, index_z);
204
205         /* We already have a node here. */
206         if(grid[index] == 1) {
207                 return;
208         }
209
210         ++number_of_nodes;
211
212         grid[index] = 1;
213 }
214
215 void VolumeMeshBuilder::add_node_with_padding(int x, int y, int z)
216 {
217         for(int px = x - params->pad_size; px < x + params->pad_size; ++px) {
218                 for(int py = y - params->pad_size; py < y + params->pad_size; ++py) {
219                         for(int pz = z - params->pad_size; pz < z + params->pad_size; ++pz) {
220                                 add_node(px, py, pz);
221                         }
222                 }
223         }
224 }
225
226 void VolumeMeshBuilder::create_mesh(vector<float3> &vertices,
227                                                                         vector<int> &indices,
228                                                                         vector<float3> &face_normals)
229 {
230         /* We create vertices in index space (is), and only convert them to object
231          * space when done. */
232         vector<int3> vertices_is;
233         vector<QuadData> quads;
234
235         generate_vertices_and_quads(vertices_is, quads);
236
237         deduplicate_vertices(vertices_is, quads);
238
239         convert_object_space(vertices_is, vertices);
240
241         convert_quads_to_tris(quads, indices, face_normals);
242 }
243
244 void VolumeMeshBuilder::generate_vertices_and_quads(
245                 vector<ccl::int3> &vertices_is,
246                 vector<QuadData> &quads)
247 {
248         /* Overallocation, we could count the number of quads and vertices to create
249          * in a pre-pass if memory becomes an issue. */
250         vertices_is.reserve(number_of_nodes*8);
251         quads.reserve(number_of_nodes*6);
252
253         for(int z = 0; z < res.z; ++z) {
254                 for(int y = 0; y < res.y; ++y) {
255                         for(int x = 0; x < res.x; ++x) {
256                                 size_t voxel_index = compute_voxel_index(res, x, y, z);
257                                 if(grid[voxel_index] == 0) {
258                                         continue;
259                                 }
260
261                                 /* Compute min and max coords of the node in index space. */
262                                 int3 min = make_int3((x - pad_offset.x)*CUBE_SIZE,
263                                                                          (y - pad_offset.y)*CUBE_SIZE,
264                                                                          (z - pad_offset.z)*CUBE_SIZE);
265
266                                 /* Maximum is just CUBE_SIZE voxels away from minimum on each axis. */
267                                 int3 max = make_int3(min.x + CUBE_SIZE, min.y + CUBE_SIZE, min.z + CUBE_SIZE);
268
269                                 int3 corners[8] = {
270                                         make_int3(min[0], min[1], min[2]),
271                                         make_int3(max[0], min[1], min[2]),
272                                         make_int3(max[0], max[1], min[2]),
273                                         make_int3(min[0], max[1], min[2]),
274                                         make_int3(min[0], min[1], max[2]),
275                                         make_int3(max[0], min[1], max[2]),
276                                         make_int3(max[0], max[1], max[2]),
277                                         make_int3(min[0], max[1], max[2]),
278                                 };
279
280                                 /* Only create a quad if on the border between an active and
281                                  * an inactive node.
282                                  */
283
284                                 voxel_index = compute_voxel_index(res, x - 1, y, z);
285                                 if(voxel_index == -1 || grid[voxel_index] == 0) {
286                                         create_quad(corners, vertices_is, quads, QUAD_X_MIN);
287                                 }
288
289                                 voxel_index = compute_voxel_index(res, x + 1, y, z);
290                                 if(voxel_index == -1 || grid[voxel_index] == 0) {
291                                         create_quad(corners, vertices_is, quads, QUAD_X_MAX);
292                                 }
293
294                                 voxel_index = compute_voxel_index(res, x, y - 1, z);
295                                 if(voxel_index == -1 || grid[voxel_index] == 0) {
296                                         create_quad(corners, vertices_is, quads, QUAD_Y_MIN);
297                                 }
298
299                                 voxel_index = compute_voxel_index(res, x, y + 1, z);
300                                 if(voxel_index == -1 || grid[voxel_index] == 0) {
301                                         create_quad(corners, vertices_is, quads, QUAD_Y_MAX);
302                                 }
303
304                                 voxel_index = compute_voxel_index(res, x, y, z - 1);
305                                 if(voxel_index == -1 || grid[voxel_index] == 0) {
306                                         create_quad(corners, vertices_is, quads, QUAD_Z_MIN);
307                                 }
308
309                                 voxel_index = compute_voxel_index(res, x, y, z + 1);
310                                 if(voxel_index == -1 || grid[voxel_index] == 0) {
311                                         create_quad(corners, vertices_is, quads, QUAD_Z_MAX);
312                                 }
313                         }
314                 }
315         }
316 }
317
318 void VolumeMeshBuilder::deduplicate_vertices(vector<int3> &vertices,
319                                                                                          vector<QuadData> &quads)
320 {
321         vector<int3> sorted_vertices = vertices;
322         std::sort(sorted_vertices.begin(), sorted_vertices.end());
323         vector<int3>::iterator it = std::unique(sorted_vertices.begin(), sorted_vertices.end());
324         sorted_vertices.resize(std::distance(sorted_vertices.begin(), it));
325
326         vector<QuadData> new_quads = quads;
327
328         for(size_t i = 0; i < vertices.size(); ++i) {
329                 for(size_t j = 0; j < sorted_vertices.size(); ++j) {
330                         if(vertices[i] != sorted_vertices[j]) {
331                                 continue;
332                         }
333
334                         for(int k = 0; k < quads.size(); ++k) {
335                                 if(quads[k].v0 == i) {
336                                         new_quads[k].v0 = j;
337                                 }
338                                 else if(quads[k].v1 == i) {
339                                         new_quads[k].v1 = j;
340                                 }
341                                 else if(quads[k].v2 == i) {
342                                         new_quads[k].v2 = j;
343                                 }
344                                 else if(quads[k].v3 == i) {
345                                         new_quads[k].v3 = j;
346                                 }
347                         }
348
349                         break;
350                 }
351         }
352
353         vertices = sorted_vertices;
354         quads = new_quads;
355 }
356
357 void VolumeMeshBuilder::convert_object_space(const vector<int3> &vertices,
358                                                                                          vector<float3> &out_vertices)
359 {
360         out_vertices.reserve(vertices.size());
361
362         for(size_t i = 0; i < vertices.size(); ++i) {
363                 float3 vertex = make_float3(vertices[i].x, vertices[i].y, vertices[i].z);
364                 vertex *= params->cell_size;
365                 vertex += params->start_point;
366
367                 out_vertices.push_back(vertex);
368         }
369 }
370
371 void VolumeMeshBuilder::convert_quads_to_tris(const vector<QuadData> &quads,
372                                                                                           vector<int> &tris,
373                                                                                           vector<float3> &face_normals)
374 {
375         int index_offset = 0;
376         tris.resize(quads.size()*6);
377         face_normals.reserve(quads.size()*2);
378
379         for(size_t i = 0; i < quads.size(); ++i) {
380                 tris[index_offset++] = quads[i].v0;
381                 tris[index_offset++] = quads[i].v2;
382                 tris[index_offset++] = quads[i].v1;
383
384                 face_normals.push_back(quads[i].normal);
385
386                 tris[index_offset++] = quads[i].v0;
387                 tris[index_offset++] = quads[i].v3;
388                 tris[index_offset++] = quads[i].v2;
389
390                 face_normals.push_back(quads[i].normal);
391         }
392 }
393
394 /* ************************************************************************** */
395
396 /* For debugging: render the created mesh using the default diffuse shader. */
397 //#define RENDER_DIFFUSE
398
399 struct VoxelAttributeGrid {
400         float *data;
401         int channels;
402 };
403
404 void MeshManager::create_volume_mesh(Scene *scene,
405                                                                          Mesh *mesh,
406                                                                          Progress& progress)
407 {
408         string msg = string_printf("Computing Volume Mesh %s", mesh->name.c_str());
409         progress.set_status("Updating Mesh", msg);
410
411         vector<VoxelAttributeGrid> voxel_grids;
412
413         /* Compute volume parameters. */
414         VolumeParams volume_params;
415         volume_params.resolution = make_int3(0, 0, 0);
416
417         foreach(Attribute& attr, mesh->attributes.attributes) {
418                 if(attr.element != ATTR_ELEMENT_VOXEL) {
419                         continue;
420                 }
421
422                 VoxelAttribute *voxel = attr.data_voxel();
423                 device_memory *image_memory = scene->image_manager->image_memory(voxel->slot);
424                 int3 resolution = make_int3(image_memory->data_width,
425                                     image_memory->data_height,
426                                             image_memory->data_depth);
427
428                 if(volume_params.resolution == make_int3(0, 0, 0)) {
429                         volume_params.resolution = resolution;
430                 }
431                 else if(volume_params.resolution != resolution) {
432                         VLOG(1) << "Can't create volume mesh, all voxel grid resolutions must be equal\n";
433                         return;
434                 }
435
436                 VoxelAttributeGrid voxel_grid;
437                 voxel_grid.data = static_cast<float*>(image_memory->host_pointer);
438                 voxel_grid.channels = image_memory->data_elements;
439                 voxel_grids.push_back(voxel_grid);
440         }
441
442         if(voxel_grids.empty()) {
443                 return;
444         }
445
446         int pad_size = 0;
447
448         foreach(Shader *shader, mesh->used_shaders) {
449                 if(!shader->has_volume) {
450                         continue;
451                 }
452
453                 if(shader->volume_interpolation_method == VOLUME_INTERPOLATION_LINEAR) {
454                         pad_size = max(1, pad_size);
455                 }
456                 else if(shader->volume_interpolation_method == VOLUME_INTERPOLATION_CUBIC) {
457                         pad_size = max(2, pad_size);
458                 }
459         }
460
461         /* Compute start point and cell size from transform. */
462         Attribute *attr = mesh->attributes.find(ATTR_STD_GENERATED_TRANSFORM);
463         const int3 resolution = volume_params.resolution;
464         float3 start_point = make_float3(0.0f, 0.0f, 0.0f);
465         float3 cell_size = make_float3(1.0f/resolution.x,
466                                                                    1.0f/resolution.y,
467                                                                    1.0f/resolution.z);
468
469         if(attr) {
470                 const Transform *tfm = attr->data_transform();
471                 const Transform itfm = transform_inverse(*tfm);
472                 start_point = transform_point(&itfm, start_point);
473                 cell_size = transform_direction(&itfm, cell_size);
474         }
475
476         volume_params.start_point = start_point;
477         volume_params.cell_size = cell_size;
478         volume_params.pad_size = pad_size;
479
480         VolumeMeshBuilder builder(&volume_params);
481         const float isovalue = mesh->volume_isovalue;
482
483         for(int z = 0; z < resolution.z; ++z) {
484                 for(int y = 0; y < resolution.y; ++y) {
485                         for(int x = 0; x < resolution.x; ++x) {
486                                 size_t voxel_index = compute_voxel_index(resolution, x, y, z);
487
488                                 for(size_t i = 0; i < voxel_grids.size(); ++i) {
489                                         const VoxelAttributeGrid &voxel_grid = voxel_grids[i];
490
491                                         if(voxel_grid.channels == 1) {
492                                                 if(voxel_grid.data[voxel_index] >= isovalue) {
493                                                         builder.add_node_with_padding(x, y, z);
494                                                         break;
495                                                 }
496                                         }
497                                         else if(voxel_grid.channels == 3) {
498                                                 voxel_index = compute_voxel_index(resolution, x*3, y, z);
499
500                                                 if(voxel_grid.data[voxel_index] >= isovalue) {
501                                                         builder.add_node_with_padding(x, y, z);
502                                                         break;
503                                                 }
504
505                                                 if(voxel_grid.data[voxel_index + 1] >= isovalue) {
506                                                         builder.add_node_with_padding(x, y, z);
507                                                         break;
508                                                 }
509
510                                                 if(voxel_grid.data[voxel_index + 2] >= isovalue) {
511                                                         builder.add_node_with_padding(x, y, z);
512                                                         break;
513                                                 }
514                                         }
515                                         else if(voxel_grid.channels == 4) {
516                                                 voxel_index = compute_voxel_index(resolution, x*4, y, z);
517
518                                                 /* check alpha first */
519                                                 if(voxel_grid.data[voxel_index + 3] < isovalue) {
520                                                         continue;
521                                                 }
522
523                                                 if(voxel_grid.data[voxel_index] >= isovalue) {
524                                                         builder.add_node_with_padding(x, y, z);
525                                                         continue;
526                                                 }
527
528                                                 if(voxel_grid.data[voxel_index + 1] >= isovalue) {
529                                                         builder.add_node_with_padding(x, y, z);
530                                                         continue;
531                                                 }
532
533                                                 if(voxel_grid.data[voxel_index + 2] >= isovalue) {
534                                                         builder.add_node_with_padding(x, y, z);
535                                                         continue;
536                                                 }
537                                         }
538                                 }
539                         }
540                 }
541         }
542
543         vector<float3> vertices;
544         vector<int> indices;
545         vector<float3> face_normals;
546         builder.create_mesh(vertices, indices, face_normals);
547
548 #ifdef RENDER_DIFFUSE
549         int shader = mesh->used_shaders[0]->id;
550 #else
551         int shader = mesh->shader[0];
552 #endif
553
554         mesh->clear(true);
555         mesh->reserve_mesh(vertices.size(), indices.size()/3);
556
557         for(size_t i = 0; i < vertices.size(); ++i) {
558                 mesh->add_vertex(vertices[i]);
559         }
560
561         for(size_t i = 0; i < indices.size(); i += 3) {
562                 mesh->add_triangle(indices[i], indices[i + 1], indices[i + 2], shader, false);
563         }
564
565         Attribute *attr_fN = mesh->attributes.add(ATTR_STD_FACE_NORMAL);
566         float3 *fN = attr_fN->data_float3();
567
568         for(size_t i = 0; i < face_normals.size(); ++i) {
569                 fN[i] = face_normals[i];
570         }
571
572         VLOG(1) << "Memory usage volume mesh: "
573                         << ((vertices.size() + face_normals.size())*sizeof(float3) + indices.size()*sizeof(int))/(1024.0*1024.0)
574                         << "Mb.";
575
576         VLOG(1) << "Memory usage volume grid: "
577                         << (resolution.x*resolution.y*resolution.z*sizeof(float))/(1024.0*1024.0)
578                         << "Mb.";
579 }
580
581 CCL_NAMESPACE_END