メッシュをスムージングしてオフセットするプログラム

概要

メッシュ操作の勉強のため、スムージングとオフセットのプログラムを作ってみました。

処理手順

1. エッジサイズが小さすぎる三角形を除外(オフセットすると良くない作用をしそうなため)
2. スムージング (法線をなるべく素直に外側に向けたいため)
3. 頂点を法線方向にオフセット。この時、元のメッシュとの交点を取って、交点があった場合はその点を起点としてオフセット、なかった場合はそのままオフセット (オフセットしたのに元よりも縮んだ箇所をなくすため)

結果

PolyHaven にある carrot_cake_4k をお借りして、オフセット処理を実施しました。 全体的にはオフセットされてますが、クラッカーや穴の側面の窪みが膨らんで干渉してしまっていていまいちな感じとなりました。

使用しているライブラリ

• OpenNURBS (Permissiveな独自ライセンス)

参考サイト

• pmpライブラリ 試したところ、Remeshが結構良い仕事をするようでしたので、使用しているアルゴリズムをここから調べました。

コード

下記、定義してあるのに使っていない関数は作りかけで、怪しい処理になってます。

// Copyright (C) Mocchi (mocchi_2003@yahoo.co.jp)
// License: Boost Software License   See LICENSE.txt for the full license.

#include "opennurbs.h"
#include <cstdio>
#include <cstdint>
#include <unordered_set>
#include <unordered_map>
#include <thread>

#include <windows.h>

float get_ieee754(ON__UINT8 p[4]) {
	return *reinterpret_cast<float *>(p);
}

struct pairhash {
	size_t operator()(const std::pair<uint32_t, uint32_t> &p) const {
		return std::hash<uint64_t>{}(p.first + static_cast<uint64_t>(p.second) * 0x100000000ULL);
	}
};

bool ConvertFromBinarySTL(ON_SimpleArray<ON__UINT8> &stl, ON_Mesh &mesh) {
	mesh.Destroy();
	if (stl.Count() < 84) return false;
	ON__UINT32 num_facets =
		static_cast<ON__UINT32>(stl[80]) + static_cast<ON__UINT32>(stl[81]) * 256 +
		static_cast<ON__UINT32>(stl[82]) * 65536 + static_cast<ON__UINT32>(stl[83]) * 16777216;

	if (stl.Count() < 84 + static_cast<int>(num_facets) * 50 - 2) return false;

	mesh.DestroyDoublePrecisionVertices();
	mesh.SetSinglePrecisionVerticesAsValid();
	ON__UINT8 *p = stl.First() + 84;
	for (ON__UINT32 i = 0; i < num_facets; ++i, p += 50) {
		int vcnt = mesh.m_V.Count();
		mesh.m_FN.AppendNew().Set(get_ieee754(p), get_ieee754(p + 4), get_ieee754(p + 8));
		mesh.m_FN.Last()->Unitize();
		mesh.m_V.AppendNew().Set(get_ieee754(p + 12), get_ieee754(p + 16), get_ieee754(p + 20));
		mesh.m_V.AppendNew().Set(get_ieee754(p + 24), get_ieee754(p + 28), get_ieee754(p + 32));
		mesh.m_V.AppendNew().Set(get_ieee754(p + 36), get_ieee754(p + 40), get_ieee754(p + 44));
		mesh.SetTriangle(mesh.m_F.Count(), vcnt, vcnt + 1, vcnt + 2);
	}
	mesh.CombineIdenticalVertices();
	mesh.Compact();
	return true;
}

struct MeshRayIntersection {
	ON_Mesh *mesh;
	ON_RTree *rtree;
	ON_3dPoint model_center;
	double rough_radius;

