3Dヒルベルト曲線をVoxlap KV6形式で出力するツール

hilbert.cpp

#include <cstdio>
#include <vector>
#include <iostream>
#include <array>
#include <cstdlib>
#include <cstdint>
#include <memory>
#include <valarray>
#include <cassert>

namespace
{
	struct voxel_model
	{
		int w, h, d;

		std::valarray<bool> b;

		voxel_model(int w, int h, int d):
		w(w), h(h), d(d) {
			b.resize(w * h * d);
			for(std::size_t i = 0; i < b.size(); i++)
				b[i] = false;
		}

		bool& operator ()(int x, int y, int z) {
			assert(x >= 0); assert(y >= 0); assert(z >= 0);
			assert(x < w); assert(y < h); assert(z < d);
			return b[((x * h) + y) * d + z];
		}

		void write(FILE *f)
		{
			auto writeI32 = [=](std::int32_t i)
			{
				fwrite(&i, 1, 4, f);
			};
			auto writeF32 = [=](float i)
			{
				fwrite(&i, 1, 4, f);
			};
			fwrite("Kvxl", 1, 4, f);
			writeI32(w);
			writeI32(h);
			writeI32(d);
			writeF32(0);
			writeF32(0);
			writeF32(0);

			std::vector<uint32_t> xoffset;
			std::vector<uint16_t> xyoffset;
			xoffset.resize(w);
			xyoffset.resize(w * h);

			struct block_data
			{
				uint32_t color;
				uint16_t zPos;
				uint8_t visFaces;
				uint8_t lighting;
			};	
			static_assert(sizeof(block_data) == 8, "");
			std::vector<block_data> blocks;

			std::size_t pos = 0;
			block_data block { 0xffffff, 0, 0, 0 };
			for(int x = 0; x < w; x++) {
				auto lastPos1 = pos;
				fprintf(stderr, "--- x = %d\n", x);
				for(int y = 0; y < h; y++) {
					auto lastPos2 = pos;
					for(int z = 0; z < d; z++) {
						if ((*this)(x, y, z)) {
							block.zPos = static_cast<uint16_t>(z);
							blocks.push_back(block);
							pos++;
							fprintf(stderr, "*");
						} else { fprintf(stderr, " "); }
					}
					fprintf(stderr, "\n");
					xyoffset[x * h + y] = pos - lastPos2;
				}
				xoffset[x] = pos - lastPos1;
			}

			writeI32(blocks.size());
			fwrite(blocks.data(), sizeof(block), blocks.size(), f);
			fwrite(xoffset.data(), 4, xoffset.size(), f);
			fwrite(xyoffset.data(), 2, xyoffset.size(), f);
		}
	};

	std::unique_ptr<voxel_model> current_model;

	struct vector
	{
		int x, y, z;

		inline vector operator +(const vector& v) {
			return vector {x + v.x, y + v.y, z + v.z};
		}
		inline vector operator -(const vector& v) {
			return vector {x - v.x, y - v.y, z - v.z};
		}
		inline vector operator *(int i) {
			return vector { x * i, y * i, z * i };
		}
		inline vector& operator += (const vector& v) {
			x += v.x; y += v.y; z += v.z;
			return *this;
		}
		inline vector operator -() const {
			return vector { -x, -y, -z };
		}
		int manhattan(const vector& o) const {
			int xx = x - o.x, yy = y - o.y, zz = z - o.z;
			if (xx < 0) xx = -xx;
			if (yy < 0) yy = -yy;
			if (zz < 0) zz = -zz;
			return xx + yy + zz;
		}
	};

	vector dirX {1, 0, 0};
	vector dirY {0, 1, 0};
	vector dirZ {0, 0, 1};
	vector pos {0, 0, 0};

	vector lastPos { -1, -1, -1 };

	void output_pos()
	{
		auto pos = ::pos + dirX + dirY + dirZ;
		std::cout << pos.x << " " << pos.y << " " << pos.z << std::endl;
		if(current_model) {
			//spos.x <<= 1;
			//pos.y <<= 1;
			//pos.z <<= 1;
			if (lastPos.x == -1) lastPos = pos;
			while(lastPos.x > pos.x) {
				(*current_model)(lastPos.x, lastPos.y, lastPos.z) = true;
				lastPos.x--;
			}
			while(lastPos.y > pos.y) {
				(*current_model)(lastPos.x, lastPos.y, lastPos.z) = true;
				lastPos.y--;
			}
			while(lastPos.z > pos.z) {
				(*current_model)(lastPos.x, lastPos.y, lastPos.z) = true;
				lastPos.z--;
			}
			while(lastPos.x < pos.x) {
				(*current_model)(lastPos.x, lastPos.y, lastPos.z) = true;
				lastPos.x++;
			}
			while(lastPos.y < pos.y) {
				(*current_model)(lastPos.x, lastPos.y, lastPos.z) = true;
				lastPos.y++;
			}
			while(lastPos.z < pos.z) {
				(*current_model)(lastPos.x, lastPos.y, lastPos.z) = true;
				lastPos.z++;
			}
			(*current_model)(pos.x, pos.y, pos.z) = true;
		}
	}

