Sampling / Interpolating multiple 3D arrays (curves)
My goal is to interpolate curves in 2D and 3D space at a predetermined sequential distance to perform PCA on multiple curves.
Suppose a data frame is from multiple 3D arrays (each of different sizes):
>>> df.curves
0 [[0.0, 0.0, 0.91452991453, 0.91452991453, 1.0]...
1 [[0.0, 0.0, 0.734693877551, 0.734693877551, 1....
2 [[0.0, 0.0, 1.0, 1.0, 1.0], [0.0, 0.6435643564...
3 [[0.0, 0.0, 0.551020408163, 0.551020408163, 1....
4 [[0.0, 0.0, 1.0, 1.0, 1.0], [0.0, 0.4389027431...
5 [[0.0, 0.0, 0.734693877551, 0.734693877551, 1....
Name: curves, dtype: object
>>> df.curves[0]
array([[ 0. , 0. , 0.73469388, 0.73469388, 1. ],
[ 0. , 0.1097561 , 0.47560976, 0.5 , 1. ],
[ 1. , 0.65036675, 0.08801956, 0.06845966, 0. ]])
What are the dimensions x
, y
, z
where all the dimensions are the same length, and the size x
and y
increase monotonically:
3D graph
I am trying to sample data equidistantly and allow comparability between curves with the same sampling rate.
A simple sampling function for a 2D curve (no y
dim) would be for a row of a data frame:
def sample2DCurve(row, res=10, method='linear'):
# coords of interpolation
xnew = np.linspace(0, 1, res)
# call scipy interpolator interp1d
# create interpolation function for 2D data
sample2D = interpolate.interp1d(row[0], row[1], kind=method)
# sample data points based on xnew
znew = sample2D(xnew)
return np.array([xnew, znew])
For 3D data, I use path interpolation:
def sample3DCurves(row, res=10, method='linear'):
#npts = row[0].size
#p = np.zeros(npts, dtype=float)
#for i in range(1, npts):
# dx = row[0][i] - row[0][i-1]
# dy = row[1][i] - row[1][i-1]
# dz = row[2][i] - row[2][i-1]
# v = np.array([dx, dy, dz])
# p[i] = p[i-1] + np.linalg.norm(v)
#==============================================================================
# edit: cleaner algebra
x, *y, z = row
# vecs between subsequently measured points
vecs = np.diff(row)
# path: cum distance along points (norm from first to ith point)
path = np.cumsum(np.linalg.norm(vecs, axis=0))
path = np.insert(path, 0, 0)
#==============================================================================
## coords of interpolation
coords = np.linspace(p[0], p[-1], res) #p[0]=0 p[-1]=max(p)
# interpolation func for each axis with the path
sampleX = interpolate.interp1d(p, row[0], kind=method)
sampleY = interpolate.interp1d(p, row[1], kind=method)
sampleZ = interpolate.interp1d(p, row[2], kind=method)
# sample each dim
xnew = sampleX(coords)
ynew = sampleY(coords)
znew = sampleZ(coords)
return np.array([xnew, ynew, znew])
As another approach in 3D, I would like to interpolate along contour lines that form circles in a plane x,y
with a uniform radius:
Circular contours around [0,0,0] in a plane x,y
with three-dimensional intersection
The value is then z
interpolated based on the intersection of the contour lines with the (linearly) interpolated curve projected on the plane x,y
.
But I am struggling to define circular lines and intersect with the projection of the curve / path vectors in the plane x,y
.
Any advice is greatly appreciated! (also in other languages - R / Matlab, etc.)
source to share
For posterity, my crude (simpler and more pythonic code is much appreciated) workaround.
It remains a question / analysis if this is a really useful preprocessing step for pivotal analysis to study the shape of curves.
def seq_sampling(row, res=10, method='linear'):
#3D sequential along x and y (isocircles):
x, y, z = row
# distance to origin for each point (support vectors lengths)
point_distance = np.linalg.norm(row[(0,2),], axis=0)
# isocircle radii
max_radius = math.sqrt(x[-1]+y[-1])
radii = np.linspace(0, max_radius, res)
# last (distance to origin) inner data points per circle (start point of segments)
start_per_radius = [np.max(np.where(point_distance <= radius)) for radius in radii]
# initialize coords
new_x = np.zeros_like(radii)
new_y = np.zeros_like(radii)
new_z = np.zeros_like(radii)
# assign first an last known coordinates
new_x[0], new_y[0] = x[0], y[0] # 0, 0
new_x[-1], new_y[-1] = x[-1], y[-1] # 1, 1
for radius, startpoint in enumerate(start_per_radius[1:-1]):
# intersect circles of radius with corresponding intersecting vectors
# based on https://math.stackexchange.com/questions/311921/get-location-of-vector-circle-intersection
# fix index count (starts with radius > 0)
radius += 1
# span line segment with point O outside and point I inside of iso-circle
endpoint = startpoint+1
O_x = x[endpoint]
O_y = y[endpoint]
I_x = x[startpoint]
I_y = y[startpoint]
# coefficients
a = (O_x-I_x)**2 + (O_y-I_y)**2
b = 2*((O_x-I_x)*(I_y) + (O_y-I_y)*(I_y))
c = (I_x)**2 + (I_y)**2 - radii[radius]**2
# !radicant cannot be zero given:
# each segment is defined by max point lying inside or on iso-circle and the next point
# as both axis are monotonically (strict monotonically y) increasing the next point lies outside of the ico-circle
# thus (in 2D) a segment is intersecting a circle by definition.
t = 2*c / (- b - math.sqrt(b**2 - 4*a*c))
#check if intersection lies on line segment / within boundaries of t [0,1]
if (t >= 0) and (t <= 1):
new_x[radius], new_y[radius] = (O_x - I_x)*t + I_x, (O_y-I_y)*t + I_y
# interpolate new_y based on projected new_y
new_z = interpolate.interp1d(y, z, kind='linear')(new_y[1:-1])
# assign first an last known coordinates
new_z = np.insert(new_z,0,z[0])
new_z = np.append(new_z,z[-1])
return np.array([new_x, new_y, new_z])
source to share