	void Initialize(ON_Mesh *mesh_, ON_RTree *rtree_) {
		mesh = mesh_;
		rtree = rtree_;

		ON_BoundingBox tbb;
		mesh_->GetTightBoundingBox(tbb);
		rough_radius = tbb.Diagonal().Length() * 0.5;
		model_center = tbb.Center();
	}
	// 計算結果
	struct Result {
		ON_3dPoint pt;
		int face_idx;
	};
	struct RI_Work {
		MeshRayIntersection *this_;
		Result *result;
		ON_RTreeCapsule cc;
	};
	bool RayIntersection(const ON_3dRay &ray, Result &result) {
		RI_Work w;
		w.this_ = this;
		w.result = &result;
		ON_RTreeCapsule &cc = w.cc;
		result.face_idx = -1;
		for (int i = 0; i < 3; ++i) cc.m_point[0][i] = ray.m_P[i];
		double len_max = rough_radius + ray.m_P.DistanceTo(model_center);
		ON_3dPoint pt = ray.m_P + ray.m_V * len_max;
		for (int i = 0; i < 3; ++i) cc.m_point[1][i] = pt[i];
		cc.m_radius = 0;
		cc.m_domain[0] = 0, cc.m_domain[1] = 1;

		rtree->Search(&cc, IntersectCallback, &w);
		return (result.face_idx >= 0);
	}
	static bool IntersectCallback(void* a_context, ON__INT_PTR a_id) {
		auto *ri = static_cast<RI_Work *>(a_context);
		MeshRayIntersection *this_ = ri->this_;
		ON_Mesh &m = *this_->mesh;
		ON_MeshFace &f = m.m_F[a_id];
		ON_RTreeCapsule &cc = ri->cc;
		ON_3dPoint ptA(m.m_V[f.vi[0]]), ptB(m.m_V[f.vi[1]]), ptC(m.m_V[f.vi[2]]);
		ON_Plane pln(ptA, ptB, ptC);
		double t;
		// 線分に収まっているかどうかの判定
		ON_Line lin(cc.m_point[0], cc.m_point[1]);
		ON_Intersect(lin, pln.plane_equation, &t);
		if (t < 0 || t > 1.0) return true;
		ON_3dPoint ptL = lin.PointAt(t);
		// 三角形と点との内外判定 http://www.sousakuba.com/Programming/gs_hittest_point_triangle.html
		ON_3dVector cpA = ON_CrossProduct(ptA - ptC, ptL - ptA);
		ON_3dVector cpB = ON_CrossProduct(ptB - ptA, ptL - ptB);
		ON_3dVector cpC = ON_CrossProduct(ptC - ptB, ptL - ptC);
		if (ON_DotProduct(cpA, cpB) < 0 || ON_DotProduct(cpA, cpC) < 0) return true;

		ri->result->pt = ptL;
		ri->result->face_idx = a_id;
		return true;
	}
};

void Mesh_NeighborVertices(const ON_Mesh &mesh, ON_ClassArray<ON_SimpleArray<int>> &neighbors_vertices, bool trace_sort = false) {
	int cnt_vertex = mesh.VertexCount();
	ON_ClassArray<std::unordered_set<int>> neighbors_vertices_;
	neighbors_vertices_.SetCapacity(cnt_vertex);
	neighbors_vertices_.SetCount(cnt_vertex);
	for (int j = 0; j < mesh.m_F.Count(); ++j) {
		auto &f = mesh.m_F[j];
		for (int i = 0; i < 3; ++i) {
			int vi_c = f.vi[i], vi_n = f.vi[(i + 1) % 3];
			if (vi_c == vi_n) continue;
			neighbors_vertices_[vi_c].insert(vi_n);
		}
	}
	neighbors_vertices.SetCapacity(cnt_vertex);
	neighbors_vertices.SetCount(cnt_vertex);
	for (int j = 0; j < cnt_vertex; ++j) {
		auto &nv_i = neighbors_vertices_[j];
		if (nv_i.size() == 0) continue;
		auto &nv_o = neighbors_vertices[j];
		nv_o.SetCapacity(nv_i.size());
		if (trace_sort) {
			std::unordered_set<int> v_traced;
			int vi_init = *nv_i.begin();
			int vi_c, vi_n = vi_init;

			for (;;) {
				vi_c = vi_n, vi_n = -1;
				v_traced.insert(vi_c);
				nv_o.Append(vi_c);
				auto &nv_c = neighbors_vertices_[vi_c];

				for (auto iter = nv_c.begin(); iter != nv_c.end(); ++iter) {
					if (nv_i.find(*iter) == nv_i.end()) continue;
					if (v_traced.find(*iter) != v_traced.end()) continue;
					vi_n = *iter;
					v_traced.insert(vi_n);
					break;
				}
				if (vi_n == -1) break;
			}
		} else {
			for (auto iter = nv_i.begin(); iter != nv_i.end(); ++iter) {
				nv_o.Append(*iter);
			}
		}
	}
}

