Skip to content

Instantly share code, notes, and snippets.

@Mathias-Fuchs
Created December 10, 2019 18:28
Show Gist options
  • Save Mathias-Fuchs/265d4885c45004bb3001c5f2b17bf289 to your computer and use it in GitHub Desktop.
Save Mathias-Fuchs/265d4885c45004bb3001c5f2b17bf289 to your computer and use it in GitHub Desktop.
some code quarry for homology computation
using System;
using System.Collections.Generic;
using System.Collections.ObjectModel;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using MathNet.Numerics;
using System.Numerics;
using MathNet.Numerics.Providers.LinearAlgebra.Mkl;
using LibmatCore;
using Rhino.Geometry;
using System.Diagnostics;
namespace LibmatGeometry
{
// https://www.ayasdi.com/blog/topology/measuring-shape/
public class Homology
{
// holds a bunch of stuff
// a tri mesh version
public Mesh _M { get; private set; }
// a custom version of the edges
MeshTopology _MT;
public List<Tuple<int, int>> myEdgeList { get; private set; }
// the most important part: the simplicial complex' boundary operators
public MathNet.Numerics.LinearAlgebra.Double.DenseMatrix d1 { get; private set; }
// much faster than sparse matrix, for some reason
public MathNet.Numerics.LinearAlgebra.Double.DenseMatrix d2 { get; private set; }
// the Betti numbers
public int[] betti { get; private set; }
public String _laprovider;
public MathNet.Numerics.LinearAlgebra.Vector<double>[] _harmonicOneChains;
// public DataFrame _U;
// #if experimental
public double[] conformalCoordinate;
// #endif
// exists on Mac only ..... no comment.
public static Vector3d CrossProduct(Vector3d a, Vector3d b)
{
return new Vector3d(
a.Y * b.Z - a.Z * b.Y,
a.Z * b.X - a.X * b.Z,
a.X * b.Y - a.Y * b.X
);
}
// sets up the matrices computing the simplicial complex
public Homology(Mesh M)
{
// sanitize mesh
Control.UseNativeMKL();
_laprovider = Control.LinearAlgebraProvider.ToString();
M.Weld(0.1);
M.Faces.ConvertQuadsToTriangles();
M.UnifyNormals();
M.Normals.ComputeNormals();
this._M = M;
// Rhino.Geometry.Collections.MeshTopologyEdgeList edgeList = M.TopologyEdges;
_MT = new MeshTopology(this._M, true);
myEdgeList = _MT._list;
// System.Diagnostics.Debug.Assert(myEdgeList.Count == v1); // doesn't hold for big meshes
// set up the boundary operator https://en.wikipedia.org/wiki/Chain_(algebraic_topology)
// https://en.wikipedia.org/wiki/Simplicial_homology
// the boundary of a one-chain is a zero-chain
d1 = new MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(this._M.Vertices.Count, myEdgeList.Count);
d2 = new MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(myEdgeList.Count, this._M.Faces.Count);
for (int edgeIndex = 0; edgeIndex < myEdgeList.Count; edgeIndex++)
{
//Rhino.IndexPair ip = edgeList.GetTopologyVertices(edgeIndex);
var myip = myEdgeList[edgeIndex];
int firstIndex = myip.Item1;
int secondIndex = myip.Item2;
d1[firstIndex, edgeIndex] = 1;
d1[secondIndex, edgeIndex] = -1;
#if extraDebug
System.Diagnostics.Debug.WriteLine("In the first loop, edge index " + edgeIndex.ToString() + " has " +
"contributed a plus one to vertex " + firstIndex.ToString() + " " +
"with coordinates (" +
M.Vertices[firstIndex].X.ToString() + ", " +
M.Vertices[firstIndex].Y.ToString() + ", " +
M.Vertices[firstIndex].Z.ToString() + ") and a minus one to vertex " + secondIndex.ToString() +
" with coordinates (" +
M.Vertices[secondIndex].X.ToString() + ", " +
M.Vertices[secondIndex].Y.ToString() + ", " +
M.Vertices[secondIndex].Z.ToString() + ")."
);
#endif
}
// ok, now we can set up the matrix d2
Rhino.Geometry.Collections.MeshFaceList m = this._M.Faces;
// we need to have a fixed order, so we have to convert the face list to an actual list
List<MeshFace> faces = m.ToList();
for (int faceIndex = 0; faceIndex < this._M.Faces.Count; faceIndex++)
{
MeshFace f = faces[faceIndex];
#if extraDebug
Debug.WriteLine("looking at face " + faceIndex.ToString() +
" with points" +
"(" + _M.Vertices[f.A].X.ToString() + " " + _M.Vertices[f.A].Y.ToString() + " " + _M.Vertices[f.A].Z.ToString() + ") " +
"(" + _M.Vertices[f.B].X.ToString() + " " + _M.Vertices[f.B].Y.ToString() + " " + _M.Vertices[f.B].Z.ToString() + ") " +
"(" + _M.Vertices[f.C].X.ToString() + " " + _M.Vertices[f.C].Y.ToString() + " " + _M.Vertices[f.C].Z.ToString() + ")." +
"(" + _M.Vertices[f.D].X.ToString() + " " + _M.Vertices[f.D].Y.ToString() + " " + _M.Vertices[f.D].Z.ToString() + ").")
;
#endif
// int[] edges = edgeList.GetEdgesForFace(faceIndex); // this one doesn't give correct edges. Bug in Rhino!!
// instead, the following returns the correct code.
var edges = new int[3];
edges[0] = _MT.GetEdgeIndex(f.A, f.B);
edges[1] = _MT.GetEdgeIndex(f.B, f.C);
edges[2] = _MT.GetEdgeIndex(f.A, f.C);
// have to assert all edges have been found
System.Diagnostics.Debug.Assert(edges[0] != -1 && edges[1] != -1 && edges[2] != -1);
foreach (int edgeIndex in edges)
{
// have to keep track of the parity of the iteration of this loop
// edgeCounter++;
// bool edgeCounterOdd = (edgeCounter % 2 == 1);
// each edge generates an entry in the matrix representing the simplicial homology boundary operator.
// however, it is clumsy to find its sign ... plur or minus one.
// We have to decide whether the orientation of the face induces the same
// orientation on the edge as the one given by its associated indexPair.
// It does induce the right one if the scalar product between the face's normal with the
// crossed product of the vector from the first to the second index with the vector from the first
// to the center of the triangle is positive. Otherwise, it doesn't.
var myip = myEdgeList[edgeIndex];
Point3d firstPointOnEdge = this._M.Vertices[myip.Item1];
Point3d secondPointOnEdge = this._M.Vertices[myip.Item2];
Vector3d firstToSecond = secondPointOnEdge - firstPointOnEdge;
// one could also use the face's third vertex here
Point3d faceCenter = new Point3d(
(_M.Vertices[f.A].X + _M.Vertices[f.B].X + _M.Vertices[f.C].X) / 3,
(_M.Vertices[f.A].Y + _M.Vertices[f.B].Y + _M.Vertices[f.C].Y) / 3,
(_M.Vertices[f.A].Z + _M.Vertices[f.B].Z + _M.Vertices[f.C].Z) / 3
);
Vector3d firstToCenter = faceCenter - firstPointOnEdge;
Vector3d crossProduct = CrossProduct(firstToSecond, firstToCenter);
Vector3d DirectionOfFaceNormal = CrossProduct(
((Point3d)_M.Vertices[f.B] - (Point3d)_M.Vertices[f.A]),
((Point3d)_M.Vertices[f.C] - (Point3d)_M.Vertices[f.A])
);
double scalarProduct = crossProduct * DirectionOfFaceNormal;
// if the scalar product is zero, they are perpendicular which must not happen, by the geometry of the situation.
System.Diagnostics.Debug.Assert(scalarProduct != 0);
bool OrientationsCoincide = (scalarProduct > 0);
// finally, we can decide about the sign of the matrix entry.
int sign = (OrientationsCoincide) ? 1 : -1;
// done! We can finally set the corresponding matrix entry.
d2[edgeIndex, faceIndex] += sign;
#if extraDebug
System.Diagnostics.Debug.WriteLine("In the second loop, face index " +
faceIndex.ToString() + " with center (" +
faceCenter.X.ToString() + ", " +
faceCenter.Y.ToString() + ", " +
faceCenter.Z.ToString() + ") has " +
"contributed a sign " + sign.ToString() + " to edge index " + edgeIndex.ToString() + " " +
"which leads from (" +
firstPointOnEdge.X.ToString() + ", " +
firstPointOnEdge.Y.ToString() + ", " +
firstPointOnEdge.Z.ToString() + ") to (" +
secondPointOnEdge.X.ToString() + ", " +
secondPointOnEdge.Y.ToString() + ", " +
secondPointOnEdge.Z.ToString() + ").");
#endif
}
}
// ok, now we have the de-Rham complex! That's already been the bulk of the work.
// the number of connected components of the mesh must be equal to v0 minus the rank of d1,
// because both are equal to the dimension of the zero-th simplicial homology group with real
// coefficients.
// System.Diagnostics.Debug.Assert(M.DisjointMeshCount == v0 - rank);
// now starts the computation-intensive part
//IRHostSession _session = RHostSession.Create("one");
//TU(d1, _session).Wait();
int d1rank = d1.Rank();
// int d1rank = qrr.Diagonal().Find(g => g < 0.00001).Item1;
int d2rank = d2.Rank();
// int d1nullity = d1.Nullity();
int d1nullity = this.myEdgeList.Count - d1rank;
betti = new int[3];
betti[0] = d1.RowCount - d1rank;
// the dimension of the kernel of d2 is the dimension of Z2 minus the rank
betti[2] = d2.ColumnCount - d2rank;
// dimension of kernel of d1 minus dimension of image of d2
betti[1] = d1nullity - d2rank;
var dSquare = d1 * d2;
double test = dSquare.FrobeniusNorm(); // has to be zero
System.Diagnostics.Debug.Assert(test < 0.01);
// GetHomologyBasis
//var projImaged1T = d1.Transpose() * (d1 * d1.Transpose()).PseudoInverse() * d1;
// var projKerneld1 = MathNet.Numerics.LinearAlgebra.Double.SparseMatrix.CreateDiagonal(this.myEdgeList.Count, this.myEdgeList.Count, 1) - projImaged1T;
//int kerneld1 = projKerneld1.Rank();
//var homBasis = projKerneld1.Range();
//int homBasisCount = homBasis.Length;
//int a = 2;
// we now need an ONB for the null space of d1.
// since ker(A) = orthogonal complement of (image of t(A)), we first build the latter image
// var s = d1.Svd();
// we need to get the d1nullity last columns of s.VT
// var kernelOfd1 = MathNet.Numerics.LinearAlgebra.
// strategy: find a matrix whose column space is the kernel of d1
// then, project each column onto the space orthogonal to the image of d2
// this will not lead outside the kernel of d1
// then, compute the .Range() of the matrix thus obtained. This is a homology basis.
// we can then find a harmonic one-form
// var imaged1T = d1.Transpose().Range();
// need to find its orthogonal complement
var TopologicalEdgeLaplacian = d2 * d2.Transpose() + d1.Transpose() * d1; // formula (3) of http://www3.cs.stonybrook.edu/~gu/publications/pdf/2003/sgp_global_conformal_param.pdf
// var projectionKernelLaplacian = MathNet.Numerics.LinearAlgebra.Double.SparseMatrix.CreateDiagonal(this.myEdgeList.Count, this.myEdgeList.Count, 1) - TopologicalEdgeLaplacian.PseudoInverse() * TopologicalEdgeLaplacian;
// var hom = projectionKernelLaplacian.Range();
// int kkk = projectionKernelLaplacian.Rank();
// int hum = TopologicalEdgeLaplacian.Nullity();
_harmonicOneChains = TopologicalEdgeLaplacian.Kernel();
// out of interest: is it really harmonic, i.e., is it closed and co-closed?
// d1 * _harmonicOneChains;
// d2.Transpose() * _harmonicOneChains;
// we can now set up the linear system of equations for the harmonic one-form.
// thre requirements are: the cotangent-weighted laplacian has to be zero
// then, it is harmonic itself, meaning that
// 1) Closedness
// secondly, it is closed meaning that for every face, its integral around the face is zero.
// 2) Duality
// its pairing with the harmonic one-chains has to be zero (duality)
// 3) Harmonicity
// var cotWeightedLaplacian = new MathNet.Numerics.LinearAlgebra.Double.SparseMatrix(this.myEdgeList.Count);
// this gives a single requirement row
// this gives a requirement for every face
int twog = _harmonicOneChains.Length;
// var Requirement = new MathNet.Numerics.LinearAlgebra.Double.SparseMatrix(_M.Faces.Count + twog + 1, myEdgeList.Count);
// to clad it with zeros
var Requirement = new MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(_M.Faces.Count + twog + _M.Vertices.Count, myEdgeList.Count);
// closedness
// gives #F requirements
int meshCounter = 0;
foreach (MeshFace f in _M.Faces)
{
Requirement[meshCounter, _MT.GetEdgeIndex(f.A, f.B)] = 1;
Requirement[meshCounter, _MT.GetEdgeIndex(f.B, f.C)] = 1;
Requirement[meshCounter, _MT.GetEdgeIndex(f.A, f.C)] = 1;
meshCounter++;
}
// 2) Duality
// gives 2g requirements ... one for each cohomology class of degree 1
for (int j = 0; j < twog; j++)
{
for (int i = 0; i < myEdgeList.Count; i++)
{
Requirement[j + _M.Faces.Count, i] = _harmonicOneChains[j][i];
}
}
// 3) Harmonicity in the sense of being in the kernel of the contangent-weighted Laplacian
// gives a requirement for any vertex
// this part of the requirement matrix corresponds to
MathNet.Numerics.LinearAlgebra.Double.DenseMatrix cotangents = _MT.EdgeCotangentWeight();
for (int j = 0; j < _MT._list.Count; j++)
{
double cotgt = cotangents[_MT._list[j].Item1, _MT._list[j].Item2];
Requirement[_MT._list[j].Item1 + _M.Faces.Count + twog, j] = -cotgt;
Requirement[_MT._list[j].Item2 + _M.Faces.Count + twog, j] = -cotgt;
}
// in total, we get #F + 2g + #V = #E + 2 requirements (because of #V - #E + #F = 2 - 2g, the Euler characteristic), so these are just two requirements too many, which is not too bad.
// int a = Requirement.Nullity();
// we wanna do var x = A.Solve(b); for this, we need b
var b = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(_M.Faces.Count + twog + _M.Vertices.Count);
// try something dual to the first basis vector
b[_M.Faces.Count] = 1;
// var firstConformalOneForm = Requirement.Solve(b);
// var n = Requirement.Nullity();
var someSolution = Requirement.QR().Solve(b);
// #if experimentalvar someSolution = Requirement.QR().Solve(b);
// var someSolution = Requirement.PseudoInverse() * b;
conformalCoordinate = new double[_M.Vertices.Count];
bool[] visited = new bool[_M.Vertices.Count];
for (int i = 0; i < _M.Vertices.Count; i++) visited[i] = false;
conformalCoordinate[0] = 0;
visited[0] = true;
// while there is any unvisited vertex left
// works only if there is only one connected component ....
while (visited.ToList().Any(g => !g)) {
foreach (var e in _MT._enum)
{
if (visited[e.Item1] & !visited[e.Item2])
{
conformalCoordinate[e.Item2] = conformalCoordinate[e.Item1] + someSolution[_MT.GetEdgeIndex(e.Item1, e.Item2)];
visited[e.Item2] = true;
}
if (visited[e.Item2] & !visited[e.Item1])
{
conformalCoordinate[e.Item1] = conformalCoordinate[e.Item2] - someSolution[_MT.GetEdgeIndex(e.Item1, e.Item2)];
visited[e.Item1] = true;
}
}
}
// #endif
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment