【算法】移动立方体实现+vtk导出表面

最新推荐文章于 2023-02-28 10:33:06 发布
最新推荐文章于 2023-02-28 10:33:06 发布
阅读量1.6w 收藏 3
点赞数 1
分类专栏： 算法 文章标签： 移动立方体
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/jin739738709/article/details/122703040
版权

参考：Marching Cubes算法 - 计算机图形学_曾经的大胖的博客-CSDN博客

参考：Polygonising a scalar field (Marching Cubes)

 

下面是源码：

导入自己的nii文件就可以生成表面和表面法向量

nii文件就类似于下图：

感兴趣区域的值为1，背景区域为0 

Vectors.h

#ifndef VECTORS_H
#define VECTORS_H
// File Name: Vectors.h
// Last Modified: 7/8/2000
// Author: Raghavendra Chandrashekara
// Email: rc99@doc.ic.ac.uk, rchandrashekara@hotmail.com
//
// Description: This file contains some useful structures.

typedef float POINT3D[3];
typedef float VECTOR3D[3];

struct POINT3DXYZ {
	float x, y, z;
	friend POINT3DXYZ operator+(const POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2);
	friend POINT3DXYZ operator-(const POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2);
	friend POINT3DXYZ operator*(const POINT3DXYZ& pt3dPoint, float fScale);
	friend POINT3DXYZ operator*(float fScale, const POINT3DXYZ& pt3dPoint);
	friend POINT3DXYZ operator/(const POINT3DXYZ& pt3dPoint, float fScale);
	friend POINT3DXYZ& operator*=(POINT3DXYZ& pt3dPoint, float fScale);
	friend POINT3DXYZ& operator/=(POINT3DXYZ& pt3dPoint, float fScale);
	friend POINT3DXYZ& operator+=(POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2);
	friend POINT3DXYZ& operator-=(POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2);
};

typedef POINT3DXYZ VECTOR3DXYZ;
#endif // VECTORS_H

Vectors.cpp

// File Name: Vectors.h
// Last Modified: 9/8/2000
// Author: Raghavendra Chandrashekara
// Email: rc99@doc.ic.ac.uk, rchandrashekara@hotmail.com
//
// Description: This is the implementation file for POINT3DXYZ class.

//#include "stdafx.h"
#include "Vectors.h"

POINT3DXYZ operator+(const POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2)
{
	POINT3DXYZ result;

	result.x = pt3dPoint1.x + pt3dPoint2.x;
	result.y = pt3dPoint1.y + pt3dPoint2.y;
	result.z = pt3dPoint1.z + pt3dPoint2.z;

	return result;
}

POINT3DXYZ operator-(const POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2)
{
	POINT3DXYZ result;

	result.x = pt3dPoint1.x - pt3dPoint2.x;
	result.y = pt3dPoint1.y - pt3dPoint2.y;
	result.z = pt3dPoint1.z - pt3dPoint2.z;

	return result;
}

POINT3DXYZ operator*(const POINT3DXYZ& pt3dPoint, float fScale)
{
	POINT3DXYZ result;

	result.x = pt3dPoint.x*fScale;
	result.y = pt3dPoint.y*fScale;
	result.z = pt3dPoint.z*fScale;

	return result;
}

POINT3DXYZ operator*(float fScale, const POINT3DXYZ& pt3dPoint)
{
	POINT3DXYZ result;

	result.x = pt3dPoint.x*fScale;
	result.y = pt3dPoint.y*fScale;
	result.z = pt3dPoint.z*fScale;

	return result;
}

POINT3DXYZ operator/(const POINT3DXYZ& pt3dPoint, float fScale)
{
	POINT3DXYZ result;

	result.x = pt3dPoint.x/fScale;
	result.y = pt3dPoint.y/fScale;
	result.z = pt3dPoint.z/fScale;
	
	return result;
}

POINT3DXYZ& operator*=(POINT3DXYZ& pt3dPoint, float fScale)
{
	pt3dPoint.x *= fScale;
	pt3dPoint.y *= fScale;
	pt3dPoint.z *= fScale;

	return pt3dPoint;
}

POINT3DXYZ& operator/=(POINT3DXYZ& pt3dPoint, float fScale)
{
	pt3dPoint.x /= fScale;
	pt3dPoint.y /= fScale;
	pt3dPoint.z /= fScale;

	return pt3dPoint;
}

POINT3DXYZ& operator+=(POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2)
{
	pt3dPoint1.x += pt3dPoint2.x;
	pt3dPoint1.y += pt3dPoint2.y;
	pt3dPoint1.z += pt3dPoint2.z;

	return pt3dPoint1;
}

POINT3DXYZ& operator-=(POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2)
{
	pt3dPoint1.x -= pt3dPoint2.x;
	pt3dPoint1.y -= pt3dPoint2.y;
	pt3dPoint1.z -= pt3dPoint2.z;
	
	return pt3dPoint1;
}

CIsoSurface.h

#ifndef CISOSURFACE_H
#define CISOSURFACE_H
// File Name: CIsoSurface.h
// Last Modified: 5/8/2000
// Author: Raghavendra Chandrashekara (basesd on source code
// provided by Paul Bourke and Cory Gene Bloyd)
// Email: rc99@doc.ic.ac.uk, rchandrashekara@hotmail.com
//
// Description: This is the interface file for the CIsoSurface class.
// CIsoSurface can be used to construct an isosurface from a scalar
// field.

#include <map>
#include <vector>
#include "Vectors.h"

struct POINT3DID {
	unsigned int newID;
	float x, y, z;
};

typedef std::map<unsigned int, POINT3DID> ID2POINT3DID;

struct TRIANGLE {
	unsigned int pointID[3];
};

typedef std::vector<TRIANGLE> TRIANGLEVECTOR;

template <class T> class CIsoSurface {
public:
	// Constructor and destructor.
	CIsoSurface();
	~CIsoSurface();
	
	// Generates the isosurface from the scalar field contained in the
	// buffer ptScalarField[].
	void GenerateSurface(const T* ptScalarField, T tIsoLevel, unsigned int nCellsX, unsigned int nCellsY,  unsigned int nCellsZ, float fCellLengthX, float fCellLengthY, float fCellLengthZ);

	// Returns true if a valid surface has been generated.
	bool IsSurfaceValid();

	// Deletes the isosurface.
	void DeleteSurface();

	// Returns the length, width, and height of the volume in which the
	// isosurface in enclosed in.  Returns -1 if the surface is not
	// valid.
	int GetVolumeLengths(float& fVolLengthX, float& fVolLengthY, float& fVolLengthZ);

	// The vertices which make up the isosurface.
	POINT3D* m_ppt3dVertices;

	// The indices of the vertices which make up the triangles.
	unsigned int* m_piTriangleIndices;

	// The number of vertices which make up the isosurface.
	unsigned int m_nVertices;

	// The number of triangles which make up the isosurface.
	unsigned int m_nTriangles;

protected:

	// The number of normals.
	unsigned int m_nNormals;

	// The normals.
	VECTOR3D* m_pvec3dNormals;

	// List of POINT3Ds which form the isosurface.
	ID2POINT3DID m_i2pt3idVertices;

	// List of TRIANGLES which form the triangulation of the isosurface.
	TRIANGLEVECTOR m_trivecTriangles;

	// Returns the edge ID.
	unsigned int GetEdgeID(unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo);

	// Returns the vertex ID.
	unsigned int GetVertexID(unsigned int nX, unsigned int nY, unsigned int nZ);

	// Calculates the intersection point of the isosurface with an
	// edge.
	POINT3DID CalculateIntersection(unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo);

	// Interpolates between two grid points to produce the point at which
	// the isosurface intersects an edge.
	POINT3DID Interpolate(float fX1, float fY1, float fZ1, float fX2, float fY2, float fZ2, T tVal1, T tVal2);
 
	// Renames vertices and triangles so that they can be accessed more
	// efficiently.
	void RenameVerticesAndTriangles();

	// Calculates the normals.
	void CalculateNormals();

	// No. of cells in x, y, and z directions.
	unsigned int m_nCellsX, m_nCellsY, m_nCellsZ;

	// Cell length in x, y, and z directions.
	float m_fCellLengthX, m_fCellLengthY, m_fCellLengthZ;

	// The buffer holding the scalar field.
	const T* m_ptScalarField;

	// The isosurface value.
	T m_tIsoLevel;

	// Indicates whether a valid surface is present.
	bool m_bValidSurface;

	// Lookup tables used in the construction of the isosurface.
	static const unsigned int m_edgeTable[256];
	static const unsigned int m_triTable[256][16];
};
#endif // CISOSURFACE_H

CIsoSurface.cpp

// File Name: CIsoSurface.cpp
// Last Modified: 5/8/2000
// Author: Raghavendra Chandrashekara (based on source code provided
// by Paul Bourke and Cory Gene Bloyd)
// Email: rc99@doc.ic.ac.uk, rchandrashekara@hotmail.com
//
// Description: This is the implementation file for the CIsoSurface class.

//#include "stdafx.h"
#include <math.h>
#include "CIsoSurface.h"

template <class T> const unsigned int CIsoSurface<T>::m_edgeTable[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
};

template <class T> const unsigned int CIsoSurface<T>::m_triTable[256][16] = {
	{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
	{3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
	{3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
	{3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
	{9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
	{1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
	{9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
	{2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
	{8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
	{9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
	{4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
	{3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
	{1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
	{4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
	{4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
	{9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
	{1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
	{5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
	{2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
	{9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
	{0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
	{2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
	{10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
	{4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
	{5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
	{5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
	{9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
	{0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
	{1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
	{10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
	{8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
	{2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
	{7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
	{9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
	{2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
	{11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
	{9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -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, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
	{1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
	{9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
	{5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
	{2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
	{0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
	{5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
	{6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
	{0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
	{3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
	{6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
	{5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
	{1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
	{10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
	{6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
	{1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
	{8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -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, -1, -1, -1, -1, -1, -1},
	{5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
	{0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -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, -1, -1, -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, -1, -1, -1},
	{10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
	{10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
	{8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
	{1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
	{3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
	{0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
	{10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
	{0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
	{3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -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, -1, -1, -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, -1, -1, -1, -1, -1, -1},
	{6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
	{0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
	{10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
	{10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
	{1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -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, -1, -1, -1, -1, -1, -1},
	{7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -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, -1, -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, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
	{7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
	{10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
	{2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
	{6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
	{7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
	{2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
	{1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
	{10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
	{10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
	{0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
	{7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
	{6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
	{8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
	{9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
	{6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
	{1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
	{4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -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, -1, -1, -1, -1, -1, -1},
	{0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
	{1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
	{8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
	{10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
	{4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
	{10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
	{5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
	{11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
	{9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
	{6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
	{7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -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, -1, -1, -1, -1, -1, -1},
	{9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
	{3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -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, -1, -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, -1, -1, -1},
	{6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
	{3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
	{0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
	{6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
	{1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -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, -1, -1, -1},
	{5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
	{9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -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, -1, -1, -1, -1, -1, -1, -1, -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, -1, -1, -1},
	{0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
	{5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
	{10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
	{11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
	{0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
	{9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -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, -1, -1, -1, -1, -1, -1},
	{8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
	{9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
	{9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
	{1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
	{9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
	{9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
	{5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
	{0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -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, -1, -1, -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, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
	{5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -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, -1, -1, -1},
	{8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
	{0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
	{9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
	{0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
	{1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -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, -1, -1, -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, -1, -1, -1, -1, -1, -1},
	{11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
	{2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -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, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
	{4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
	{4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
	{0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
	{3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
	{3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
	{0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
	{9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
	{1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
	{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}
};

template <class T> CIsoSurface<T>::CIsoSurface()
{
	printf("[jhq] CIsoSurface<T>::CIsoSurface\n");
	m_fCellLengthX = 0;
	m_fCellLengthY = 0;
	m_fCellLengthZ = 0;
	m_nCellsX = 0;
	m_nCellsY = 0;
	m_nCellsZ = 0;
	m_nTriangles = 0;
	m_nNormals = 0;
	m_nVertices = 0;
	m_ppt3dVertices = NULL;
	m_piTriangleIndices = NULL;
	m_pvec3dNormals = NULL;
	m_ptScalarField = NULL;
	m_tIsoLevel = 0;
	m_bValidSurface = false;
}

template <class T> CIsoSurface<T>::~CIsoSurface()
{
	DeleteSurface();
}
#include<iostream>
#include<bitset>
//#define TEST
template <class T> void CIsoSurface<T>::GenerateSurface(const T* ptScalarField, T tIsoLevel, unsigned int nCellsX, unsigned int nCellsY, unsigned int nCellsZ, float fCellLengthX, float fCellLengthY, float fCellLengthZ)
{

	printf("[jhq] CIsoSurface<T>::GenerateSurface\n");


	if (m_bValidSurface)
		DeleteSurface();

	m_tIsoLevel = tIsoLevel;
	m_nCellsX = nCellsX;
	m_nCellsY = nCellsY;
	m_nCellsZ = nCellsZ;
	m_fCellLengthX = fCellLengthX;
	m_fCellLengthY = fCellLengthY;
	m_fCellLengthZ = fCellLengthZ;
	m_ptScalarField = ptScalarField;

	unsigned int nPointsInXDirection = (m_nCellsX + 1);
	unsigned int nPointsInSlice = nPointsInXDirection * (m_nCellsY + 1);

	// Generate isosurface.
	for (unsigned int z = 0; z < m_nCellsZ; z++) {
		for (unsigned int y = 0; y < m_nCellsY; y++) {
			for (unsigned int x = 0; x < m_nCellsX; x++) {
				// Calculate table lookup index from those vertices which are below the isolevel.

				unsigned int tableIndex = 0;
				if (m_ptScalarField[z*nPointsInSlice + y * nPointsInXDirection + x] < m_tIsoLevel)
					tableIndex |= 1;
				if (m_ptScalarField[z*nPointsInSlice + (y + 1)*nPointsInXDirection + x] < m_tIsoLevel)
					tableIndex |= 2;
				if (m_ptScalarField[z*nPointsInSlice + (y + 1)*nPointsInXDirection + (x + 1)] < m_tIsoLevel)
					tableIndex |= 4;
				if (m_ptScalarField[z*nPointsInSlice + y * nPointsInXDirection + (x + 1)] < m_tIsoLevel)
					tableIndex |= 8;
				if (m_ptScalarField[(z + 1)*nPointsInSlice + y * nPointsInXDirection + x] < m_tIsoLevel)
					tableIndex |= 16;
				if (m_ptScalarField[(z + 1)*nPointsInSlice + (y + 1)*nPointsInXDirection + x] < m_tIsoLevel)
					tableIndex |= 32;
				if (m_ptScalarField[(z + 1)*nPointsInSlice + (y + 1)*nPointsInXDirection + (x + 1)] < m_tIsoLevel)
					tableIndex |= 64;
				if (m_ptScalarField[(z + 1)*nPointsInSlice + y * nPointsInXDirection + (x + 1)] < m_tIsoLevel)
					tableIndex |= 128;
#ifdef TEST
					std::cout << "tableIndex: " << std::bitset<8>(tableIndex) << std::endl;
					printf("[jhq] x: %d y: %d z: %d tableIndex: %d\n", x, y, z, tableIndex);
#endif

				// Now create a triangulation of the isosurface in this cell.
				if (m_edgeTable[tableIndex] != 0) {

					if (m_edgeTable[tableIndex] & 8) {
						POINT3DID pt = CalculateIntersection(x, y, z, 3);
						unsigned int id = GetEdgeID(x, y, z, 3);
						m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt));
#ifdef TEST
						printf("[jhq] 3 pt.x: %f pt.y: %f pt.z: %f id: %d\n", pt.x, pt.y, pt.z, id);
#endif
					}
					if (m_edgeTable[tableIndex] & 1) {
						POINT3DID pt = CalculateIntersection(x, y, z, 0);
						unsigned int id = GetEdgeID(x, y, z, 0);
						m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt));
#ifdef TEST
						printf("[jhq] 0  pt.x: %f pt.y: %f pt.z: %f id: %d\n", pt.x, pt.y, pt.z, id);
#endif
					}
					if (m_edgeTable[tableIndex] & 256) {
						POINT3DID pt = CalculateIntersection(x, y, z, 8);
						unsigned int id = GetEdgeID(x, y, z, 8);
						m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt));
#ifdef TEST
						printf("[jhq] 8 pt.x: %f pt.y: %f pt.z: %f id: %d\n", pt.x, pt.y, pt.z, id);
#endif
					}

					if (x == m_nCellsX - 1) {
						if (m_edgeTable[tableIndex] & 4) {
							POINT3DID pt = CalculateIntersection(x, y, z, 2);
							unsigned int id = GetEdgeID(x, y, z, 2);
							m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt));
#ifdef TEST
							printf("[jhq] 2 pt.x: %f pt.y: %f pt.z: %f id: %d\n", pt.x, pt.y, pt.z, id);
#endif
						}
						if (m_edgeTable[tableIndex] & 2048) {
							POINT3DID pt = CalculateIntersection(x, y, z, 11);
							unsigned int id = GetEdgeID(x, y, z, 11);
							m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt));
#ifdef TEST
							printf("[jhq] 11 pt.x: %f pt.y: %f pt.z: %f id: %d\n", pt.x, pt.y, pt.z, id);
#endif
						}
					}
					if (y == m_nCellsY - 1) {
						if (m_edgeTable[tableIndex] & 2) {
							POINT3DID pt = CalculateIntersection(x, y, z, 1);
							unsigned int id = GetEdgeID(x, y, z, 1);
							m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt));
#ifdef TEST
							printf("[jhq] 1 pt.x: %f pt.y: %f pt.z: %f id: %d\n", pt.x, pt.y, pt.z, id);
#endif
						}
						if (m_edgeTable[tableIndex] & 512) {
							POINT3DID pt = CalculateIntersection(x, y, z, 9);
							unsigned int id = GetEdgeID(x, y, z, 9);
							m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt));
#ifdef TEST
							printf("[jhq] 9 pt.x: %f pt.y: %f pt.z: %f id: %d\n", pt.x, pt.y, pt.z, id);
#endif
						}
					}
					if (z == m_nCellsZ - 1) {
						if (m_edgeTable[tableIndex] & 16) {
							POINT3DID pt = CalculateIntersection(x, y, z, 4);
							unsigned int id = GetEdgeID(x, y, z, 4);
							m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt));
#ifdef TEST
							printf("[jhq] 4 pt.x: %f pt.y: %f pt.z: %f id: %d\n", pt.x, pt.y, pt.z, id);
#endif
						}
						if (m_edgeTable[tableIndex] & 128) {
							POINT3DID pt = CalculateIntersection(x, y, z, 7);
							unsigned int id = GetEdgeID(x, y, z, 7);
							m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt));
#ifdef TEST
							printf("[jhq] 7 pt.x: %f pt.y: %f pt.z: %f id: %d\n", pt.x, pt.y, pt.z, id);
#endif
						}
					}
					if ((x == m_nCellsX - 1) && (y == m_nCellsY - 1))
						if (m_edgeTable[tableIndex] & 1024) {
							POINT3DID pt = CalculateIntersection(x, y, z, 10);
							unsigned int id = GetEdgeID(x, y, z, 10);
							m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt));
#ifdef TEST
							printf("[jhq] 10 pt.x: %f pt.y: %f pt.z: %f id: %d\n", pt.x, pt.y, pt.z, id);
#endif
						}
					if ((x == m_nCellsX - 1) && (z == m_nCellsZ - 1))
						if (m_edgeTable[tableIndex] & 64) {
							POINT3DID pt = CalculateIntersection(x, y, z, 6);
							unsigned int id = GetEdgeID(x, y, z, 6);
							m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt));
#ifdef TEST
							printf("[jhq] 6 pt.x: %f pt.y: %f pt.z: %f id: %d\n", pt.x, pt.y, pt.z, id);
#endif
						}
					if ((y == m_nCellsY - 1) && (z == m_nCellsZ - 1))
						if (m_edgeTable[tableIndex] & 32) {
							POINT3DID pt = CalculateIntersection(x, y, z, 5);
							unsigned int id = GetEdgeID(x, y, z, 5);
							m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt));
#ifdef TEST
							printf("[jhq] 5 pt.x: %f pt.y: %f pt.z: %f id: %d\n", pt.x, pt.y, pt.z, id);
#endif
						}

					for (unsigned int i = 0; m_triTable[tableIndex][i] != -1; i += 3) {
						TRIANGLE triangle;
						unsigned int pointID0, pointID1, pointID2;
						pointID0 = GetEdgeID(x, y, z, m_triTable[tableIndex][i]);
						pointID1 = GetEdgeID(x, y, z, m_triTable[tableIndex][i + 1]);
						pointID2 = GetEdgeID(x, y, z, m_triTable[tableIndex][i + 2]);
						triangle.pointID[0] = pointID0;
						triangle.pointID[1] = pointID1;
						triangle.pointID[2] = pointID2;
						m_trivecTriangles.push_back(triangle);
#ifdef TEST
						printf("[jhq] m_triTable: %d %d %d \n", m_triTable[tableIndex][i], m_triTable[tableIndex][i + 1], m_triTable[tableIndex][i + 2]);
						printf("[jhq] triangle: %d %d %d \n", pointID0, pointID1, pointID2);
#endif
					}
				}
			}
		}

	}


	RenameVerticesAndTriangles();
	CalculateNormals();
	m_bValidSurface = true;

}

template <class T> bool CIsoSurface<T>::IsSurfaceValid()
{
	return m_bValidSurface;
}

template <class T> void CIsoSurface<T>::DeleteSurface()
{
	m_fCellLengthX = 0;
	m_fCellLengthY = 0;
	m_fCellLengthZ = 0;
	m_nCellsX = 0;
	m_nCellsY = 0;
	m_nCellsZ = 0;
	m_nTriangles = 0;
	m_nNormals = 0;
	m_nVertices = 0;
	if (m_ppt3dVertices != NULL) {
		delete[] m_ppt3dVertices;
		m_ppt3dVertices = NULL;
	}
	if (m_piTriangleIndices != NULL) {
		delete[] m_piTriangleIndices;
		m_piTriangleIndices = NULL;
	}
	if (m_pvec3dNormals != NULL) {
		delete[] m_pvec3dNormals;
		m_pvec3dNormals = NULL;
	}
	m_ptScalarField = NULL;
	m_tIsoLevel = 0;
	m_bValidSurface = false;
}

template <class T> int CIsoSurface<T>::GetVolumeLengths(float& fVolLengthX, float& fVolLengthY, float& fVolLengthZ)
{
	if (IsSurfaceValid()) {
		fVolLengthX = m_fCellLengthX*m_nCellsX;
		fVolLengthY = m_fCellLengthY*m_nCellsY;
		fVolLengthZ = m_fCellLengthZ*m_nCellsZ;
		return 1;
	}
	else
		return -1;
}

template <class T> unsigned int CIsoSurface<T>::GetEdgeID(unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo)
{
	switch (nEdgeNo) {
	case 0:
		return GetVertexID(nX, nY, nZ) + 1;
	case 1:
		return GetVertexID(nX, nY + 1, nZ);
	case 2:
		return GetVertexID(nX + 1, nY, nZ) + 1;
	case 3:
		return GetVertexID(nX, nY, nZ);
	case 4:
		return GetVertexID(nX, nY, nZ + 1) + 1;
	case 5:
		return GetVertexID(nX, nY + 1, nZ + 1);
	case 6:
		return GetVertexID(nX + 1, nY, nZ + 1) + 1;
	case 7:
		return GetVertexID(nX, nY, nZ + 1);
	case 8:
		return GetVertexID(nX, nY, nZ) + 2;
	case 9:
		return GetVertexID(nX, nY + 1, nZ) + 2;
	case 10:
		return GetVertexID(nX + 1, nY + 1, nZ) + 2;
	case 11:
		return GetVertexID(nX + 1, nY, nZ) + 2;
	default:
		// Invalid edge no.
		return -1;
	}
}

template <class T> unsigned int CIsoSurface<T>::GetVertexID(unsigned int nX, unsigned int nY, unsigned int nZ)
{
	return 3*(nZ*(m_nCellsY + 1)*(m_nCellsX + 1) + nY*(m_nCellsX + 1) + nX);
}

template <class T> POINT3DID CIsoSurface<T>::CalculateIntersection(
	unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo)
{
	float x1, y1, z1, x2, y2, z2;
	unsigned int v1x = nX, v1y = nY, v1z = nZ;
	unsigned int v2x = nX, v2y = nY, v2z = nZ;
	
	switch (nEdgeNo)
	{
	case 0:
		v2y += 1;
		break;
	case 1:
		v1y += 1;
		v2x += 1;
		v2y += 1;
		break;
	case 2:
		v1x += 1;
		v1y += 1;
		v2x += 1;
		break;
	case 3:
		v1x += 1;
		break;
	case 4:
		v1z += 1;
		v2y += 1;
		v2z += 1;
		break;
	case 5:
		v1y += 1;
		v1z += 1;
		v2x += 1;
		v2y += 1;
		v2z += 1;
		break;
	case 6:
		v1x += 1;
		v1y += 1;
		v1z += 1;
		v2x += 1;
		v2z += 1;
		break;
	case 7:
		v1x += 1;
		v1z += 1;
		v2z += 1;
		break;
	case 8:
		v2z += 1;
		break;
	case 9:
		v1y += 1;
		v2y += 1;
		v2z += 1;
		break;
	case 10:
		v1x += 1;
		v1y += 1;
		v2x += 1;
		v2y += 1;
		v2z += 1;
		break;
	case 11:
		v1x += 1;
		v2x += 1;
		v2z += 1;
		break;
	}

	x1 = v1x*m_fCellLengthX;
	y1 = v1y*m_fCellLengthY;
	z1 = v1z*m_fCellLengthZ;
	x2 = v2x*m_fCellLengthX;
	y2 = v2y*m_fCellLengthY;
	z2 = v2z*m_fCellLengthZ;

	unsigned int nPointsInXDirection = (m_nCellsX + 1);
	unsigned int nPointsInSlice = nPointsInXDirection*(m_nCellsY + 1);
	T val1 = m_ptScalarField[v1z*nPointsInSlice + v1y*nPointsInXDirection + v1x];
	T val2 = m_ptScalarField[v2z*nPointsInSlice + v2y*nPointsInXDirection + v2x];
	POINT3DID intersection = Interpolate(x1, y1, z1, x2, y2, z2, val1, val2);
	
	return intersection;
}

bool mycompare(float fX1, float fY1, float fZ1, float fX2, float fY2, float fZ2) {
	if (fX2 < fX1)
		return true;
	else if (fX2 > fX2)
		return false;

	if (fY2 < fY1)
		return true;
	else if (fY2 > fY2)
		return false;

	if (fZ2 < fZ1)
		return true;
	else if (fZ2 > fZ1)
		return false;

	return false;
}

template <class T> POINT3DID CIsoSurface<T>::Interpolate(
	float fX1, float fY1, float fZ1, float fX2, float fY2, float fZ2, T tVal1, T tVal2)
{
#if 1
	POINT3DID interpolation;
	float mu;

	mu = float((m_tIsoLevel - tVal1))/(tVal2 - tVal1);
	interpolation.x = fX1 + mu*(fX2 - fX1);
	interpolation.y = fY1 + mu*(fY2 - fY1);
	interpolation.z = fZ1 + mu*(fZ2 - fZ1);

	return interpolation;
#else
	POINT3DID interpolation;
	float mu;

	mu = float((m_tIsoLevel - tVal1))/(tVal2 - tVal1);

	if (mycompare(fX1, fY1, fZ1, fX2, fY2, fZ2)) {
		float x=fX1, y=fY1, z=fZ1;
		fX1 = fX2;fY1 = fY2;fZ1 = fZ2;
		fX2 = x;fY2 = y;fZ2 = z;
	}
	if (std::abs(tVal1 - tVal2) > 0.00001) {
		interpolation.x = fX1 + mu*(fX2 - fX1);
		interpolation.y = fY1 + mu*(fY2 - fY1);
		interpolation.z = fZ1 + mu*(fZ2 - fZ1);
	}
	else {
		interpolation.x = fX1;
		interpolation.y = fY1;
		interpolation.z = fZ1;
	}


	return interpolation;
#endif
}

template <class T> void CIsoSurface<T>::RenameVerticesAndTriangles()
{
	printf("[jhq] CIsoSurface<T>::RenameVerticesAndTriangles\n");
	unsigned int nextID = 0;
	ID2POINT3DID::iterator mapIterator = m_i2pt3idVertices.begin();
	TRIANGLEVECTOR::iterator vecIterator = m_trivecTriangles.begin();

	// Rename vertices.
	while (mapIterator != m_i2pt3idVertices.end()) {
		(*mapIterator).second.newID = nextID;
		nextID++;
		mapIterator++;
	}

	// Now rename triangles.
	while (vecIterator != m_trivecTriangles.end()) {
		for (unsigned int i = 0; i < 3; i++) {
			unsigned int newID = m_i2pt3idVertices[(*vecIterator).pointID[i]].newID;
			(*vecIterator).pointID[i] = newID;
		}
		vecIterator++;
	}

	// Copy all the vertices and triangles into two arrays so that they
	// can be efficiently accessed.
	// Copy vertices.
	mapIterator = m_i2pt3idVertices.begin();
	m_nVertices = m_i2pt3idVertices.size();
	m_ppt3dVertices = new POINT3D[m_nVertices];
	for (unsigned int i = 0; i < m_nVertices; i++, mapIterator++) {
		m_ppt3dVertices[i][0] = (*mapIterator).second.x;
		m_ppt3dVertices[i][1] = (*mapIterator).second.y;
		m_ppt3dVertices[i][2] = (*mapIterator).second.z;
	}
	// Copy vertex indices which make triangles.
	vecIterator = m_trivecTriangles.begin();
	m_nTriangles = m_trivecTriangles.size();
	m_piTriangleIndices = new unsigned int[m_nTriangles*3];
	for (int i = 0; i < m_nTriangles; i++, vecIterator++) {
		m_piTriangleIndices[i*3] = (*vecIterator).pointID[0];
		m_piTriangleIndices[i*3+1] = (*vecIterator).pointID[1];
		m_piTriangleIndices[i*3+2] = (*vecIterator).pointID[2];
	}

	m_i2pt3idVertices.clear();
	m_trivecTriangles.clear();
}

template <class T> void CIsoSurface<T>::CalculateNormals()
{
	printf("[jhq] CIsoSurface<T>::CalculateNormals\n");
	m_nNormals = m_nVertices;
	m_pvec3dNormals = new VECTOR3D[m_nNormals];
	
	// Set all normals to 0.
	for (unsigned int i = 0; i < m_nNormals; i++) {
		m_pvec3dNormals[i][0] = 0;
		m_pvec3dNormals[i][1] = 0;
		m_pvec3dNormals[i][2] = 0;
	}

	// Calculate normals.
	for (int i = 0; i < m_nTriangles; i++) {
		VECTOR3D vec1, vec2, normal;
		unsigned int id0, id1, id2;
		id0 = m_piTriangleIndices[i*3];
		id1 = m_piTriangleIndices[i*3+1];
		id2 = m_piTriangleIndices[i*3+2];
		vec1[0] = m_ppt3dVertices[id1][0] - m_ppt3dVertices[id0][0];
		vec1[1] = m_ppt3dVertices[id1][1] - m_ppt3dVertices[id0][1];
		vec1[2] = m_ppt3dVertices[id1][2] - m_ppt3dVertices[id0][2];
		vec2[0] = m_ppt3dVertices[id2][0] - m_ppt3dVertices[id0][0];
		vec2[1] = m_ppt3dVertices[id2][1] - m_ppt3dVertices[id0][1];
		vec2[2] = m_ppt3dVertices[id2][2] - m_ppt3dVertices[id0][2];
		normal[0] = vec1[2]*vec2[1] - vec1[1]*vec2[2];
		normal[1] = vec1[0]*vec2[2] - vec1[2]*vec2[0];
		normal[2] = vec1[1]*vec2[0] - vec1[0]*vec2[1];
		m_pvec3dNormals[id0][0] += normal[0];
		m_pvec3dNormals[id0][1] += normal[1];
		m_pvec3dNormals[id0][2] += normal[2];
		m_pvec3dNormals[id1][0] += normal[0];
		m_pvec3dNormals[id1][1] += normal[1];
		m_pvec3dNormals[id1][2] += normal[2];
		m_pvec3dNormals[id2][0] += normal[0];
		m_pvec3dNormals[id2][1] += normal[1];
		m_pvec3dNormals[id2][2] += normal[2];
	}

	// Normalize normals.
	for (int i = 0; i < m_nNormals; i++) {
		float length = sqrt(m_pvec3dNormals[i][0]*m_pvec3dNormals[i][0] + m_pvec3dNormals[i][1]*m_pvec3dNormals[i][1] + m_pvec3dNormals[i][2]*m_pvec3dNormals[i][2]);
		m_pvec3dNormals[i][0] /= length;
		m_pvec3dNormals[i][1] /= length;
		m_pvec3dNormals[i][2] /= length;
	}
}

template class CIsoSurface<short>;
template class CIsoSurface<unsigned short>;
template class CIsoSurface<float>;

main.cpp

#include"CIsoSurface.h"
#include<iostream>

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include <itkNiftiImageIOFactory.h>
#include "itkCastImageFilter.h"
#include "vtkAutoInit.h" 
VTK_MODULE_INIT(vtkRenderingOpenGL2);
VTK_MODULE_INIT(vtkInteractionStyle);

#include "vtkActor.h"
#include "vtkCamera.h"
#include "vtkCellArray.h"   //基元数组类
#include "vtkFloatArray.h"  //浮点型数组类
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkPolyData.h"    //多边形数据类型类
#include "vtkPolyDataMapper.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkRenderer.h"
#include <vtkXMLPolyDataWriter.h>
#include "vtkSmoothPolyDataFilter.h"
#include "vtkPolyDataNormals.h"
#include "vtkWindowedSincPolyDataFilter.h"
#include "vtkMaskPoints.h"
#include "vtkArrowSource.h"
#include "vtkGlyph3D.h"
#include "vtkProperty.h"
constexpr unsigned int Dimension = 3;
typedef itk::Image<short, Dimension>				ShortImageType;
typedef itk::Image<float, Dimension>				FloatImageType;
typedef itk::Image<char, Dimension>					CCharImageType;
typedef itk::Image<unsigned char, Dimension>		UCharImageType;

#if 1
int main(int argc, char ** argv) {

	itk::NiftiImageIOFactory::RegisterOneFactory();

	using CharReaderType = itk::ImageFileReader<FloatImageType>;
	CharReaderType::Pointer char_reader = CharReaderType::New();
	char_reader->SetFileName("test.nii");
	char_reader->Update();
	FloatImageType::Pointer itkImg = char_reader->GetOutput();

	FloatImageType::PixelType* data = itkImg->GetBufferPointer();

	itk::Size<3> size = itkImg->GetLargestPossibleRegion().GetSize();
	std::cout << "size: " << size[0] << " " << size[1] << " " << size[2] << std::endl;

	int totalsize = size[0] * size[1] * size[2];
	std::cout << "totalsize: " << totalsize <<std::endl;

	FloatImageType::SpacingType spacing = itkImg->GetSpacing();
	std::cout << "spacing: " << spacing[0] << " " << spacing[1] << " " << spacing[2] << std::endl;

	itk::Matrix<itk::SpacePrecisionType,3,3> matrix = itkImg->GetDirection();
	std::cout << "matrix: " << matrix[0][0] << " " << matrix[0][1] << " " << matrix[0][2] << std::endl;
	std::cout << "matrix: " << matrix[1][0] << " " << matrix[1][1] << " " << matrix[1][2] << std::endl;
	std::cout << "matrix: " << matrix[2][0] << " " << matrix[2][1] << " " << matrix[2][2] << std::endl;

	itk::Point<itk::SpacePrecisionType,3> origin = itkImg->GetOrigin();
	std::cout << "origin: " << origin[0] << " " << origin[1] << " " << origin[2] << std::endl;

	FloatImageType::PixelType* iter = data;
	for (int i = 0; i < totalsize; i++) {
		*iter *= -1;
		iter++;
	}
	
	CIsoSurface<float> mc;
	mc.GenerateSurface(data, 0, size[0]-1, size[1]-1, size[2]-1, spacing[0], spacing[1], spacing[2]);
	auto m_piTriangleIndices = mc.m_piTriangleIndices;
	auto m_ppt3dVertices = mc.m_ppt3dVertices;
	auto m_nVertices = mc.m_nVertices;
	auto m_nTriangles = mc.m_nTriangles;

	std::cout << "m_nVertices: " << m_nVertices << std::endl;
	//for (unsigned int i = 0; i < m_nVertices; i++) {
	//	std::cout << i << " " << m_ppt3dVertices[i][0] << " " << m_ppt3dVertices[i][1] << " " << m_ppt3dVertices[i][2] << std::endl;
	//}

	std::cout << "m_nTriangles: " << m_nTriangles << std::endl;
	//for (unsigned int i = 0; i < m_nTriangles; i++) {
	//	std::cout << i << " " << m_piTriangleIndices[i*3] << " " << m_piTriangleIndices[i*3+1] << " " << m_piTriangleIndices[i*3+2] << std::endl;
	//}


	// We'll create the building blocks of polydata including data attributes.
	//实例化一个多边形数据对象cube
	vtkPolyData *cube = vtkPolyData::New();
	//实例化一个点对象points
	vtkPoints *points = vtkPoints::New();
	//单元对象
	vtkCellArray *polys = vtkCellArray::New();

	itk::Vector<itk::SpacePrecisionType, 3> vector;
	for (unsigned int i = 0; i < m_nVertices; i++) {
		for (int j = 0; j < 3; j++) vector[j] = m_ppt3dVertices[i][j];
		vector = matrix * vector;

		float pt[3] = {vector[0]+origin[0],vector[1]+origin[1],vector[2]+origin[2]};
		points->InsertPoint(i, pt);
	}
	for (unsigned int i = 0; i < m_nTriangles; i++) {
		vtkIdType triangle[3] = { m_piTriangleIndices[i * 3],m_piTriangleIndices[i * 3 + 1],m_piTriangleIndices[i * 3 + 2] };
		polys->InsertNextCell(3, triangle);
	}

	//与vtkPolyData型数据对象进行关联
	cube->SetPoints(points);//进行点关联
	points->Delete();
	cube->SetPolys(polys);//进行面关联
	polys->Delete();

vtkSmartPointer<vtkPolyDataNormals> normFilter =
        vtkSmartPointer<vtkPolyDataNormals>::New();
    normFilter->SetInputData(cube);
    normFilter->SetComputePointNormals(1);//开启点法向量计算
    normFilter->SetComputeCellNormals(0); //关闭单元法向量计算
    normFilter->SetAutoOrientNormals(1);
    normFilter->SetSplitting(0);
    normFilter->Update();

    vtkSmartPointer<vtkMaskPoints> mask =
        vtkSmartPointer<vtkMaskPoints>::New();
    mask->SetInputData(normFilter->GetOutput());
    mask->SetMaximumNumberOfPoints(3000);
    mask->RandomModeOn();
    mask->Update();

    vtkSmartPointer<vtkArrowSource> arrow =
        vtkSmartPointer<vtkArrowSource>::New();
    arrow->Update(); //一定要更新 否则数据没有添加进来，程序会报错

    vtkSmartPointer<vtkGlyph3D> glyph =
        vtkSmartPointer<vtkGlyph3D>::New();
    glyph->SetInputData(mask->GetOutput());
    glyph->SetSourceData(arrow->GetOutput());//每一点用箭头代替
    glyph->SetVectorModeToUseNormal();//设置向量显示模式和法向量一致
    glyph->SetScaleFactor(5); //设置伸缩比例
    glyph->Update();

    vtkSmartPointer<vtkPolyDataMapper> mapper =
        vtkSmartPointer<vtkPolyDataMapper>::New();
    mapper->SetInputData(cube);
    vtkSmartPointer<vtkPolyDataMapper> normMapper =
        vtkSmartPointer<vtkPolyDataMapper>::New();
    normMapper->SetInputData(normFilter->GetOutput());
    vtkSmartPointer<vtkPolyDataMapper> glyphMapper =
        vtkSmartPointer<vtkPolyDataMapper>::New();
    glyphMapper->SetInputData(glyph->GetOutput());

    vtkSmartPointer<vtkActor> actor =
        vtkSmartPointer<vtkActor>::New();
    actor->SetMapper(mapper);
    vtkSmartPointer<vtkActor> normActor =
        vtkSmartPointer<vtkActor>::New();
    normActor->SetMapper(normMapper);
    vtkSmartPointer<vtkActor> glyphActor =
        vtkSmartPointer<vtkActor>::New();
    glyphActor->SetMapper(glyphMapper);
    glyphActor->GetProperty()->SetColor(1, 0, 0);

    double origView[4] = { 0, 0, 0.33, 1 };
    double normView[4] = { 0.33, 0, 0.66, 1 };
    double glyphView[4] = { 0.66, 0, 1, 1 };
    vtkSmartPointer<vtkRenderer> origRender =
        vtkSmartPointer<vtkRenderer>::New();
    origRender->SetViewport(origView);
    origRender->AddActor(actor);
    origRender->SetBackground(1, 0, 0);
    vtkSmartPointer<vtkRenderer> normRender =
        vtkSmartPointer<vtkRenderer>::New();
    normRender->SetViewport(normView);
    normRender->AddActor(normActor);
    normRender->SetBackground(0, 1, 0);
    vtkSmartPointer<vtkRenderer> glyphRender =
        vtkSmartPointer<vtkRenderer>::New();
    glyphRender->SetViewport(glyphView);
    glyphRender->AddActor(glyphActor);
    glyphRender->AddActor(normActor);
    glyphRender->SetBackground(0, 0, 1);

    vtkSmartPointer<vtkRenderWindow> rw =
        vtkSmartPointer<vtkRenderWindow>::New();
    rw->AddRenderer(origRender);
    rw->AddRenderer(normRender);
    rw->AddRenderer(glyphRender);
    rw->SetWindowName("Calculating Point Norm & Cell Norm");
    rw->SetSize(960, 320);
    rw->Render();

    vtkSmartPointer<vtkRenderWindowInteractor> rwi =
        vtkSmartPointer<vtkRenderWindowInteractor>::New();
    rwi->SetRenderWindow(rw);
    rwi->Initialize();
    rwi->Start();

    return 0;
}

#else
int main(int argc, char ** argv) {

	itk::NiftiImageIOFactory::RegisterOneFactory();

	using CharReaderType = itk::ImageFileReader<FloatImageType>;
	CharReaderType::Pointer char_reader = CharReaderType::New();
	char_reader->SetFileName("test.nii");
	char_reader->Update();
	FloatImageType::Pointer itkImg = char_reader->GetOutput();

	FloatImageType::PixelType* data = itkImg->GetBufferPointer();

	itk::Size<3> size = itkImg->GetLargestPossibleRegion().GetSize();
	std::cout << "size: " << size[0] << " " << size[1] << " " << size[2] << std::endl;

	int totalsize = size[0] * size[1] * size[2];
	std::cout << "totalsize: " << totalsize <<std::endl;

	FloatImageType::SpacingType spacing = itkImg->GetSpacing();
	std::cout << "spacing: " << spacing[0] << " " << spacing[1] << " " << spacing[2] << std::endl;

	itk::Matrix<itk::SpacePrecisionType,3,3> matrix = itkImg->GetDirection();
	std::cout << "matrix: " << matrix[0][0] << " " << matrix[0][1] << " " << matrix[0][2] << std::endl;
	std::cout << "matrix: " << matrix[1][0] << " " << matrix[1][1] << " " << matrix[1][2] << std::endl;
	std::cout << "matrix: " << matrix[2][0] << " " << matrix[2][1] << " " << matrix[2][2] << std::endl;

	itk::Point<itk::SpacePrecisionType,3> origin = itkImg->GetOrigin();
	std::cout << "origin: " << origin[0] << " " << origin[1] << " " << origin[2] << std::endl;

	FloatImageType::PixelType* iter = data;
	for (int i = 0; i < totalsize; i++) {
		*iter *= -1;
		iter++;
	}
	
	//移动立方体MC
	CIsoSurface<float> mc;
	mc.GenerateSurface(data, 0, size[0]-1, size[1]-1, size[2]-1, spacing[0], spacing[1], spacing[2]);



	auto m_piTriangleIndices = mc.m_piTriangleIndices;
	auto m_ppt3dVertices = mc.m_ppt3dVertices;
	auto m_nVertices = mc.m_nVertices;
	auto m_nTriangles = mc.m_nTriangles;

	std::cout << "m_nVertices: " << m_nVertices << std::endl;
	//for (unsigned int i = 0; i < m_nVertices; i++) {
	//	std::cout << i << " " << m_ppt3dVertices[i][0] << " " << m_ppt3dVertices[i][1] << " " << m_ppt3dVertices[i][2] << std::endl;
	//}

	std::cout << "m_nTriangles: " << m_nTriangles << std::endl;
	//for (unsigned int i = 0; i < m_nTriangles; i++) {
	//	std::cout << i << " " << m_piTriangleIndices[i*3] << " " << m_piTriangleIndices[i*3+1] << " " << m_piTriangleIndices[i*3+2] << std::endl;
	//}


	// We'll create the building blocks of polydata including data attributes.
	//实例化一个多边形数据对象cube
	vtkPolyData *cube = vtkPolyData::New();
	//实例化一个点对象points
	vtkPoints *points = vtkPoints::New();
	//单元对象
	vtkCellArray *polys = vtkCellArray::New();

	itk::Vector<itk::SpacePrecisionType, 3> vector;
	for (unsigned int i = 0; i < m_nVertices; i++) {
		for (int j = 0; j < 3; j++) vector[j] = m_ppt3dVertices[i][j];
		vector = matrix * vector;

		float pt[3] = {vector[0]+origin[0],vector[1]+origin[1],vector[2]+origin[2]};
		points->InsertPoint(i, pt);
	}
	for (unsigned int i = 0; i < m_nTriangles; i++) {
		vtkIdType triangle[3] = { m_piTriangleIndices[i * 3],m_piTriangleIndices[i * 3 + 1],m_piTriangleIndices[i * 3 + 2] };
		polys->InsertNextCell(3, triangle);
	}

	//与vtkPolyData型数据对象进行关联
	cube->SetPoints(points);//进行点关联
	points->Delete();
	cube->SetPolys(polys);//进行面关联
	polys->Delete();


#if 0
	vtkSmartPointer<vtkWindowedSincPolyDataFilter> smooth =
			vtkSmartPointer<vtkWindowedSincPolyDataFilter>::New();
	smooth->SetInputData(cube);
    smooth->SetNumberOfIterations(20);
    smooth->Update();
	cube = smooth->GetOutput();
#endif


	//多边形数据送入图像界面进行绘制
	vtkPolyDataMapper *cubeMapper = vtkPolyDataMapper::New();
	cubeMapper->SetInputData(cube);//vtkPolyData数据输出给映射器输入
	cubeMapper->SetScalarRange(0, 7);//设置标量数据范围0~7

	//实例化一个演员cubeActor
	vtkActor *cubeActor = vtkActor::New();
	cubeActor->SetMapper(cubeMapper);

	vtkCamera* camera = vtkCamera::New();
	camera->SetPosition(1, 1, 1);//设置相机位置为（1，1，1）
	camera->SetFocalPoint(0, 0, 0);

	vtkRenderer* renderer = vtkRenderer::New();
	renderer->AddActor(cubeActor);
	renderer->SetActiveCamera(camera);
	renderer->ResetCamera();
	renderer->SetBackground(0, 0, 1);

	vtkRenderWindow *renWin = vtkRenderWindow::New();
	renWin->AddRenderer(renderer);
	renWin->SetSize(800, 800);

	vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
	iren->SetRenderWindow(renWin);

	// interact with data
	renWin->Render();//绘制舞台上的东西
	iren->Start();//开始交互

	vtkSmartPointer<vtkXMLPolyDataWriter> writer =
        vtkSmartPointer<vtkXMLPolyDataWriter>::New();
    writer->SetFileName("test.vtp");
    writer->SetInputData(cube);
    writer->Write();

	// Clean up
	cube->Delete();
	cubeMapper->Delete();
	cubeActor->Delete();
	camera->Delete();
	renderer->Delete();
	renWin->Delete();
	iren->Delete();


	//typedef itk::ImageFileWriter<FloatImageType> WriterType;
	//WriterType::Pointer writer = WriterType::New();
	//writer->SetFileName("output.nii");
	//writer->SetInput(itkImg);
	//writer->Update();




	return 0;
}

#endif

JinSu_
关注
  • 1
    点赞
  • 踩
  • 3
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
VTK移动立方体法处理CT图像建立三维模型
今天写代码了吗
09-04 842
//移动立方体 #include "vtkDICOMImageReader.h" #include "vtkRenderWindowInteractor.h" #include "vtkRenderer.h" #include "vtkRenderWindow.h" #include "vtkMarchingCubes.h" #include "vtkStripper.h" #include "vtkActor.h" #include "vtkPolyDataMapper.h" #include "vtk
Deprecated-paraview-glance-cornerstone:将Cornerstone.js集成到ParaView Glance中的演示。 当前版本是VTK.js和OHIF的集成
05-10
（未维护）ParaView Glance基石 此回购演示了与ParaView Glance的基础集成。 有关的更多信息，请参见 Glance。 注意：此回购已到期。 它基于Glance v1，而开发已经转移到Glance v2。 保留此仓库仅用于演示目的。 功能性： 添加基石注释（椭圆，角度，探针，长度） 查看基石注释状态 基于切片的注释 局限性： 没有叠加支持 不支持在2D视图上叠加几何 建立/运行 # Install prerequsistes npm install # Run development server # You can see the output at http://localhost:9999/ npm run start # Build project for distribution npm run build
三维数据分析、线激光切片数据分析——VTK显示
qq_29912325的博客
06-02 888
1.项目介绍 利用三维点云对线激光采集的数据进行分析，并对其中凹陷部分进行测量。如下为单线分析。主要目标是测量凹槽到下端平面的距离。 2.运用到的知识点 2.1.最小二乘三维平面拟合 数据点云是通过双目相机对打在目标上的线激光匹配得到的结果，也就是说激光的切面不在XOZ或XOY平面内，为了方便数据分析，将三维点云投影到线激光平面。 通过最小二乘拟合，可以通过采样数据点得到平面方程：z=Ax+By+C z = Ax + By + C z=Ax+By+C 公式推导参考：最小二乘解 //m为采样数据数，大于等于
VTK笔记——医学图像等值面提取（Marching Cubes）
热门推荐
Andy的专栏
01-20 1万+
等值面可视化的三种技术，体绘制，移动立方体（Marching Cubes）和部分立方体（Dividing Cubes）. 这篇笔记主要记录的是移动立方体（Marching Cubes）.
练习：VTK 四视图显示（三视图+三正交切片）
LS_JOKER的博客
01-05 2828
主要功能如下：根据三正交切片的位置联动显示其他窗口切片，在三正交切片调整窗宽窗位后其他三视图同步窗宽窗位，并显示切片序号。版本：vtk9.2.0+vs2022。
VTK框选模型表面拾取面片——仅选中前表面
qq_37366618的博客
10-13 2346
VTK框选表面拾取面片——仅选中前表面
Windows10+QT5.15.2+PCL1.13.1+VTK9.2.0
最新发布
01-16
PCL 1.13.1提供了大量用于点云处理的算法，包括滤波、分割、特征提取、表面重建等。这些功能使开发者能够有效地对点云数据进行预处理、分析和建模。PCL库的API设计简洁，易于理解和使用，且与VTK有良好的兼容性，...
PCL1.12+VTK9.1+QT6编译部署
06-10
VTK（Visualization Toolkit）则是一个强大的可视化库，广泛用于科学计算和数据分析中的图形渲染。QT是一个跨平台的开发框架，常用于创建用户界面。本文将详细介绍如何在您的系统上编译和部署PCL 1.12、VTK 9.1以及...
vtk_cube.rar_vtk_vtk纹理立方体
09-21
在本案例中，“vtk_cube.rar_vtk_vtk纹理立方体”是一个针对VTK初学者的教程，通过VC++（Visual C++）环境，演示如何在VTK中创建并显示一个带有纹理的立方体。下面我们将深入探讨VTK库、纹理映射以及如何使用VC++...
pyqt+vtk.rar
12-17
配合本站vtk总结文章网址：https://blog.csdn.net/qq_40155090/article/details/111310107（PyQT5+VTK8加载3D模型（obj、stl、3ds）汇总，以及vtk双相机、获取vtk二维坐标点） 使用
opencv+VTK+PCL.zip
01-04
7.1.1版本的VTK增强了对大数据集的处理能力，支持多种数据类型和高级渲染技术，如体绘制、表面渲染和透明度控制等。VTK不仅提供了强大的可视化工具，还允许用户自定义数据处理和显示流程。 PCL则专注于3D点云处理，...
vtk教程第十二章 应用
qq_39966730的博客
02-28 1448
我们已经描述了一个广泛的可视化技术工具包的设计和实现。在本章中，我们将研究几个案例研究，以展示如何使用这些工具来深入了解重要的应用领域。这些领域包括医学成像、金融可视化、建模、计算流体动力学、有限元分析和算法可视化。对于每一种情况，我们简要描述了问题域，以及我们期望通过可视化获得哪些信息。然后我们设计一种方法来显示结果。很多时候，我们会使用特定于应用程序的工具来扩展可视化工具包的功能。最后，我们给出了一个示例程序并显示了生成的图像。我们所经历的可视化设计过程在每种情况下都是相似的。首先，我们读取或生成特定于
VTK笔记-裁剪分割-裁剪闭合曲面形成闭合多边形数据-vtkClipClosedSurface类
黑山老妖的博客
06-26 3582
vtkClipClosedSurface   vtkClipClosedSurface类使用平面集合剪裁闭合曲面。   vtkClipClosedSurface将使用一组剪裁平面剪裁闭合的多边形曲面。它将通过在剪切输入数据的地方创建新的多边形面来生成新的闭合曲面。   非流体表面不应用作vtkClipClosedSurface的输入。输入曲面不应具有开放边，并且不能具有由两个以上面共享的任何边。  vtkFeatureEdges过滤器可用于验证数据集是否满足这些条件。此外，输入曲面不应自相交，这意味着曲面的
VTK修炼之道55：图形基本操作进阶_表面重建技术（等值面提取）
沈春旭的博客
02-05 5758
1.等值面提取 等值面（线）提取是一种常用的可视化技术，常应用于医学、地质、气象等领域。例如，在医学图像处理中，由于CT、MRI等图像分辨率越来越高，虽然体绘制技术可以清晰地对数据内部结构进行可视化，但是其计算量和效率却制约了其使用。此时可通过等值面提取技术，仅提取感兴趣的一个或者几个组织轮廓，并生成网格模型以供后续的处理和研究。 根据数据类型的不同，VTK中提供了多个等值面提取类，其
vtk表面提取参数研究
技术随笔
12-19 2466
分析vtkMarchingCubes的参数，解决组织表面提取可能存在的平滑效果不好的问题。
获取物体外轮廓，射线探测法。
qq_23158477的博客
05-11 519
ExtractOutsideSurface.cxx #include <vtkActor.h> #include <vtkAppendPolyData.h> #include <vtkCellLocator.h> #include <vtkNamedColors.h> #include <vtkNew.h> #include <vtkPolyDataConnectivityFilter.h> #include <vtkPolyDa
python 使用vtk库生成三维模型
weixin_46579211的博客
06-27 4989
最近在学习python的vtk库，发现相关资料还是比较少的，在这里总结一下本人最近学习的一些经验 使用了The Cancer Imaging Archived的医学图像数据集的标签文件 代码参考了下面的链接 https://kitware.github.io/vtk-examples/site/Python/CompositeData/MultiBlockDataSet/ https://kitware.github.io/vtk-examples/site/Python/Medical/Gen..
VTK学习笔记（五）vtkRotationalExtrusionFilter
yxmlhc的专栏
11-08 732
1. 头文件中的翻译 扫描多边形数据创建“裙子”从自由边和线，线从顶点。 vtkRotationalExtrusionFilter是一个建模过滤器。它以多边形数据作为输入，并在输出时生成多边形数据。输入数据集围绕z轴旋转以创建新的多边形原语。这些原语形成一个“裙子”或扫面。例如，扫线会得到一个圆柱体壳，扫圆会得到一个环面。 这个过滤器有许多控制参数。你可以通过“Capping”实例变量来控制2D对象(即多边形或三角形带)的扫描是否被生成的几何体所覆盖。此外，您还可以...
vtk教程第九章 高级算法
qq_39966730的博客
02-23 2124
我们再次回到可视化算法。本章描述的算法要么实现起来比较复杂，要么在3D可视化应用中应用较少。我们将算法分类为标量、向量、张量或建模算法。9.1标量算法正如我们所看到的，标量算法通常涉及通过查找表映射标量值，或者创建等高线或曲面。在本节中，我们将研究另一种轮廓算法，分割立方体，该算法使用密集点云生成轮廓曲面。我们也描述地毯情节。地毯图不是真正的3D可视化技术，但被广泛用于可视化许多类型的标量数据。最后，裁剪是与轮廓相关的另一个重要算法，其中单元格作为标量值的函数被切割成碎片。分割立方体是一种类似于移动立方体的
python+vtk实现ct医学影像dicom文件体绘制和面绘制三维重建源码
09-07
要使用Python和VTK实现CT医学影像DICOM文件的体绘制和面绘制三维重建，你可以参考以下源码： ``` python import vtk # 创建一个渲染窗口并设置交互方式 renWin = vtk.vtkRenderWindow() iren = vtk.vtkRenderWindowInteractor() iren.SetRenderWindow(renWin) # 读取DICOM文件 reader = vtk.vtkDICOMImageReader() reader.SetDirectoryName("path/to/dicom/files") reader.Update() # 创建体绘制的体素数据集 volumeMapper = vtk.vtkFixedPointVolumeRayCastMapper() volumeMapper.SetInputConnection(reader.GetOutputPort()) # 设置体绘制的颜色和透明度传输函数 volumeProperty = vtk.vtkVolumeProperty() volumeProperty.ShadeOn() volumeProperty.SetColor(vtk.vtkColorTransferFunction()) volumeProperty.SetScalarOpacity(vtk.vtkPiecewiseFunction()) # 创建体绘制的可视化对象 volume = vtk.vtkVolume() volume.SetMapper(volumeMapper) volume.SetProperty(volumeProperty) # 创建面绘制的等值面数据集 contourFilter = vtk.vtkMarchingCubes() contourFilter.SetInputConnection(reader.GetOutputPort()) contourFilter.SetValue(0, thresholdValue) # 设置阈值，提取等值面 # 创建面绘制的Mapper和Actor contourMapper = vtk.vtkPolyDataMapper() contourMapper.SetInputConnection(contourFilter.GetOutputPort()) contourActor = vtk.vtkActor() contourActor.SetMapper(contourMapper) # 创建渲染器和渲染窗口 renderer = vtk.vtkRenderer() renWin.AddRenderer(renderer) renderer.AddActor(volume) renderer.AddActor(contourActor) renderer.SetBackground(0, 0, 0) # 设置背景颜色为黑色 # 设置相机视角 camera = renderer.GetActiveCamera() camera.SetPosition(0, 0, -1) # 设置相机位置 camera.SetFocalPoint(0, 0, 0) # 设置焦点 camera.SetViewUp(0, -1, 0) # 设置视角 # 激活渲染器和交互操作 renderer.ResetCamera() renWin.Render() iren.Start() ``` 请注意，上述代码只提供了一个基本的框架，实际使用时需要根据具体需求进行调整。同时，你需要将代码中的"path/to/dicom/files"替换为实际的DICOM文件路径，并根据需要设置体绘制和面绘制的参数。 希望以上内容对你有所帮助！
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值