// Reference : M. Botsch, L. Kobbelt,
//    "A Remeshing Approach to Multiresolution Modeling".
//    Eurographics Symposium on Geometry Processing 2004, Pages 189–196
int Remesh_SplitEdges(ON_Mesh &mesh, double edge_length_to_split) {
	int splitted_edges_count = 0;
	std::unordered_map<std::pair<uint32_t, uint32_t>, int, pairhash> edge_splitted;
	for (int k = 0; k < 10; ++k) {
		int splitted_edges_count_prev = splitted_edges_count;
		int face_count = mesh.m_F.Count();
		for (int j = 0; j < face_count; ++j) {
			auto &f = mesh.m_F[j];
			for (int i = 0; i < 3; ++i) {
				int in = (i + 1) % 3;
				std::pair<int, int> edge(f.vi[i], f.vi[in]);
				if (edge.first > edge.second) std::swap(edge.first, edge.second);
				int new_vertex = -1;
				auto iter = edge_splitted.find(edge);
				if (iter != edge_splitted.end()) {
					// 分割済のエッジのため、このFaceは分割する。
					new_vertex = iter->second;
				} else {
					ON_3dPoint p1 = mesh.Vertex(edge.first);
					ON_3dVector edge_dir = mesh.Vertex(edge.second) - p1;
					double edge_length = edge_dir.Length();
					if (edge_length >= edge_length_to_split) {
						new_vertex = mesh.m_V.Count();
						mesh.SetVertex(new_vertex, p1 + edge_dir * 0.5);
						edge_splitted.insert(std::make_pair(edge, new_vertex));
						++splitted_edges_count;
					}
				}


				if (new_vertex >= 0) {
					int prev_vertex = f.vi[in];
					int mate_vertex = f.vi[(i + 2) % 3];
					f.vi[in] = new_vertex;

					auto &nf = mesh.m_F.AppendNew();
					nf.vi[0] = new_vertex;
					nf.vi[1] = prev_vertex;
					nf.vi[2] = mate_vertex;

					// 分割したらこのFaceはそれ以上触らない。
					break;
				}
			}
		}
		if (splitted_edges_count_prev == splitted_edges_count) break;
	}
	return splitted_edges_count;
}

int Remesh_CollapseEdges(ON_Mesh &mesh, double edge_length_to_collapse) {
	int cnt_vtx_orig = mesh.VertexCount();
	int collapsed_edges_count = 0;
	ON_SimpleArray<ON_2dex> edges;
	mesh.GetMeshEdges(edges);
	for (int k = 0; k < 10; ++k) {
		int collapsed_edges_count_prev = collapsed_edges_count;


		for (int j = 0; j < edges.Count(); ++j) {
			ON_3dPoint p1 = mesh.Vertex(edges[j].i);
			ON_3dVector edge_dir = mesh.Vertex(edges[j].j) - p1;
			if (edge_dir.IsZero() || edge_dir.Length() > edge_length_to_collapse) continue;
			ON_3dPoint pm = p1 + edge_dir * 0.5;
			mesh.SetVertex(edges[j].i, pm);
			mesh.SetVertex(edges[j].j, pm);
			++collapsed_edges_count;
		}

		if (collapsed_edges_count_prev == collapsed_edges_count) break;
	}
	mesh.CombineIdenticalVertices(true, true);
	std::printf("Remesh_CollapseEdges: vertices %d => %d\n", cnt_vtx_orig, mesh.VertexCount());
	return collapsed_edges_count;
}

