OpenGL : Laplacian Smoothing & Taubin Smoothing

Laplacian Smoothing

void Mesh::Laplacian(int iter)
{
	vector<Vec3<double>> newPos;

	for (int it = 0; it < iter; it++) {
		for (auto v : _vertices) {
			auto pos = v->_pos;
			Vec3<double> p;
			for (auto nv : v->_nbVertices)
			{
				p += (nv->_pos - pos) / 2.0f;
			}
			p /= v->_nbVertices.size();
			pos += p;
			newPos.push_back(pos);
		}

		for (int i = 0; i < _vertices.size(); i++)
		{
			auto v = _vertices[i];
			v->_pos = newPos[i];
		}
		newPos.clear();
	}
	computeNormal();
}

Taubin Smoothing

void Mesh::Taubin(int iter)
{
	vector<Vec3<double>> newPos;

	for (int it = 0; it < iter; it++) {
		for (auto v : _vertices) {
			auto pos = v->_pos;
			Vec3<double> p;
			for (auto nv : v->_nbVertices)
			{
				p += (nv->_pos - pos) / 2.0f;
			}
			p /= v->_nbVertices.size();
			pos += p;
			newPos.push_back(pos);
		}

		for (int i = 0; i < _vertices.size(); i++)
		{
			auto v = _vertices[i];
			v->_pos = newPos[i];
		}
		newPos.clear();
	}
	for (int it = 0; it < iter; it++) {
		for (auto v : _vertices) {
			auto pos = v->_pos;
			Vec3<double> p;
			for (auto nv : v->_nbVertices)
			{
				p += (nv->_pos - pos) / 2.0f;
			}
			p /= v->_nbVertices.size();
			p *= -1;
			pos += p;
			newPos.push_back(pos);
		}

		for (int i = 0; i < _vertices.size(); i++)
		{
			auto v = _vertices[i];
			v->_pos = newPos[i];
		}
		newPos.clear();
	}
	computeNormal();
}

 

void Mesh::withCot(int iter)
{
	vector<Vec3<double>> newPos;
	vector<Vec3<double>> overlapPos;
	
	for (int it = 0; it < iter; it++) 
	{
		for (auto v : _vertices) 
		{
			auto pos = v->_pos;
			Vec3<double> p;
			double total_w = 0;
			for (auto nv : v->_nbVertices)
			{
				for (auto fv : v->_nbFaces)
				{
					if (fv->getIndex(nv) != -1)
					{
						int index_v = fv->getIndex(v);
						int index_nv = fv->getIndex(nv);
						overlapPos.push_back(fv->v(3 - index_v - index_nv)->_pos);
					}
				}
				//a,b 계산
				double a1 = sqrt(pow(v->_pos.x() - overlapPos[0].x(), 2) + pow(v->_pos.y() - overlapPos[0].y(), 2) + pow(v->_pos.z() - overlapPos[0].z(), 2));
				double b1 = sqrt(pow(nv->_pos.x() - overlapPos[0].x(), 2) + pow(nv->_pos.y() - overlapPos[0].y(), 2) + pow(nv->_pos.z() - overlapPos[0].z(), 2));
				//가중치 계산
				double theta1 = (v->_pos - overlapPos[0]).Dot(nv->_pos - overlapPos[0]);
				double w1 = acos(theta1 / (a1 * b1));

				double a2 = sqrt(pow(v->_pos.x() - overlapPos[1].x(), 2) + pow(v->_pos.y() - overlapPos[1].y(), 2) + pow(v->_pos.z() - overlapPos[1].z(), 2));
				double b2 = sqrt(pow(nv->_pos.x() - overlapPos[1].x(), 2) + pow(nv->_pos.y() - overlapPos[1].y(), 2) + pow(nv->_pos.z() - overlapPos[1].z(), 2));

				double theta2 = (v->_pos - overlapPos[1]).Dot(nv->_pos - overlapPos[1]);
				double w2 = acos(theta2 / (a2 * b2));

				//가중치를 nv pos에 곱함
				p += (nv->_pos-v->_pos) * 0.5 * ((1 / tan(w1)) + (1 / tan(w2)));
				total_w += 0.5 * ((1 / tan(w1)) + (1 / tan(w2)));

				//overlap 비움
				overlapPos.clear();
			}
			p *= (1 / total_w);
			pos += p;
			newPos.push_back(pos);
		}

		for (int i = 0; i < _vertices.size(); i++)
		{
			auto v = _vertices[i];
			v->_pos = newPos[i];
		}
		newPos.clear();
	}
	computeNormal();
}