Commit 3df117cac075736484a0b39e6d412d809567560c

Authored by David Mayerich
1 parent 1ae6fbe0

edited spherical harmonics code and added distance field calculator to Python NWT code

Showing 2 changed files with 88 additions and 0 deletions   Show diff stats
python/network.py
... ... @@ -7,6 +7,7 @@ Created on Sat Sep 16 16:34:49 2017
7 7  
8 8 import struct
9 9 import numpy as np
  10 +import scipy as sp
10 11 import networkx as nx
11 12 import matplotlib.pyplot as plt
12 13 import math
... ... @@ -293,3 +294,40 @@ def aabb(nodeList, edgeList):
293 294 if upper[c] < i.p[c]:
294 295 upper[c] = i.p[c]
295 296 return lower, upper
  297 +
  298 +#calculate the distance field at a given resolution
  299 +# R (triple) resolution along each dimension
  300 +def distancefield(nodeList, edgeList, R=(100, 100, 100)):
  301 +
  302 + #get a list of all node positions in the network
  303 + P = []
  304 + for e in edgeList:
  305 + for p in e.points:
  306 + P.append(p)
  307 +
  308 + #turn that list into a Numpy array so that we can create a KD tree
  309 + P = np.array(P)
  310 +
  311 + #generate a KD-Tree out of the network point array
  312 + tree = sp.spatial.cKDTree(P)
  313 +
  314 + plt.scatter(P[:, 0], P[:, 1])
  315 +
  316 + #specify the resolution of the ouput grid
  317 + R = (200, 200, 200)
  318 +
  319 + #generate a meshgrid of the appropriate size and resolution to surround the network
  320 + lower, upper = aabb(nodeList, edgeList) #get the space occupied by the network
  321 + x = np.linspace(lower[0], upper[0], R[0]) #get the grid points for uniform sampling of this space
  322 + y = np.linspace(lower[1], upper[1], R[1])
  323 + z = np.linspace(lower[2], upper[2], R[2])
  324 + X, Y, Z = np.meshgrid(x, y, z)
  325 + #Z = 150 * numpy.ones(X.shape)
  326 +
  327 +
  328 + Q = np.stack((X, Y, Z), 3)
  329 +
  330 +
  331 + D, I = tree.query(Q)
  332 +
  333 + return D
... ...
stim/math/spharmonics.h
... ... @@ -389,6 +389,56 @@ namespace stim {
389 389 }
390 390 }
391 391 }
  392 +
  393 + /// Generate spherical harmonic coefficients based on a set of N samples
  394 + /*void fit(std::vector<stim::vec3<T> > sph_pts, unsigned int L, bool norm = true)
  395 + {
  396 + //std::vector<T> coeffs;
  397 +
  398 + //generate a matrix for fitting
  399 + int B = L*(L+2)+1; //calculate the matrix size
  400 + stim::matrix<T> mat(B, B); //allocate space for the matrix
  401 +
  402 +
  403 +
  404 + std::vector<T> sums;
  405 + //int B = l*(l+2)+1;
  406 + coeffs.resize(B);
  407 + sums.resize(B);
  408 + //stim::matrix<T> mat(B, B);
  409 + for(int i = 0; i < sph_pts.size(); i++)
  410 + {
  411 + mcBegin(l,m);
  412 + mcSample(sph_pts[i][1], sph_pts[i][2], 1.0);
  413 + for(int j = 0; j < B; j++)
  414 + {
  415 + sums[j] += C[j];
  416 + // sums[j] += C[j]*sums[j];
  417 + }
  418 + mcEnd();
  419 + }
  420 + for(int i = 0; i < B; i++)
  421 + {
  422 + for(int j = 0; j < B; j++)
  423 + {
  424 + mat(i,j) = sums[i]*sums[j];
  425 + }
  426 + }
  427 +
  428 + if(mat.det() == 0)
  429 + {
  430 + std::cerr << " matrix not solvable " << std::endl;
  431 + }
  432 + else
  433 + {
  434 + //for(int i = 0; i <
  435 + }
  436 + }*/
  437 +
  438 +
  439 +
  440 +
  441 +
392 442 }; //end class sph_harmonics
393 443  
394 444  
... ...