From 3df117cac075736484a0b39e6d412d809567560c Mon Sep 17 00:00:00 2001
From: David Mayerich <mayerich@uh.edu>
Date: Wed, 20 Dec 2017 15:02:07 -0600
Subject: [PATCH] edited spherical harmonics code and added distance field calculator to Python NWT code

---
 python/network.py       | 38 ++++++++++++++++++++++++++++++++++++++
 stim/math/spharmonics.h | 50 ++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 88 insertions(+), 0 deletions(-)

diff --git a/python/network.py b/python/network.py
index b14d2b7..dff541f 100644
--- a/python/network.py
+++ b/python/network.py
@@ -7,6 +7,7 @@ Created on Sat Sep 16 16:34:49 2017
 
 import struct
 import numpy as np
+import scipy as sp
 import networkx as nx
 import matplotlib.pyplot as plt
 import math
@@ -293,3 +294,40 @@ def aabb(nodeList, edgeList):
             if upper[c] < i.p[c]:
                 upper[c] = i.p[c]
     return lower, upper
+
+#calculate the distance field at a given resolution
+#   R (triple) resolution along each dimension
+def distancefield(nodeList, edgeList, R=(100, 100, 100)):
+    
+    #get a list of all node positions in the network
+    P = []
+    for e in edgeList:
+        for p in e.points:
+            P.append(p)
+            
+    #turn that list into a Numpy array so that we can create a KD tree
+    P = np.array(P)
+    
+    #generate a KD-Tree out of the network point array
+    tree = sp.spatial.cKDTree(P)
+    
+    plt.scatter(P[:, 0], P[:, 1])
+    
+    #specify the resolution of the ouput grid
+    R = (200, 200, 200)
+    
+    #generate a meshgrid of the appropriate size and resolution to surround the network
+    lower, upper = aabb(nodeList, edgeList)           #get the space occupied by the network
+    x = np.linspace(lower[0], upper[0], R[0])   #get the grid points for uniform sampling of this space
+    y = np.linspace(lower[1], upper[1], R[1])
+    z = np.linspace(lower[2], upper[2], R[2])
+    X, Y, Z = np.meshgrid(x, y, z)
+    #Z = 150 * numpy.ones(X.shape)
+    
+    
+    Q = np.stack((X, Y, Z), 3)
+    
+    
+    D, I = tree.query(Q)
+    
+    return D
diff --git a/stim/math/spharmonics.h b/stim/math/spharmonics.h
index 434e0cb..a09f3eb 100644
--- a/stim/math/spharmonics.h
+++ b/stim/math/spharmonics.h
@@ -389,6 +389,56 @@ namespace stim {
 				}
 			}
 		}
+
+		/// Generate spherical harmonic coefficients based on a set of N samples
+		/*void fit(std::vector<stim::vec3<T> > sph_pts, unsigned int L, bool norm = true)
+		{
+			//std::vector<T> coeffs;
+
+			//generate a matrix for fitting
+			int B = L*(L+2)+1;					//calculate the matrix size
+			stim::matrix<T> mat(B, B);			//allocate space for the matrix
+
+
+
+			std::vector<T> sums;
+			//int B = l*(l+2)+1;
+			coeffs.resize(B);
+			sums.resize(B);
+			//stim::matrix<T> mat(B, B);
+			for(int i = 0; i < sph_pts.size(); i++)
+			{
+				mcBegin(l,m);
+				mcSample(sph_pts[i][1], sph_pts[i][2], 1.0);
+				for(int j = 0; j < B; j++)
+				{
+					sums[j] += C[j];
+					//      sums[j] += C[j]*sums[j];
+				}       
+				mcEnd();
+			}
+			for(int i = 0; i < B; i++)
+			{
+				for(int j = 0; j < B; j++)
+				{
+					mat(i,j) = sums[i]*sums[j];
+				}
+			}
+
+			if(mat.det() == 0)
+			{
+				std::cerr << " matrix not solvable " << std::endl;
+			}
+			else
+			{
+				//for(int i = 0; i <
+			}
+		}*/
+
+
+
+
+
 	};		//end class sph_harmonics
 
 
--
libgit2 0.21.4