1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
#include "HarmonicMapper.hpp"
#include <iostream>
#include <math.h>
#include <map>

namespace MeshKit {

HarmonicMapper::HarmonicMapper(MKCore* core, vector<Vertex> &v, vector<Face> &t, vector<Edge> &e, vector<set<int> > a)
{
  mk_core = core;
  vtx.insert(vtx.begin(), v.begin(), v.end());
  tri.insert(tri.begin(), t.begin(), t.end());
  edges.insert(edges.begin(), e.begin(), e.end());
  adj.insert(adj.begin(), a.begin(), a.end());
  
}

void HarmonicMapper::execute()
{
	_iterative_map(0.001);

}

void HarmonicMapper::getUV(vector<Vertex> &v)
{
	int count = -1;
	for (vector<Vertex>::iterator it = vtx.begin(); it != vtx.end(); it++){<--- Prefer prefix ++/-- operators for non-primitive types.
		count++;
		if (it->onBoundary)  continue;
		v[count].uv[0] = it->uv[0];
		v[count].uv[1] = it->uv[1];
		//std::cout << "HarmonicMapper index = " << count << "\tuv = {" << it->uv[0] << "," << it->uv[1] << "}\n";
	}	

}

//matrix method
void HarmonicMapper::_map()<--- The function '_map' is never used.
{
	int n_interior = 0, n_boundary = 0;
	for (vector<Vertex>::iterator it = vtx.begin(); it != vtx.end(); it++){<--- Prefer prefix ++/-- operators for non-primitive types.
		if (it->onBoundary){
			it->id = n_boundary; 			
			n_boundary++;
		}
		else{
			it->id = n_interior;
			n_interior++;
		}
	}
	
	//using Armadillo to solve linear equations
	mat A = zeros<mat>(n_interior, n_interior);
	mat B = zeros<mat>(n_interior, n_boundary);
	//setting up the matrix A

	for (vector<Vertex>::iterator it = vtx.begin(); it != vtx.end(); it++){<--- Prefer prefix ++/-- operators for non-primitive types.
		if (it->onBoundary) continue;
		int index_i = it->id;
		double sw = 0.0;
		for (set<int>::iterator it_adj = adj[it->index].begin(); it_adj != adj[it->index].end(); it_adj++){<--- Prefer prefix ++/-- operators for non-primitive types.
			int index_vtx = -1;
			if (edges[*it_adj].connect[0]->index == it->index)
				index_vtx = edges[*it_adj].connect[1]->index;
			else
				index_vtx = edges[*it_adj].connect[0]->index;
			int index_j = vtx[index_vtx].id;
			double w = edges[*it_adj].e;
			if (vtx[index_vtx].onBoundary)
				B(index_i, index_j) = w;
			else
				A(index_i, index_j) = -1.0*w;
			
			sw += w;
		}
		A(index_i, index_i) = sw;
	}
	
	mat b = zeros<mat>(n_interior, 2);
	
	for (vector<Vertex>::iterator it = vtx.begin(); it != vtx.end(); it++){<--- Prefer prefix ++/-- operators for non-primitive types.
		if (!it->onBoundary)  continue;
		b(it->id,0) = it->uv[0];
		b(it->id,1) = it->uv[1];
	}
	mat c = B * b;//n_interior * 2
    mat uv_coord = solve(A, c);
	
	int count = -1;
	for (vector<Vertex>::iterator it = vtx.begin(); it != vtx.end(); it++){<--- Prefer prefix ++/-- operators for non-primitive types.
		if (it->onBoundary)  continue;
		count++;
		it->uv[0] = uv_coord(count, 0);
		it->uv[1] = uv_coord(count, 1);
	}	
	

}

//iterative method
void HarmonicMapper::_iterative_map(double epsilon)
{
	//move interior vertices to its center of neighbors
	//boundary nodes are set to (0,0)
	vector<int> interior;
	for (std::vector<Vertex>::iterator it = vtx.begin(); it != vtx.end(); it++)<--- Prefer prefix ++/-- operators for non-primitive types.
		if (!it->onBoundary){
			interior.push_back(it->index);
			it->uv[0] = 0.0;
			it->uv[1] = 0.0;
		}
				

	while(true){
		double error = -1.0e+10;		
		
		
		for (std::vector<int>::iterator it = interior.begin(); it != interior.end(); it++){<--- Prefer prefix ++/-- operators for non-primitive types.
			Vector2D uv;
			uv[0] = 0.0;
			uv[1] = 0.0;
			double weight = 0.0;
			for (set<int>::iterator it_e = adj[*it].begin(); it_e != adj[*it].end(); it_e++){<--- Prefer prefix ++/-- operators for non-primitive types.
				int adj_v = -1;
				if (edges[*it_e].connect[0]->index == (*it))
					adj_v = edges[*it_e].connect[1]->index;
				else
					adj_v = edges[*it_e].connect[0]->index;
				uv[0] += edges[*it_e].e*vtx[adj_v].uv[0];
				uv[1] += edges[*it_e].e*vtx[adj_v].uv[1];
				weight += edges[*it_e].e;
			}
			Vector2D pre_uv;
			pre_uv[0] = vtx[*it].uv[0];
			pre_uv[1] = vtx[*it].uv[1];
			vtx[*it].uv[0] = uv[0]/weight;
			vtx[*it].uv[1] = uv[1]/weight;
			double v_err = sqrt(pow(pre_uv[0]-vtx[*it].uv[0],2)+pow(pre_uv[1]-vtx[*it].uv[1],2));
			error = (v_err > error)? v_err : error;
		}
		
		if (error < epsilon) break;
	}
}


HarmonicMapper::~HarmonicMapper()
{
  std::cout << "It is over now in smoothing" << endl;
}

}