void Remesh_OptimizeValence(ON_Mesh &mesh) {
#if 0
	ON_SimpleArray<int> valences(mesh.VertexCount());
	valences.SetCount(mesh.VertexCount());
	for (int i = 0; i < valences.Count(); ++i) {
		valences[i] = 0;
	}

	for (int j = 0; j < mesh.m_F.Count(); ++j) {
		auto &f = mesh.m_F[j];
		for (int i = 0; i < 3; ++i) {
			int vi_c = f.vi[i], vi_n = f.vi[(i + 1) % 3];
			valences[vi_c]++;
			valences[vi_n]++;
		}
	}

	std::unordered_map<std::pair<int, int>, std::pair<int, int>> edges_to_flip;
	for (int j = 0; j < mesh.m_F.Count(); ++j) {
		auto &f = mesh.m_F[j];
		for (int i = 0; i < 3; ++i) {
			int vi_c = f.vi[i], vi_n = f.vi[(i + 1) % 3];
			valences[vi_c]++;
			valences[vi_n]++;
		}
	}
#endif
}

void Mesh_CalculateVertexArea_BaryCentric(const ON_Mesh &mesh, ON_SimpleArray<double> &bary_centric_area) {
	int cnt_vertex = mesh.VertexCount();
	bary_centric_area.SetCapacity(cnt_vertex);
	bary_centric_area.SetCount(cnt_vertex);

	for (int i = 0; i < cnt_vertex; ++i) {
		bary_centric_area[i] = 0;
	}
	for (int j = 0; j < mesh.m_F.Count(); ++j) {
		ON_3dPoint vtx[3];
		auto &f = mesh.m_F[j];
		for (int i = 0; i < 3; ++i) {
			vtx[i] = mesh.Vertex(f.vi[i]);
		}
		ON_3dVector dir1 = vtx[1] - vtx[0], dir2 = vtx[2] - vtx[0];
		double area = ON_CrossProduct(dir1, dir2).Length() * 0.5;
		double area_3 = area / 3.0;
		for (int i = 0; i < 3; ++i) {
			bary_centric_area[f.vi[i]] += area_3;
		}
	}
}

void Mesh_CalculateVertexArea_Voronoi(const ON_Mesh &mesh, ON_SimpleArray<double> &voronoi_area) {
	int cnt_vertex = mesh.VertexCount();
	voronoi_area.SetCapacity(cnt_vertex);
	voronoi_area.SetCount(cnt_vertex);

	for (int i = 0; i < cnt_vertex; ++i) {
		voronoi_area[i] = 0;
	}

	// ・各三角形に対して、ボロノイ点を求め、三角形の頂点毎にボロノイ面積を求める。
	// ・ボロノイ面積を頂点毎に合算する。
	ON_3dPoint vtx[3];
	ON_3dVector dir_edge[3];

	for (int j = 0; j < mesh.m_F.Count(); ++j) {
		// 中点の平面を求める。
		ON_Plane pln[3];
		auto &f = mesh.m_F[j];
		for (int i = 0; i < 3; ++i) {
			vtx[i] = mesh.Vertex(f.vi[i]);
		}
		for (int i = 0; i < 3; ++i) {
			int in = (i + 1) % 3;
			ON_3dPoint &ps = vtx[i], &pe = vtx[i + 1];
			ON_3dVector dir = pe - ps;
			pln[i].CreateFromNormal(ps + dir * 0.5, dir);

			dir_edge[i] = dir;
		}
		ON_Plane pln_tri(vtx[0], vtx[1], vtx[2]);
		ON_3dPoint pt_voronoi[3];
		for (int i = 0; i < 3; ++i) {
			ON_Intersect(pln[i], pln[(i+1)%3], pln_tri, pt_voronoi[i]);
		}

		for (int i = 0; i < 3; ++i) {
			ON_3dVector dir_voronoi = pt_voronoi[(i+2)%3] - vtx[i];
			double &area = voronoi_area[f.vi[i]];
			area += ON_CrossProduct(dir_edge[i], dir_voronoi).Length() * 0.5;
		}
	}
}

