Still busy, but this will delay you some more
Consider the following points:
>>> a = np.array([[1,1], [2,2], [3,3], [4,4]], dtype='float64')
>>> b = affine(a, angle=-22.5)[0]
>>> c = b + [10,10]
>>>
>>> a
array([[ 1.00000000, 1.00000000],
[ 2.00000000, 2.00000000],
[ 3.00000000, 3.00000000],
[ 4.00000000, 4.00000000]])
>>> b
array([[-1.95984445, -0.81179415],
[-0.65328148, -0.27059805],
[ 0.65328148, 0.27059805],
[ 1.95984445, 0.81179415]])
>>> c
array([[ 8.04015555, 9.18820585],
[ 9.34671852, 9.72940195],
[ 10.65328148, 10.27059805],
[ 11.95984445, 10.81179415]])
This function (although only lightly tested) should translate and rotate the points from one form to another.
def warp(src=None, dest=None):
""" move the src points based on the dest pnt affine parameters
"""
a_slope, a_inter, a_cent = reg_params(src)
b_slope, b_inter, b_cent = reg_params(dest)
d_slope = b_slope - a_slope
results = affine(src, d_slope)
new_pnts = results[0] + b_cent
return new_pnts, results
def reg_params(src=None):
""" determine the correlation properties of a point set
"""
a_cent = cent_pnt(src)
a_t = translate(src, a_cent)
xt = a_t[:,0]
yt = a_t[:,1]
rho = np.corrcoef(x=xt, y=yt, rowvar=1, bias=0, ddof=None)
arr = np.vstack([xt, np.ones(len(xt))]).T
slope, inter = np.linalg.lstsq(arr, yt)[0]
deg_slope = np.degrees(np.arctan2(slope,1.0))
return deg_slope, inter, a_cent
I left in 3 commented lines for examining outputs of the covariance and x, y variance calculations if one wants to examine them. The other line, just shows an alternate method of producing the arr variable.
So if we put this together (with reg_params duplicated), then we have
def cent_pnt(pnts):
"""Return the centre of a point array"""
pnts = np.asarray(pnts)
return np.average(pnts, axis=0).squeeze()
def translate(pnts, cent=None):
"""Translate the points about the origin
:Requires:
:--------
: pnts - 2D sequence of point coordinates (N, 2)
: cent - a list, tuple, ndarray of xc, yc
:
:Returns:
:-------
: p_trans - input points translated about the origin
:
:-------
"""
pnts = np.asarray(pnts)
if cent is None:
p_trans = pnts - cent_pnt(pnts)
else:
p_trans = pnts - np.asarray(cent)
return p_trans
def rot_matrix(rotation=0):
"""Return the rotation matrix given points and rotation angle
:Requires:
:--------
: pnts - 2D sequence of point coordinates (N, 2)
: rotation - angle in degrees
:
:Returns:
:-------
: rot_mat - rotation matrix for 2D transform
:
:Notes:
:-----
: for homogenous coordinates use
: rot_mat = np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]])
:--------
"""
rad = np.deg2rad(rotation)
c = np.cos(rad)
s = np.sin(rad)
rot_mat = np.array([[c, -s], [s, c]])
return rot_mat
def affine(pnts, angle=0, cent=None):
"""Produce an affine transformation of a point set, given an angle of
: rotation and a point to rotate about.
:Requires:
:--------
: pnts - a series of 2D points a list or ndarray (N,2)
: angle - rotation angle, in degrees, relative to the x-axis.
: cent - center of rotation, either provided or determined from the
: center of the point cloud. [x,y], np.array([x,y])
:
:Returns:
: p_rot - the rotated points.
:
:Notes:
:-----
: The functions that are used are: cent_pnt, translate, rot_matrix.
: Once the points are translated, and the rotation matrix is obtained,
: the points are transposed to facilite broadcasting, the dot product
: determined and the resultant transposed back to the input shape.
:
"""
p = np.asarray(pnts)
if cent is None:
cent = cent_pnt(pnts)
p_trans = translate(p, cent=None)
rot_mat = rot_matrix(rotation=angle)
p_rot = np.dot(rot_mat, p_trans.T).T
return p_rot, cent, rot_mat
Notice I did a slice on the warp in this example, because I didn't need the rest of the returned values from the warp. Looks pretty good when you compare the from-to results with the inputs.
>>> warp(a, b)[0]
array([[-1.95984445, -0.81179415],
[-0.65328148, -0.27059805],
[ 0.65328148, 0.27059805],
[ 1.95984445, 0.81179415]])
>>> warp(a, c)[0]
array([[ 8.04015555, 9.18820585],
[ 9.34671852, 9.72940195],
[ 10.65328148, 10.27059805],
[ 11.95984445, 10.81179415]])
>>>
>>>
>>>
>>> warp(b, c)[0]
array([[ 8.04015555, 9.18820585],
[ 9.34671852, 9.72940195],
[ 10.65328148, 10.27059805],
[ 11.95984445, 10.81179415]])
>>> warp(c, a)[0]
array([[ 1.00000000, 1.00000000],
[ 2.00000000, 2.00000000],
[ 3.00000000, 3.00000000],
[ 4.00000000, 4.00000000]])
I hope this gets the idea across Neil. This is somewhat simplified and I will still look for the fuller code but the basics entail determining the nature of the parameters of the input points (center, slope etc), then translating, rotating and translating to new center.