	void forward()
	{
		pos += dirZ;
		output_pos();
	}

	void rotate_inner(int& x, int& y, bool b)
	{
		int ox = x, oy = y;
		if (b) {
			x = oy; y = -ox;
		} else {
			x = -oy; y = ox;
		}
	}

	vector local_to_global(vector v)
	{
		return dirX * v.x + dirY * v.y + dirZ * v.z;
	}

	void rotate_global(vector v, vector& rotated)
	{
		if (v.x == 0) {
			if (v.y == 0) {
				if (v.z > 0) {
					rotate_inner(rotated.x, rotated.y, false);
				} else if (v.z < 0) {
					rotate_inner(rotated.x, rotated.y, true);
				} else {
					std::terminate();
				}
			} else if (v.y > 0) {
				rotate_inner(rotated.x, rotated.z, false);
			} else if (v.y < 0) {
				rotate_inner(rotated.x, rotated.z, true);
			} else {
				std::terminate();
			}
		} else if (v.x > 0) {
			rotate_inner(rotated.y, rotated.z, false);
		} else if (v.x < 0) {
			rotate_inner(rotated.y, rotated.z, true);
		} else {
			std::terminate();
		}
	}
	void mirror_x()
	{
		dirX.x = -dirX.x;
		dirY.x = -dirY.x;
		dirZ.x = -dirZ.x;
	}
	void mirror_y()
	{
		dirX.y = -dirX.y;
		dirY.y = -dirY.y;
		dirZ.y = -dirZ.y;
	}
	void mirror_z()
	{
		dirX.z = -dirX.z;
		dirY.z = -dirY.z;
		dirZ.z = -dirZ.z;
	}

	void rotate_global(vector v)
	{
		rotate_global(v, dirX);
		rotate_global(v, dirY);
		rotate_global(v, dirZ);
	}

	void rotate_local(vector v)
	{
		rotate_global(local_to_global(v));
	}

	void production(int level)
	{
		if (level < 0) {
			output_pos();
			return;
		}

		auto save = []
		{ 
			return std::array<vector, 3> {dirX, dirY, dirZ};
		};
		auto restore = [](const std::array<vector, 3>& d)
		{
			dirX = d[0]; dirY = d[1]; dirZ = d[2];
		};

		auto initialPos = pos;
		auto initialDir = save();

		dirX = initialDir[0];
		dirY = initialDir[2];
		dirZ = initialDir[1];

		production(level - 1);
		auto half = 2 << level;//pos.manhattan(initialPos) + 1;
		auto full = half * 2;

		// 1		
		pos = initialPos + initialDir[1] * half;

		dirX = initialDir[2];
		dirY = initialDir[1];
		dirZ = initialDir[0];

		production(level - 1);

		// 2
		pos = initialPos + (initialDir[1] + initialDir[0]) * half;

		dirX = initialDir[0];
		dirY = initialDir[1];
		dirZ = initialDir[2];

		production(level - 1);

		// 3
		pos = initialPos + (initialDir[2] + initialDir[1] + initialDir[0]) * half;

		dirX = -initialDir[1];
		dirY = -initialDir[2];
		dirZ = initialDir[0];

		production(level - 1);

		// 4
		pos = initialPos + initialDir[0] * full + (initialDir[1] + initialDir[2]) * half;

		dirX = -initialDir[1];
		dirY = initialDir[2];
		dirZ = -initialDir[0];

		production(level - 1);

		// 5
		pos = initialPos + (initialDir[2] + initialDir[1] + initialDir[0]) * half;

		dirX = initialDir[0];
		dirY = initialDir[1];
		dirZ = initialDir[2];

		production(level - 1);

		// 6
		pos = initialPos + initialDir[2] * full + (initialDir[1] + initialDir[0]) * half;

		dirX = -initialDir[2];
		dirY = initialDir[1];
		dirZ = -initialDir[0];

		production(level - 1);

		// 7
		pos = initialPos + initialDir[2] * full + (initialDir[1]) * half;

		dirX = initialDir[0];
		dirY = -initialDir[2];
		dirZ = -initialDir[1];

		production(level - 1);

		restore(initialDir);
	}
}

int main()
{
	int level = 4;
	current_model.reset(new voxel_model(4<<level, 4<<level, 4<<level));
	production(level);

	auto *f = std::fopen("/tmp/output.kv6", "wb");
	current_model->write(f);
	std::fclose(f);
}