// http://rodolphe-vaillant.fr/entry/33/curvature-of-a-triangle-mesh-definition-and-computation
void Mesh_CalculateCurvature_LaplaceOperator(const ON_Mesh &mesh, ON_SimpleArray<ON_SurfaceCurvature> &kappa) {
	int cnt_vertex = mesh.VertexCount();
	kappa.SetCapacity(cnt_vertex);
	kappa.SetCount(cnt_vertex);

	ON_ClassArray<ON_SimpleArray<int>> neighbors_vertices;
	Mesh_NeighborVertices(mesh, neighbors_vertices, true);

	ON_SimpleArray<double> vertex_area;
	Mesh_CalculateVertexArea_BaryCentric(mesh, vertex_area);

	for (int j = 0; j < cnt_vertex; ++j) {
		ON_3dPoint v_i = mesh.Vertex(j);
		auto &nv = neighbors_vertices[j];
		double sum_rad = 0;
		ON_3dVector v_laplace(0, 0, 0);
		for (int i = 0; i < nv.Count(); ++i) {
			//       v_i
			//        *_
			// dir_b /: ~-_  dir_b
			//      /  :   ~-_
			// v_p *_   :    _* v_n
			//       ~-_ :_-~  
			//   dir_a  ~*   dir_a
			//            v_c
			ON_3dPoint v_jp = mesh.Vertex(i > 0 ? i - 1 : nv.Count() - 1);
			ON_3dPoint v_jc = mesh.Vertex(i);
			ON_3dPoint v_jn = mesh.Vertex(i < nv.Count() - 1 ? i+1 : 0);
			ON_3dVector dir_a[2] = { v_jc - v_jp, v_jc - v_jn };
			ON_3dVector dir_b[2] = { v_i - v_jp, v_i - v_jn };

			//      ^
			//    b/|
			//    / | m       tanθ = m/l,  cotθ = l/m
			//   /θ|
			//  *---+->
			//    l    a
			// l = dot(a,b)
			// m = |b - l*a|
			double cots[2];
			for (int h = 0; h < 2; ++h) {
				dir_a[h].Unitize(), dir_b[h].Unitize();
				double l = ON_DotProduct(dir_a[h], dir_b[h]);
				double m = (dir_b[h] - dir_a[h] * l).Length();
				cots[h] = (m < ON_ZERO_TOLERANCE) ? 0.0 : l / m;
			}
			ON_3dVector dir_c = v_jc - v_i;
			v_laplace += dir_c * (cots[0] + cots[1]);
			dir_c.Unitize();
			double rad = std::atan2(ON_CrossProduct(dir_b[0], dir_c).Length(), ON_DotProduct(dir_b[0], dir_c));
			sum_rad += rad;
		}
		double area = vertex_area[j];
		double k_g = (ON_PI * 2.0 - sum_rad) / area;
		double k_m = v_laplace.Length() / (4.0 * area);
		if (ON_DotProduct(v_laplace, ON_3dVector(mesh.m_N[j])) > 0) k_m *= -1.0;

		double k_root = std::sqrt(k_m * k_m - k_g);
		auto &ks = kappa[j];
		ks.k1 = k_m - k_root;
		ks.k2 = k_m + k_root;
	}

}

