eic-opticks/analytic/bezier.py

0001 #!/usr/bin/env python
0002 #
0003 # Copyright (c) 2019 Opticks Team. All Rights Reserved.
0004 #
0005 # This file is part of Opticks
0006 # (see https://bitbucket.org/simoncblyth/opticks).
0007 #
0008 # Licensed under the Apache License, Version 2.0 (the "License");
0009 # you may not use this file except in compliance with the License.
0010 # You may obtain a copy of the License at
0011 #
0012 #   http://www.apache.org/licenses/LICENSE-2.0
0013 #
0014 # Unless required by applicable law or agreed to in writing, software
0015 # distributed under the License is distributed on an "AS IS" BASIS,
0016 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0017 # See the License for the specific language governing permissions and
0018 # limitations under the License.
0019 #
0020
0021 """
0022
0023
0024 Interactive Javascript to pick Bezier points
0025
0026 * http://blog.ivank.net/interpolation-with-cubic-splines.html
0027 * https://www.desmos.com/calculator/cahqdxeshd
0028
0029
0030 Cubic Bezier polynomial interpolating from P0 to P3 as u
0031 goes from 0 to 1 controlled by P1, P2::
0032
0033     B(u) = P0*(1-u)**3 + P1*3*u*(1-u)**2 + P2*3*u**2*(1-u) + P3*u**3
0034
0035 Or more in spirit of Bezier decide on begin/end points and
0036 control points
0037
0038 ::
0039
0040     (z1, rr1)
0041     (cz1, crr1)
0042     (cz2, crr2)
0043     (z2, rr2)
0044
0045
0046 """
0047 import logging
0048
0049 log = logging.getLogger(__name__)
0050 import pylab as plt
0051
0052 import numpy as np
0053 from sympy import symbols, Poly
0054
0055 def make_bezier(t,xs):
0056     if len(xs) == 2:
0057         bz = xs[0]*(1-t) + xs[1]*t
0058     elif len(xs) == 3:
0059         bz = xs[0]*(1-t)**2 + xs[1]*2*t*(1-t) + xs[2]*t**2
0060     elif len(xs) == 4:
0061         bz = xs[0]*(1-t)**3 + xs[1]*3*t*(1-t)**2 + xs[2]*3*t**2*(1-t) +  xs[3]*t**3
0062     else:
0063         assert 0, xs
0064     pass
0065     return bz
0066
0067 class Bezier(object):
0068     """
0069     ety is the cubic in t
0070     but for surface-of-revolution need cubic in z
0071     where:
0072
0073           t = (z - z1)/(z2 - z1)
0074
0075     """
0076     @classmethod
0077     def ZCoeff(cls, xy):
0078         bz = cls(xy)
0079         return bz.zco_
0080
0081     def __init__(self, xy):
0082         assert len(xy) == 4
0083         xy = np.asarray(xy, dtype=np.float32)
0084         x_ = xy[:,0]
0085         y_ = xy[:,1]
0086
0087         t = symbols("t")
0088         x = symbols("x0,x1,x2,x3")
0089         y = symbols("y0,y1,y2,y3")
0090
0091         etx = make_bezier(t,x)
0092         ety = make_bezier(t,y)
0093
0094         etx_ = etx.subs([[x[0],x_[0]],[x[1],x_[1]],[x[2],x_[2]],[x[3],x_[3]]])
0095         ety_ = ety.subs([[y[0],y_[0]],[y[1],y_[1]],[y[2],y_[2]],[y[3],y_[3]]])
0096
0097         xco = Poly(etx,t).all_coeffs()
0098         yco = Poly(ety,t).all_coeffs()
0099         xco_ = Poly(etx_,t).all_coeffs()
0100         yco_ = Poly(ety_,t).all_coeffs()
0101
0102         log.info("xco %r %r %r  " % ( xco, xco_, x_ ))
0103         log.info("yco %r %r %r  " % ( yco, yco_, y_ ))
0104
0105         self.x_ = x_
0106         self.y_ = y_
0107         self.t = t
0108         self.etx_ = etx_
0109         self.ety_ = ety_
0110
0111         z,z1,z2 = symbols("z,z1,z2")
0112
0113         z1_ = x_[0]
0114         z2_ = x_[-1]
0115
0116         self.z = z
0117         self.z1_ = z1_
0118         self.z2_ = z2_
0119
0120         etz_ = ety_.subs(t, (z-z1)/(z2-z1) ).subs( [[z1,z1_],[z2,z2_]] )
0121         zco_ = Poly(etz_,z).all_coeffs()
0122
0123         self.xco = xco
0124         self.yco = yco
0125         self.xco_ = xco_
0126         self.yco_ = yco_
0127
0128         self.etz_ = etz_
0129         self.zco_ = zco_
0130
0131
0132     def zval(self, zz):
0133         return map(lambda _:self.etz_.subs(self.z,_), zz )
0134
0135     def xval(self, tt):
0136         return list(map(lambda _:self.etx_.subs(self.t,_), tt ))
0137     def yval(self, tt):
0138         return list(map(lambda _:self.ety_.subs(self.t,_), tt ))
0139
0140     def plot(self, ax, mx=50, my=50):
0141
0142         tt = np.linspace(0,1,100)
0143
0144         ax.set_xlim(bz.x_[0]-mx,bz.x_[3]+mx)
0145         ax.set_ylim(   0, max(bz.y_)+my)
0146         ax.plot(bz.xval(tt), bz.yval(tt), 'k')
0147
0148         mk = "ob og og or".split()
0149         for i in range(4):
0150             ax.plot( bz.x_[i], bz.y_[i], mk[i] )
0151         pass
0152
0153         ax.plot((bz.x_[3],bz.x_[2]),(bz.y_[3],bz.y_[2]), 'r-')
0154         ax.plot((bz.x_[0],bz.x_[1]),(bz.y_[0],bz.y_[1]), 'b-')
0155
0156
0157
0158 if __name__ == '__main__':
0159
0160
0161     plt.ion()
0162     plt.close()
0163
0164     xy = [[-100,30],[-50,80],[50,30],[100,100]]
0165
0166     bz = Bezier(xy)
0167     fig = plt.figure()
0168     ax = fig.add_subplot(111)
0169     bz.plot(ax)
0170     plt.show()
0171
0172     print(Bezier.ZCoeff(xy))
0173
0174
0175
0176
0177