eic-opticks/analytic/catmull_rom_spline_2.py

0001 #!/usr/bin/env python
0002 """
0003 catmull_rom_spline_2.py
0004 ==========================
0005
0006 This 2nd implementation factors out the calculation of the CatMull Rom interpolation
0007 weights : as they are the same for every segment.
0008
0009 This means that the interpolation reduces to a single matrix multiply
0010 for each segment.
0011
0012 """
0013
0014 import numpy as np, matplotlib as mp
0015 SIZE = np.array([1280, 720])
0016
0017 def catmull_rom_Mu_( u, a=0.5 ):
0018     """
0019     Derivation of the interpolation basis functions:
0020
0021     * https://www.cs.cmu.edu/~fp/courses/graphics/asst5/catmullRom.pdf
0022     * ~/opticks_refs/derivation_catmullRom_spline.pdf
0023
0024     """
0025     return np.array([
0026              -a*u + 2*a*u*u -a*u*u*u,
0027              1+(a-3)*u*u+(2-a)*u*u*u,
0028              a*u+(3-2*a)*u*u+(a-2)*u*u*u,
0029              -a*u*u+a*u*u*u ])
0030
0031
0032 def catmull_rom_Mu( N=100, a=0.5 ):
0033     """
0034     :param N: number of interpolated points for each segment
0035     :param a: "tension"  0.5 is standard for Catmull-Rom
0036     :return interpolation weights: shape (N, 4)
0037
0038     Every segment (which is controlled by a group of 4 control points)
0039     needs to use the same Mu for its interpolation so invoke this once only
0040     before the loop over segments.
0041     """
0042     return catmull_rom_Mu_(np.linspace(0,1,N), a ).T
0043
0044 def circle_points(n=16):
0045     angles = np.arange(0, 2*np.pi, 2*np.pi/n) # NB does not get to 2pi to avoid repeated point
0046     points = np.zeros( (len(angles), 3) )
0047     points[:,0] = np.cos(angles)
0048     points[:,1] = np.sin(angles)
0049     return points
0050
0051 def square_points():
0052     return np.array([ [0,0,0], [1,0,0], [0,1,0], [1,1,0] ])
0053
0054 class Plot(object):
0055     def __init__(self, points, xpp, title):
0056         fig, ax = mp.pyplot.subplots(figsize=SIZE/100.)
0057         fig.suptitle(title)
0058         ax.set_aspect('equal')
0059         ax.plot( xpp[:,0], xpp[:,1], label="xpp"  )
0060         ax.scatter( points[:,0], points[:,1], label="points" )
0061         ax.legend()
0062         fig.show()
0063
0064 if __name__ == '__main__':
0065     points = circle_points(16)
0066     assert len(points) >= 4
0067     numseg = len(points)-3     # sliding window of 4 control points for each segment
0068     looped = True
0069     if looped: numseg += 3
0070
0071     N = 100                           # number of points to interpolate each segment into
0072     Mu = catmull_rom_Mu( N, a=0.5 )   # interpolation weights to apply to each group of 4 control points
0073     print(Mu.shape)
0074
0075     xpp = np.zeros( (numseg*N,3) )
0076     for i in range(numseg):
0077         ## p: (4,3) sliding sub-array of 4 control points that wraps around at the end
0078         p = np.take( points, np.arange(i,i+4), mode='wrap', axis=0 )
0079         xpp[i*N:(i+1)*N] = np.dot( Mu, p )
0080     pass
0081     Plot(points, xpp, "analytic/catmull_rom_spline_2.py" )
0082