Python - 3D affine transformation

Question

I am working with two similar shaped, but not yet identical volumes with a volume grid inside. I am aiming to match my first volume (green) to my second volume (red). Both have a ConvexHull (http://scipy.github.io/devdocs/generated/scipy.spatial.ConvexHull.html) with inner vertices. I created multiple marker (see Figure 1) for both volumes to calculate a transformation matrix(https://community.esri.com/thread/183601-numpy-linalglstsq-coordinate-translations).

The Volume Grid Data structure of my original volume is:

array([[ 0.025, -0.055, -0.03 ],
       [-0.01 , -0.05 , -0.03 ],
       [-0.005, -0.05 , -0.03 ],
       ..., 
       [-0.01 , -0.03 ,  0.1  ],
       [-0.01 , -0.025,  0.1  ],
       [-0.015, -0.02 ,  0.1  ]])

 ,     with the shape of `(12163, 3)`

The Volume Grid Data structure of my original volume is:

array([[ 0.   , -0.055, -0.065],
       [ 0.005, -0.055, -0.065],
       [-0.005, -0.05 , -0.065],
       ..., 
       [-0.005, -0.02 ,  0.08 ],
       [ 0.   , -0.02 ,  0.08 ],
       [ 0.005, -0.02 ,  0.08 ]])

 ,     with the shape of `(14629, 3)`

The coordinates of the original markers which should be transformed are:

array([[-0.00307161, -0.01828496,  0.03521746],
       [-0.065     , -0.01828496,  0.03521746],
       [ 0.06      , -0.01828496,  0.03521746],
       [-0.00307161, -0.01828496,  0.1       ],
       [-0.00307161,  0.075     ,  0.03521746],
       [-0.00307161, -0.01828496, -0.03      ]])

The template marker are:

array([[ 0.00038417, -0.02389603,  0.00802208],
       [-0.07      , -0.02389603,  0.00802208],
       [ 0.07      , -0.02389603,  0.00802208],
       [ 0.00038417, -0.02389603,  0.08      ],
       [ 0.00038417,  0.07      ,  0.00802208],
       [ 0.00038417, -0.02389603, -0.065     ]])

I take the coordinate points of my markers to calculate the transformation matrix like:

print 'Calculating the transformation matrix..\n'

n = orig_marker.shape[0]
pad = lambda x: np.hstack([x, np.ones((x.shape[0], 1))])
unpad = lambda x: x[:,:-1]
trans_mat, res, rank, s = np.linalg.lstsq(pad(orig_marker), pad(temp_marker))


transform = lambda x: unpad(np.dot(pad(x), trans_mat))
trans_mat[np.abs(trans_mat) < 1e-10] = 0  # set really small values to zero
print 'trans matrix is:\n', trans_mat
trans_mat_inv = np.linalg.inv(trans_mat)

Out [1]:  trans matrix is [[  3.29770822e-02   1.06840729e-02   1.71325156e-03   0.00000000e+00]
     [ -7.56419706e-03   9.51696607e-03   3.51349962e-02   0.00000000e+00]
     [  5.32353680e-03   2.91946064e-01   8.44071139e-01   0.00000000e+00]
     [  1.96037928e-04  -3.51253282e-02  -3.05335725e-02   1.00000000e+00]]

After that I am applying my transformation matrix to my volume grid points by:

# apply rotation and scale
transformed_points = np.dot(orig_points, trans_mat[:3, :3].T)
# apply translation
transformed_points += trans_mat[:3, 3]
x_t, y_t, z_t = transformed_points.T

, where orig_points and temp_points are the volume grids of my volumes x_t, y_t, z_t are the coordinates of my transformed volume grids.

Since I apply rotation, scaling and translation my volume grids should match. unfortunately my volume grid still looks like Figure 2:

I am almost sure my approach is correct. I think the mistake might be in the calculation of the transformation matrix.

Can anyone see what went wrong or where I've made a mistake?

With my own self made translation the result came to look like following:

Since the result isn't accurate I would prefer a correct calculation of my transformation matrix.