void Remesh_TangentialSmoothing(ON_Mesh &mesh, double damping_factor, bool update_normals) {
	int cnt_vertex = mesh.VertexCount();
	ON_SimpleArray<double> vertex_area;

	Mesh_CalculateVertexArea_BaryCentric(mesh, vertex_area);
//	Mesh_CalculateVertexArea_Voronoi(mesh, vertex_area);

	// ・隣接点を記録する。
	ON_ClassArray<ON_SimpleArray<int>> neighbors_vertices;
	Mesh_NeighborVertices(mesh, neighbors_vertices);

	// ・ボロノイ面積で重み付けした重心を求める。
	ON_SimpleArray<ON_3dPoint> centroid(mesh.VertexCount());
	centroid.SetCount(mesh.VertexCount());
	for (int j = 0; j < cnt_vertex; ++j) {
		ON_3dPoint vtx_c = mesh.Vertex(j);
		auto &nv = neighbors_vertices[j];
		ON_3dVector vtx_nsum(0, 0, 0);
		double area_sum = 0;
		for (int i = 0; i < nv.Count(); ++i) {
			int idx = nv[i];
			double va = vertex_area[idx];
			ON_3dPoint vtx_n = mesh.Vertex(idx);
			area_sum += va;
			vtx_nsum += vtx_n * va;
		}
		centroid[j] = vtx_nsum / area_sum;
	}

	// ・重心、法線、damping_factorを使って頂点の位置をずらす。
	for (int j = 0; j < cnt_vertex; ++j) {
		ON_3dPoint vtx = mesh.Vertex(j);
		ON_3dVector dir_diff = centroid[j] - vtx;
		ON_3dVector nrm = mesh.m_N[j];
		for (int i = 0; i < 3; ++i) {
			dir_diff[i] *= (1.0 - nrm[i] * nrm[i]) * damping_factor;
		}
		vtx += dir_diff;
		mesh.SetVertex(j, vtx);
	}
	if (update_normals) {
		mesh.CombineIdenticalVertices(true, true);
		mesh.ComputeVertexNormals();
	}
}

void Remesh_AreaEqualizing(ON_Mesh &mesh, double edge_length) {
	ON_ClassArray<std::unordered_set<int>> neighbors_vertices;

	// edgeを構成する頂点番号は若い順にソートされている
	std::unordered_map<std::pair<uint32_t, uint32_t>, int, pairhash> edge_work; // first: edge, second: new_vertex_id;

	// ・edge_length * 4/3 より長いエッジを2つに分割する
	std::printf("split_edges...\n");
	Remesh_SplitEdges(mesh, edge_length * 4.0 / 3.0);

	std::printf("collapse_edges...\n");
	// ・edge_length * 4/5 より短いエッジを消去する
	Remesh_CollapseEdges(mesh, edge_length * 4.0 / 5.0);
	mesh.Compact();

	// ・valence値が非境界で6、境界で4に近づくようにエッジをフリップ
	//   (エッジを共有する2つの三角形のそのエッジを構成しない残りの頂点同士を
	//    結んだエッジを作成し、元のエッジをそれに置き換える)
	Remesh_OptimizeValence(mesh);

	mesh.ComputeVertexNormals();

	// スムージング
	std::printf("tangential smoothing...\n");
	Remesh_TangentialSmoothing(mesh, 0.2, false);
}

template <typename F> void Mesh_Offset(const ON_Mesh &mesh_i, F &offset, ON_Mesh &mesh_o) {
	mesh_o = mesh_i;
	mesh_o.ComputeVertexNormals();

	// まずはオフセットする。
	for (int i = 0; i < mesh_i.VertexCount(); ++i) {
		ON_3dPoint vtx = mesh_i.Vertex(i);
		auto &nrm = mesh_o.m_N[i];
		vtx += nrm * offset(i);
		mesh_o.SetVertex(i, vtx);
	}

	ON_SimpleArray<ON_2dex> edges;
	mesh_o.GetMeshEdges(edges);

	ON_SimpleArray<bool> vertex_moved(mesh_o.VertexCount());
	vertex_moved.SetCount(mesh_o.VertexCount());

	for (int i = 0; i < vertex_moved.Count(); ++i) {
		vertex_moved[i] = false;
	}

	// オフセットに伴って向きが反転したエッジは、両端の頂点をエッジ中点に移動してエッジ長さ0にする。
	for (int j = 0; j < 10; ++j) {
		bool reversed_exist = false;
		for (int i = 0; i < edges.Count(); ++i) {
			auto &e = edges[i];
			if (vertex_moved[e.i] || vertex_moved[e.j]) continue;
			ON_3dPoint vi = mesh_o.Vertex(e.i), vj = mesh_o.Vertex(e.j);
			ON_3dVector d_o = vj - vi;
			if (d_o.IsZero()) continue;
			ON_3dVector d_i = mesh_i.Vertex(e.j) - mesh_i.Vertex(e.i);
			if (ON_DotProduct(d_i, d_o) < 0) {
				ON_3dPoint vm = (vi + vj) * 0.5;
				mesh_o.SetVertex(e.i, vm);
				mesh_o.SetVertex(e.j, vm);
				reversed_exist = true;
				vertex_moved[e.i] = true;
				vertex_moved[e.j] = true;
			}
		}

		if (!reversed_exist) break;

	}

	// 同一頂点を削除する。
	mesh_o.CombineIdenticalVertices(true, true);
	std::printf("shrink vertices:%d => %d\n", mesh_i.VertexCount(), mesh_o.VertexCount());

	mesh_o.Compact();

}

