mirror of
https://github.com/schrodinger/pymol-open-source.git
synced 2026-06-04 20:04:21 +08:00
738 lines
23 KiB
C++
738 lines
23 KiB
C++
/**
|
|
* @file Marching cubes
|
|
*
|
|
* Based on https://github.com/ilastik/marching_cubes
|
|
*
|
|
* License: BSD-3-Clause
|
|
*
|
|
* Copyright 2018, The ilastik development team
|
|
* Copyright 2020, Schrodinger, Inc.
|
|
*/
|
|
|
|
#include "marching_cubes.h"
|
|
|
|
#include <algorithm>
|
|
#include <cassert>
|
|
#include <cmath>
|
|
#include <unordered_map>
|
|
#include <vector>
|
|
|
|
#ifdef PYMOL_OPENMP
|
|
#include <omp.h>
|
|
#endif
|
|
|
|
namespace mc
|
|
{
|
|
|
|
static const unsigned short EDGE_TABLE[256] = {0x0, 0x109, 0x203, 0x30a, 0x406,
|
|
0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
|
|
0x190, 0x99, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, 0xb9f,
|
|
0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, 0x230, 0x339, 0x33, 0x13a, 0x636, 0x73f,
|
|
0x435, 0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, 0x3a0,
|
|
0x2a9, 0x1a3, 0xaa, 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6,
|
|
0xfaa, 0xea3, 0xda9, 0xca0, 0x460, 0x569, 0x663, 0x76a, 0x66, 0x16f, 0x265,
|
|
0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60, 0x5f0, 0x4f9,
|
|
0x7f3, 0x6fa, 0x1f6, 0xff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa,
|
|
0x8f3, 0xbf9, 0xaf0, 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55, 0x15c,
|
|
0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, 0x7c0, 0x6c9, 0x5c3,
|
|
0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3,
|
|
0x9c9, 0x8c0, 0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0xcc,
|
|
0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0, 0x950, 0x859, 0xb53, 0xa5a,
|
|
0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x55, 0x35f, 0x256, 0x55a, 0x453, 0x759,
|
|
0x650, 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, 0x3f5,
|
|
0xff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, 0xb60, 0xa69, 0x963, 0x86a, 0xf66,
|
|
0xe6f, 0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x66, 0x76a, 0x663, 0x569, 0x460,
|
|
0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af,
|
|
0x7a6, 0xaa, 0x1a3, 0x2a9, 0x3a0, 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f,
|
|
0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33, 0x339, 0x230, 0xe90,
|
|
0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795, 0x49f, 0x596,
|
|
0x29a, 0x393, 0x99, 0x190, 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905,
|
|
0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0};
|
|
|
|
static const signed char TRIANGLE_TABLE[256][16] = {
|
|
{-1},
|
|
{0, 8, 3, -1},
|
|
{0, 1, 9, -1},
|
|
{1, 8, 3, 9, 8, 1, -1},
|
|
{1, 2, 10, -1},
|
|
{0, 8, 3, 1, 2, 10, -1},
|
|
{9, 2, 10, 0, 2, 9, -1},
|
|
{2, 8, 3, 2, 10, 8, 10, 9, 8, -1},
|
|
{3, 11, 2, -1},
|
|
{0, 11, 2, 8, 11, 0, -1},
|
|
{1, 9, 0, 2, 3, 11, -1},
|
|
{1, 11, 2, 1, 9, 11, 9, 8, 11, -1},
|
|
{3, 10, 1, 11, 10, 3, -1},
|
|
{0, 10, 1, 0, 8, 10, 8, 11, 10, -1},
|
|
{3, 9, 0, 3, 11, 9, 11, 10, 9, -1},
|
|
{9, 8, 10, 10, 8, 11, -1},
|
|
{4, 7, 8, -1},
|
|
{4, 3, 0, 7, 3, 4, -1},
|
|
{0, 1, 9, 8, 4, 7, -1},
|
|
{4, 1, 9, 4, 7, 1, 7, 3, 1, -1},
|
|
{1, 2, 10, 8, 4, 7, -1},
|
|
{3, 4, 7, 3, 0, 4, 1, 2, 10, -1},
|
|
{9, 2, 10, 9, 0, 2, 8, 4, 7, -1},
|
|
{2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1},
|
|
{8, 4, 7, 3, 11, 2, -1},
|
|
{11, 4, 7, 11, 2, 4, 2, 0, 4, -1},
|
|
{9, 0, 1, 8, 4, 7, 2, 3, 11, -1},
|
|
{4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1},
|
|
{3, 10, 1, 3, 11, 10, 7, 8, 4, -1},
|
|
{1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1},
|
|
{4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1},
|
|
{4, 7, 11, 4, 11, 9, 9, 11, 10, -1},
|
|
{9, 5, 4, -1},
|
|
{9, 5, 4, 0, 8, 3, -1},
|
|
{0, 5, 4, 1, 5, 0, -1},
|
|
{8, 5, 4, 8, 3, 5, 3, 1, 5, -1},
|
|
{1, 2, 10, 9, 5, 4, -1},
|
|
{3, 0, 8, 1, 2, 10, 4, 9, 5, -1},
|
|
{5, 2, 10, 5, 4, 2, 4, 0, 2, -1},
|
|
{2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1},
|
|
{9, 5, 4, 2, 3, 11, -1},
|
|
{0, 11, 2, 0, 8, 11, 4, 9, 5, -1},
|
|
{0, 5, 4, 0, 1, 5, 2, 3, 11, -1},
|
|
{2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1},
|
|
{10, 3, 11, 10, 1, 3, 9, 5, 4, -1},
|
|
{4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1},
|
|
{5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1},
|
|
{5, 4, 8, 5, 8, 10, 10, 8, 11, -1},
|
|
{9, 7, 8, 5, 7, 9, -1},
|
|
{9, 3, 0, 9, 5, 3, 5, 7, 3, -1},
|
|
{0, 7, 8, 0, 1, 7, 1, 5, 7, -1},
|
|
{1, 5, 3, 3, 5, 7, -1},
|
|
{9, 7, 8, 9, 5, 7, 10, 1, 2, -1},
|
|
{10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1},
|
|
{8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1},
|
|
{2, 10, 5, 2, 5, 3, 3, 5, 7, -1},
|
|
{7, 9, 5, 7, 8, 9, 3, 11, 2, -1},
|
|
{9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1},
|
|
{2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1},
|
|
{11, 2, 1, 11, 1, 7, 7, 1, 5, -1},
|
|
{9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1},
|
|
{5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
|
|
{11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
|
|
{11, 10, 5, 7, 11, 5, -1},
|
|
{10, 6, 5, -1},
|
|
{0, 8, 3, 5, 10, 6, -1},
|
|
{9, 0, 1, 5, 10, 6, -1},
|
|
{1, 8, 3, 1, 9, 8, 5, 10, 6, -1},
|
|
{1, 6, 5, 2, 6, 1, -1},
|
|
{1, 6, 5, 1, 2, 6, 3, 0, 8, -1},
|
|
{9, 6, 5, 9, 0, 6, 0, 2, 6, -1},
|
|
{5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1},
|
|
{2, 3, 11, 10, 6, 5, -1},
|
|
{11, 0, 8, 11, 2, 0, 10, 6, 5, -1},
|
|
{0, 1, 9, 2, 3, 11, 5, 10, 6, -1},
|
|
{5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1},
|
|
{6, 3, 11, 6, 5, 3, 5, 1, 3, -1},
|
|
{0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1},
|
|
{3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1},
|
|
{6, 5, 9, 6, 9, 11, 11, 9, 8, -1},
|
|
{5, 10, 6, 4, 7, 8, -1},
|
|
{4, 3, 0, 4, 7, 3, 6, 5, 10, -1},
|
|
{1, 9, 0, 5, 10, 6, 8, 4, 7, -1},
|
|
{10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1},
|
|
{6, 1, 2, 6, 5, 1, 4, 7, 8, -1},
|
|
{1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1},
|
|
{8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1},
|
|
{7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
|
|
{3, 11, 2, 7, 8, 4, 10, 6, 5, -1},
|
|
{5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1},
|
|
{0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1},
|
|
{9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
|
|
{8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1},
|
|
{5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
|
|
{0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
|
|
{6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1},
|
|
{10, 4, 9, 6, 4, 10, -1},
|
|
{4, 10, 6, 4, 9, 10, 0, 8, 3, -1},
|
|
{10, 0, 1, 10, 6, 0, 6, 4, 0, -1},
|
|
{8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1},
|
|
{1, 4, 9, 1, 2, 4, 2, 6, 4, -1},
|
|
{3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1},
|
|
{0, 2, 4, 4, 2, 6, -1},
|
|
{8, 3, 2, 8, 2, 4, 4, 2, 6, -1},
|
|
{10, 4, 9, 10, 6, 4, 11, 2, 3, -1},
|
|
{0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1},
|
|
{3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1},
|
|
{6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
|
|
{9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1},
|
|
{8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
|
|
{3, 11, 6, 3, 6, 0, 0, 6, 4, -1},
|
|
{6, 4, 8, 11, 6, 8, -1},
|
|
{7, 10, 6, 7, 8, 10, 8, 9, 10, -1},
|
|
{0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1},
|
|
{10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1},
|
|
{10, 6, 7, 10, 7, 1, 1, 7, 3, -1},
|
|
{1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1},
|
|
{2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
|
|
{7, 8, 0, 7, 0, 6, 6, 0, 2, -1},
|
|
{7, 3, 2, 6, 7, 2, -1},
|
|
{2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1},
|
|
{2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
|
|
{1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
|
|
{11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1},
|
|
{8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
|
|
{0, 9, 1, 11, 6, 7, -1},
|
|
{7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1},
|
|
{7, 11, 6, -1},
|
|
{7, 6, 11, -1},
|
|
{3, 0, 8, 11, 7, 6, -1},
|
|
{0, 1, 9, 11, 7, 6, -1},
|
|
{8, 1, 9, 8, 3, 1, 11, 7, 6, -1},
|
|
{10, 1, 2, 6, 11, 7, -1},
|
|
{1, 2, 10, 3, 0, 8, 6, 11, 7, -1},
|
|
{2, 9, 0, 2, 10, 9, 6, 11, 7, -1},
|
|
{6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1},
|
|
{7, 2, 3, 6, 2, 7, -1},
|
|
{7, 0, 8, 7, 6, 0, 6, 2, 0, -1},
|
|
{2, 7, 6, 2, 3, 7, 0, 1, 9, -1},
|
|
{1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1},
|
|
{10, 7, 6, 10, 1, 7, 1, 3, 7, -1},
|
|
{10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1},
|
|
{0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1},
|
|
{7, 6, 10, 7, 10, 8, 8, 10, 9, -1},
|
|
{6, 8, 4, 11, 8, 6, -1},
|
|
{3, 6, 11, 3, 0, 6, 0, 4, 6, -1},
|
|
{8, 6, 11, 8, 4, 6, 9, 0, 1, -1},
|
|
{9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1},
|
|
{6, 8, 4, 6, 11, 8, 2, 10, 1, -1},
|
|
{1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1},
|
|
{4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1},
|
|
{10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
|
|
{8, 2, 3, 8, 4, 2, 4, 6, 2, -1},
|
|
{0, 4, 2, 4, 6, 2, -1},
|
|
{1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1},
|
|
{1, 9, 4, 1, 4, 2, 2, 4, 6, -1},
|
|
{8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1},
|
|
{10, 1, 0, 10, 0, 6, 6, 0, 4, -1},
|
|
{4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
|
|
{10, 9, 4, 6, 10, 4, -1},
|
|
{4, 9, 5, 7, 6, 11, -1},
|
|
{0, 8, 3, 4, 9, 5, 11, 7, 6, -1},
|
|
{5, 0, 1, 5, 4, 0, 7, 6, 11, -1},
|
|
{11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1},
|
|
{9, 5, 4, 10, 1, 2, 7, 6, 11, -1},
|
|
{6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1},
|
|
{7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1},
|
|
{3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
|
|
{7, 2, 3, 7, 6, 2, 5, 4, 9, -1},
|
|
{9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1},
|
|
{3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1},
|
|
{6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
|
|
{9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1},
|
|
{1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
|
|
{4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
|
|
{7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1},
|
|
{6, 9, 5, 6, 11, 9, 11, 8, 9, -1},
|
|
{3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1},
|
|
{0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1},
|
|
{6, 11, 3, 6, 3, 5, 5, 3, 1, -1},
|
|
{1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1},
|
|
{0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
|
|
{11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
|
|
{6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1},
|
|
{5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1},
|
|
{9, 5, 6, 9, 6, 0, 0, 6, 2, -1},
|
|
{1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
|
|
{1, 5, 6, 2, 1, 6, -1},
|
|
{1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
|
|
{10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1},
|
|
{0, 3, 8, 5, 6, 10, -1},
|
|
{10, 5, 6, -1},
|
|
{11, 5, 10, 7, 5, 11, -1},
|
|
{11, 5, 10, 11, 7, 5, 8, 3, 0, -1},
|
|
{5, 11, 7, 5, 10, 11, 1, 9, 0, -1},
|
|
{10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1},
|
|
{11, 1, 2, 11, 7, 1, 7, 5, 1, -1},
|
|
{0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1},
|
|
{9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1},
|
|
{7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
|
|
{2, 5, 10, 2, 3, 5, 3, 7, 5, -1},
|
|
{8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1},
|
|
{9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1},
|
|
{9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
|
|
{1, 3, 5, 3, 7, 5, -1},
|
|
{0, 8, 7, 0, 7, 1, 1, 7, 5, -1},
|
|
{9, 0, 3, 9, 3, 5, 5, 3, 7, -1},
|
|
{9, 8, 7, 5, 9, 7, -1},
|
|
{5, 8, 4, 5, 10, 8, 10, 11, 8, -1},
|
|
{5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1},
|
|
{0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1},
|
|
{10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
|
|
{2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1},
|
|
{0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
|
|
{0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
|
|
{9, 4, 5, 2, 11, 3, -1},
|
|
{2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1},
|
|
{5, 10, 2, 5, 2, 4, 4, 2, 0, -1},
|
|
{3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
|
|
{5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1},
|
|
{8, 4, 5, 8, 5, 3, 3, 5, 1, -1},
|
|
{0, 4, 5, 1, 0, 5, -1},
|
|
{8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1},
|
|
{9, 4, 5, -1},
|
|
{4, 11, 7, 4, 9, 11, 9, 10, 11, -1},
|
|
{0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1},
|
|
{1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1},
|
|
{3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
|
|
{4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1},
|
|
{9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
|
|
{11, 7, 4, 11, 4, 2, 2, 4, 0, -1},
|
|
{11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1},
|
|
{2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1},
|
|
{9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
|
|
{3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
|
|
{1, 10, 2, 8, 7, 4, -1},
|
|
{4, 9, 1, 4, 1, 7, 7, 1, 3, -1},
|
|
{4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1},
|
|
{4, 0, 3, 7, 4, 3, -1},
|
|
{4, 8, 7, -1},
|
|
{9, 10, 8, 10, 11, 8, -1},
|
|
{3, 0, 9, 3, 9, 11, 11, 9, 10, -1},
|
|
{0, 1, 10, 0, 10, 8, 8, 10, 11, -1},
|
|
{3, 1, 10, 11, 3, 10, -1},
|
|
{1, 2, 11, 1, 11, 9, 9, 11, 8, -1},
|
|
{3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1},
|
|
{0, 2, 11, 8, 0, 11, -1},
|
|
{3, 2, 11, -1},
|
|
{2, 3, 8, 2, 8, 10, 10, 8, 9, -1},
|
|
{9, 10, 2, 0, 9, 2, -1},
|
|
{2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1},
|
|
{1, 10, 2, -1},
|
|
{1, 3, 8, 9, 1, 8, -1},
|
|
{0, 9, 1, -1},
|
|
{0, 3, 8, -1},
|
|
{-1},
|
|
};
|
|
|
|
constexpr unsigned EDGES_PER_CELL = 3;
|
|
|
|
static size_t vertexId(size_t x, size_t y, size_t z, size_t xDim, size_t yDim)
|
|
{
|
|
return EDGES_PER_CELL * (z * yDim * xDim + y * xDim + x);
|
|
}
|
|
|
|
inline size_t edgeId2z(size_t eid, size_t xDim, size_t yDim)
|
|
{
|
|
return (eid / EDGES_PER_CELL) / (yDim * xDim);
|
|
}
|
|
|
|
static size_t edgeId(
|
|
size_t x, size_t y, size_t z, size_t edgeNumber, size_t xDim, size_t yDim)
|
|
{
|
|
switch (edgeNumber) {
|
|
case 0:
|
|
return vertexId(x, y, z, xDim, yDim) + 1;
|
|
case 1:
|
|
return vertexId(x, y + 1, z, xDim, yDim);
|
|
case 2:
|
|
return vertexId(x + 1, y, z, xDim, yDim) + 1;
|
|
case 3:
|
|
return vertexId(x, y, z, xDim, yDim);
|
|
case 4:
|
|
return vertexId(x, y, z + 1, xDim, yDim) + 1;
|
|
case 5:
|
|
return vertexId(x, y + 1, z + 1, xDim, yDim);
|
|
case 6:
|
|
return vertexId(x + 1, y, z + 1, xDim, yDim) + 1;
|
|
case 7:
|
|
return vertexId(x, y, z + 1, xDim, yDim);
|
|
case 8:
|
|
return vertexId(x, y, z, xDim, yDim) + 2;
|
|
case 9:
|
|
return vertexId(x, y + 1, z, xDim, yDim) + 2;
|
|
case 10:
|
|
return vertexId(x + 1, y + 1, z, xDim, yDim) + 2;
|
|
case 11:
|
|
return vertexId(x + 1, y, z, xDim, yDim) + 2;
|
|
default:
|
|
assert(false);
|
|
return -1;
|
|
}
|
|
}
|
|
|
|
inline float mix(float v0, float v1, float a)
|
|
{
|
|
return v0 * (1 - a) + v1 * a;
|
|
}
|
|
|
|
inline Point mix(const Point& v0, const Point& v1, float a)
|
|
{
|
|
return {
|
|
mix(v0.x, v1.x, a),
|
|
mix(v0.y, v1.y, a),
|
|
mix(v0.z, v1.z, a),
|
|
};
|
|
}
|
|
|
|
static Point normalize(Point p)
|
|
{
|
|
auto length = std::sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]);
|
|
if (length > 0.) {
|
|
p[0] /= length;
|
|
p[1] /= length;
|
|
p[2] /= length;
|
|
} else {
|
|
p[0] = 1;
|
|
}
|
|
return p;
|
|
}
|
|
|
|
/**
|
|
* Calculate the position and (optionally) the normal on the given edge at the
|
|
* given isolevel.
|
|
*/
|
|
static Point calculateIntersection(const Field& volume, float isoLevel,
|
|
size_t x, size_t y, size_t z, size_t edgeNumber, Point* normal = nullptr)
|
|
{
|
|
auto x1 = x;
|
|
auto x2 = x;
|
|
auto y1 = y;
|
|
auto y2 = y;
|
|
auto z1 = z;
|
|
auto z2 = z;
|
|
|
|
switch (edgeNumber) {
|
|
case 0:
|
|
y2 += 1;
|
|
break;
|
|
case 1:
|
|
y1 += 1;
|
|
x2 += 1;
|
|
y2 += 1;
|
|
break;
|
|
case 2:
|
|
x1 += 1;
|
|
y1 += 1;
|
|
x2 += 1;
|
|
break;
|
|
case 3:
|
|
x1 += 1;
|
|
break;
|
|
case 4:
|
|
z1 += 1;
|
|
y2 += 1;
|
|
z2 += 1;
|
|
break;
|
|
case 5:
|
|
y1 += 1;
|
|
z1 += 1;
|
|
x2 += 1;
|
|
y2 += 1;
|
|
z2 += 1;
|
|
break;
|
|
case 6:
|
|
x1 += 1;
|
|
y1 += 1;
|
|
z1 += 1;
|
|
x2 += 1;
|
|
z2 += 1;
|
|
break;
|
|
case 7:
|
|
x1 += 1;
|
|
z1 += 1;
|
|
z2 += 1;
|
|
break;
|
|
case 8:
|
|
z2 += 1;
|
|
break;
|
|
case 9:
|
|
y1 += 1;
|
|
y2 += 1;
|
|
z2 += 1;
|
|
break;
|
|
case 10:
|
|
x1 += 1;
|
|
y1 += 1;
|
|
x2 += 1;
|
|
y2 += 1;
|
|
z2 += 1;
|
|
break;
|
|
case 11:
|
|
x1 += 1;
|
|
x2 += 1;
|
|
z2 += 1;
|
|
break;
|
|
default:
|
|
assert(false);
|
|
}
|
|
|
|
auto const pos1 = volume.get_point(x1, y1, z1);
|
|
auto const pos2 = volume.get_point(x2, y2, z2);
|
|
auto const val1 = volume.get(x1, y1, z1);
|
|
auto const val2 = volume.get(x2, y2, z2);
|
|
assert(val1 != val2);
|
|
auto const frac = (isoLevel - val1) / (val2 - val1);
|
|
|
|
if (normal) {
|
|
auto grad1 = volume.get_gradient(x1, y1, z1);
|
|
auto grad2 = volume.get_gradient(x2, y2, z2);
|
|
*normal = normalize(mix(grad1, grad2, frac));
|
|
}
|
|
|
|
return mix(pos1, pos2, frac);
|
|
}
|
|
|
|
struct IdPoint {
|
|
size_t number = 0;
|
|
Point point;
|
|
Point normal;
|
|
};
|
|
|
|
struct Triangle {
|
|
size_t pointId[3];
|
|
};
|
|
|
|
/**
|
|
* The marching cubes algorithm as described here:
|
|
* http://paulbourke.net/geometry/polygonise/
|
|
*
|
|
* @param volume The data field
|
|
* @param isoLevel Minimum isoLevel, all values >= isoLevel will contribute to
|
|
* the mesh
|
|
* @param gradient_normals Compute normals based on field gradient. If false,
|
|
* then don't compute normals.
|
|
* @return The iso-surface mesh
|
|
*/
|
|
Mesh march(const Field& volume, float isoLevel, bool gradient_normals)
|
|
{
|
|
auto const xDim = volume.xDim();
|
|
auto const yDim = volume.yDim();
|
|
auto const zDim = volume.zDim();
|
|
|
|
#ifndef PYMOL_OPENMP
|
|
using isocheck_bool_t = bool;
|
|
#else
|
|
// writes to std::vector<bool> are not thread-safe
|
|
using isocheck_bool_t = char;
|
|
#endif
|
|
|
|
// pre-compute isovalue check for better performance
|
|
std::vector<isocheck_bool_t> isocheck(xDim * yDim * zDim);
|
|
|
|
#pragma omp parallel for
|
|
for (int z = 0; z < zDim; ++z) {
|
|
for (size_t y = 0; y < yDim; ++y) {
|
|
auto const offset = xDim * y + xDim * yDim * z;
|
|
for (size_t x = 0; x < xDim; ++x) {
|
|
isocheck[x + offset] = volume.get(x, y, z) < isoLevel;
|
|
}
|
|
}
|
|
}
|
|
|
|
auto const get_isocheck = [&](size_t x, size_t y, size_t z) -> bool {
|
|
return isocheck[x + xDim * y + xDim * yDim * z];
|
|
};
|
|
|
|
auto const xEnd = xDim - 1;
|
|
auto const yEnd = yDim - 1;
|
|
auto const zEnd = zDim - 1;
|
|
|
|
// With OpenMP, we use one triangles vector per thread, and one vertexMap per
|
|
// z-index, and let OpenMP distribute the z-index runs across threads.
|
|
std::vector<std::vector<Triangle>> trianglesVec(1);
|
|
std::vector<std::unordered_map<size_t, IdPoint>> vertexMapVec(1);
|
|
|
|
#ifndef PYMOL_OPENMP
|
|
#define omp_get_thread_num() 0
|
|
#define vertexMappingGet(eid) vertexMapVec[0][eid]
|
|
#else
|
|
#define vertexMappingGet(eid) vertexMapVec[edgeId2z(eid, xDim, yDim)][eid]
|
|
|
|
trianglesVec.resize(omp_get_max_threads());
|
|
vertexMapVec.resize(zDim);
|
|
|
|
#pragma omp parallel for
|
|
#endif
|
|
for (int z = 0; z < zEnd; ++z) {
|
|
auto& triangles = trianglesVec[omp_get_thread_num()];
|
|
|
|
for (size_t y = 0; y < yEnd; ++y) {
|
|
for (size_t x = 0; x < xEnd; ++x) {
|
|
size_t tableIndex = 0;
|
|
if (get_isocheck(x, y, z))
|
|
tableIndex |= 1;
|
|
if (get_isocheck(x, y + 1, z))
|
|
tableIndex |= 2;
|
|
if (get_isocheck(x + 1, y + 1, z))
|
|
tableIndex |= 4;
|
|
if (get_isocheck(x + 1, y, z))
|
|
tableIndex |= 8;
|
|
if (get_isocheck(x, y, z + 1))
|
|
tableIndex |= 16;
|
|
if (get_isocheck(x, y + 1, z + 1))
|
|
tableIndex |= 32;
|
|
if (get_isocheck(x + 1, y + 1, z + 1))
|
|
tableIndex |= 64;
|
|
if (get_isocheck(x + 1, y, z + 1))
|
|
tableIndex |= 128;
|
|
|
|
if (EDGE_TABLE[tableIndex] == 0) {
|
|
continue;
|
|
}
|
|
|
|
auto const storeEdgeVertex = [&](size_t edgeNumber) {
|
|
if ((EDGE_TABLE[tableIndex] & (1 << edgeNumber)) != 0) {
|
|
auto eid = edgeId(x, y, z, edgeNumber, xDim, yDim);
|
|
auto& edge = vertexMappingGet(eid);
|
|
edge.point = calculateIntersection(volume, isoLevel, x, y, z,
|
|
edgeNumber, gradient_normals ? &edge.normal : nullptr);
|
|
}
|
|
};
|
|
|
|
storeEdgeVertex(3);
|
|
storeEdgeVertex(0);
|
|
storeEdgeVertex(8);
|
|
|
|
if (x == xEnd - 1) {
|
|
storeEdgeVertex(2);
|
|
storeEdgeVertex(11);
|
|
|
|
if (y == yEnd - 1) {
|
|
storeEdgeVertex(10);
|
|
}
|
|
}
|
|
|
|
if (y == yEnd - 1) {
|
|
storeEdgeVertex(1);
|
|
storeEdgeVertex(9);
|
|
|
|
if (z == zEnd - 1) {
|
|
storeEdgeVertex(5);
|
|
}
|
|
}
|
|
|
|
if (z == zEnd - 1) {
|
|
storeEdgeVertex(4);
|
|
storeEdgeVertex(7);
|
|
|
|
if (x == xEnd - 1) {
|
|
storeEdgeVertex(6);
|
|
}
|
|
}
|
|
|
|
auto const* tri_table_row = TRIANGLE_TABLE[tableIndex];
|
|
for (size_t i = 0; tri_table_row[i] != -1; i += 3) {
|
|
auto pointId0 = edgeId(x, y, z, tri_table_row[i], xDim, yDim);
|
|
auto pointId1 = edgeId(x, y, z, tri_table_row[i + 1], xDim, yDim);
|
|
auto pointId2 = edgeId(x, y, z, tri_table_row[i + 2], xDim, yDim);
|
|
triangles.push_back({pointId0, pointId1, pointId2});
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
Mesh mesh;
|
|
for (auto const& vertexMap : vertexMapVec) {
|
|
mesh.vertexCount += vertexMap.size();
|
|
}
|
|
for (auto const& triangles : trianglesVec) {
|
|
mesh.faceCount += triangles.size();
|
|
}
|
|
|
|
mesh.faces.reset(new size_t[mesh.faceCount * 3]);
|
|
mesh.vertices.reset(new Point[mesh.vertexCount]);
|
|
if (gradient_normals) {
|
|
mesh.normals.reset(new Point[mesh.vertexCount]);
|
|
}
|
|
|
|
size_t index = 0;
|
|
for (auto& vertexMap : vertexMapVec) {
|
|
for (auto& pair : vertexMap) {
|
|
mesh.vertices[index] = pair.second.point;
|
|
if (gradient_normals) {
|
|
mesh.normals[index] = pair.second.normal;
|
|
}
|
|
|
|
// index for faces
|
|
pair.second.number = index++;
|
|
}
|
|
}
|
|
|
|
index = 0;
|
|
for (auto const& triangles : trianglesVec) {
|
|
for (const auto& triangle : triangles) {
|
|
mesh.faces[index++] = vertexMappingGet(triangle.pointId[0]).number;
|
|
mesh.faces[index++] = vertexMappingGet(triangle.pointId[1]).number;
|
|
mesh.faces[index++] = vertexMappingGet(triangle.pointId[2]).number;
|
|
}
|
|
}
|
|
|
|
return mesh;
|
|
}
|
|
|
|
/**
|
|
* Calculate triangle-based normals
|
|
*/
|
|
void calculateNormals(Mesh& mesh)
|
|
{
|
|
size_t vertexCount = mesh.vertexCount;
|
|
const Point* vertices = mesh.vertices.get();
|
|
size_t triangleCount = mesh.faceCount;
|
|
const size_t* triangles = mesh.faces.get();
|
|
|
|
mesh.normals.reset(new Point[vertexCount]);
|
|
auto normals = mesh.normals.get();
|
|
|
|
#pragma omp parallel for
|
|
for (int i = 0; i < vertexCount; ++i) {
|
|
normals[i] = {0, 0, 0};
|
|
}
|
|
|
|
#pragma omp parallel for
|
|
for (int i = 0; i < triangleCount; ++i) {
|
|
size_t const id0 = triangles[i * 3];
|
|
size_t const id1 = triangles[i * 3 + 1];
|
|
size_t const id2 = triangles[i * 3 + 2];
|
|
Point const vec1{
|
|
vertices[id1][0] - vertices[id0][0],
|
|
vertices[id1][1] - vertices[id0][1],
|
|
vertices[id1][2] - vertices[id0][2],
|
|
};
|
|
Point const vec2{
|
|
vertices[id2][0] - vertices[id0][0],
|
|
vertices[id2][1] - vertices[id0][1],
|
|
vertices[id2][2] - vertices[id0][2],
|
|
};
|
|
Point const normal{
|
|
vec1[2] * vec2[1] - vec1[1] * vec2[2],
|
|
vec1[0] * vec2[2] - vec1[2] * vec2[0],
|
|
vec1[1] * vec2[0] - vec1[0] * vec2[1],
|
|
};
|
|
#pragma omp critical
|
|
{
|
|
normals[id0][0] += normal[0];
|
|
normals[id0][1] += normal[1];
|
|
normals[id0][2] += normal[2];
|
|
normals[id1][0] += normal[0];
|
|
normals[id1][1] += normal[1];
|
|
normals[id1][2] += normal[2];
|
|
normals[id2][0] += normal[0];
|
|
normals[id2][1] += normal[1];
|
|
normals[id2][2] += normal[2];
|
|
}
|
|
}
|
|
|
|
#pragma omp parallel for
|
|
for (int i = 0; i < vertexCount; ++i) {
|
|
normals[i] = normalize(normals[i]);
|
|
}
|
|
}
|
|
|
|
Point Field::get_gradient(size_t x, size_t y, size_t z) const
|
|
{
|
|
size_t xx[] = {x == 0 ? x : x - 1, std::min(x + 1, xDim() - 1)};
|
|
size_t yy[] = {y == 0 ? y : y - 1, std::min(y + 1, yDim() - 1)};
|
|
size_t zz[] = {z == 0 ? z : z - 1, std::min(z + 1, zDim() - 1)};
|
|
return {
|
|
(get(xx[0], y, z) - get(xx[1], y, z)) / std::max<int>(1, xx[1] - xx[0]),
|
|
(get(x, yy[0], z) - get(x, yy[1], z)) / std::max<int>(1, yy[1] - yy[0]),
|
|
(get(x, y, zz[0]) - get(x, y, zz[1])) / std::max<int>(1, zz[1] - zz[0]),
|
|
};
|
|
}
|
|
|
|
} // namespace mc
|