struct StopWatch {
	LARGE_INTEGER freq, c1, c2;
	double msec;
	StopWatch() {
		::QueryPerformanceFrequency(&freq);
		msec = 0;
	}
	void tic() {
		::QueryPerformanceCounter(&c1);
	}
	void toc() {
		::QueryPerformanceCounter(&c2);
		msec = static_cast<double>(c2.QuadPart - c1.QuadPart) * 1000.0 / static_cast<double>(freq.QuadPart);
	}
};

int main(int argc, char *argv[]) {
	if (argc == 1) return 0;

	ON_Mesh mesh_i, mesh_w, mesh_o;

	{
		FILE *fp = std::fopen(argv[1], "rb");
		std::fseek(fp, 0, SEEK_END);
		size_t sz = std::ftell(fp);
		std::rewind(fp);
		ON_SimpleArray<ON__UINT8> buf(sz);
		buf.SetCount(sz);
		std::fread(buf.First(), 1, sz, fp);
		std::fclose(fp);
		bool rc = ConvertFromBinarySTL(buf, mesh_i);
	}

	StopWatch sw;

	mesh_w = mesh_i;
	sw.tic();
	Remesh_CollapseEdges(mesh_w, 0.05);
	mesh_w.ComputeVertexNormals();
	mesh_w.UnitizeVertexNormals();
	sw.toc();

	std::printf("Remesh_CollapseEdges: %f msec\n", sw.msec);

	sw.tic();
	for (int i = 0; i < 3; ++i)
		Remesh_TangentialSmoothing(mesh_o, 0.5, true);
	sw.toc();

	std::printf("Remesh_TangentialSmoothing: %f msec\n", sw.msec);


	double offset = 0.005;
	sw.tic();
	ON_SimpleArray<double> offset_per_vertex(mesh_w.VertexCount());
	offset_per_vertex.SetCount(mesh_w.VertexCount());
	std::thread threads[4];
	for (int th = 0; th < 4; ++th){
		ON_RTree rtree;
		rtree.CreateMeshFaceTree(&mesh_i);
		MeshRayIntersection mri;
		mri.Initialize(&mesh_i, &rtree);

		threads[th] = std::thread([&, th]{
			for (int i = th; i < mesh_w.VertexCount(); i += 4) {
				double &offset_target = offset_per_vertex[i];
				ON_3dRay ray;
				ray.m_P = mesh_w.Vertex(i);
				ray.m_V = mesh_w.m_N[i];
				MeshRayIntersection::Result result;
				if (!mri.RayIntersection(ray, result)) offset_target = offset;
				double dist = ray.m_P.DistanceTo(result.pt);
				if (dist > 1.0) dist = 1.0;
				offset_target = dist + offset;
			}
		});
	}
	for (int i = 0; i < 4; ++i) threads[i].join();

	sw.toc();
	std::printf("MeshRayIntersection: %f msec\n", sw.msec);

	sw.tic();
	Mesh_Offset(mesh_w, [&](int vtx_idx) {
		return offset_per_vertex[vtx_idx];
	}, mesh_o);
	sw.toc();

	std::printf("Mesh_Offset: %f msec\n", sw.msec);

	ONX_Model model;
	{
		auto &obj = model.m_object_table.AppendNew();
		obj.m_object = &mesh_o;
		obj.m_bDeleteObject = false;
	}

	ON_String filename = argv[1];
	filename += "_offset.3dm";
	model.Write(filename, 4);

	return 0;
}