In [1]:
import logging


from skimage.transform import AffineTransform
import numpy as np
from numpy.linalg import inv


def _pix2unit3(x, offset, fov, pixel_size, cam_rotation=None, im_flip=None):
    """
    transform a point from pixel coordinates to NIS stage coordinates,
    taking into account offsets, fov, camera rotation or image flipping
    
    Parameters
    ----------
    x: 4-tuple
        point to transform, in pixels  
    offset: array-like
        extra offset to add to transformed bounding boxes (in units)
        (center of image or center of first tile in large image )
    fov: array-like
        field-of-view size (in units) 
    pixel_size: scalar
        pixel size in units
    cam_rotation: 2x2 mat, optional
        camera rotation matrix as provided by NIS
    im_flip: array-like
        array of 1,-1 indicating whether to flip coordinates in a dimension or not
    
    Returns
    -------
    x_tr: array-like
        transformed point, in units
    """
    
    logger = logging.getLogger(__name__)
    
    # default: no camera rotation
    if cam_rotation is None:
        cam_rotation = np.array([[1,0], [0,1]], dtype=float)
    
    # augmented rotation matrix and inverse
    cam_rot_tr = AffineTransform(np.array([[cam_rotation[0,0], cam_rotation[0,1],  0],
                                           [cam_rotation[1,0], cam_rotation[1,1],  0],
                                           [0,                                 0,  1]])
                                 )
    cam_rot_tr_i = AffineTransform(inv(cam_rot_tr.params))
    
    # default image flip: along y
    im_flip_t = AffineTransform(scale=[1,-1] if im_flip is None else im_flip)
    
    x = np.array(x, dtype=float)
    x_tr = (im_flip_t+cam_rot_tr_i)(x * pixel_size - np.array(fov, dtype=float)/2)

    res = np.squeeze(np.array(offset, dtype=float) + x_tr)
    logger.debug('transformed point {} (pixels) to {} (units)'.format(x, res))
    return res
    
    
def bbox_pix2unit3(bbox, offset, fov, pixel_size, cam_rotation=None, im_flip=None):
    """
    transform a bounding box from pixel coordinates to NIS stage coordinates,
    taking into account offsets, fov, camera rotation or image flipping
    
    Parameters
    ----------
    bbox: 4-tuple
        ymin, xmin, ymax, xmax (as output by skimages regionprops, in pixels)  
    offset: array-like
        extra offset to add to transformed bounding boxes (in units)
        (center of image or center of first tile in large image )
    fov: array-like
        field-of-view size (in units) 
    pixel_size: scalar
        pixel size in units
    cam_rotation: 2x2 mat, optional
        camera rotation matrix as provided by NIS
    im_flip: array-like
        array of 1,-1 indicating whether to flip coordinates in a dimension or not
    
    Returns
    -------
    bbox_tr: 4-tuple
        transformed bounding box (ymin, xmin, ymax, xmax - in units)
    """
    
    logger = logging.getLogger(__name__)
      
    # transform bbox
    (ymin, xmin, ymax, xmax) = bbox    
    bbox_tr = np.apply_along_axis(lambda x: _pix2unit3(x, offset, fov, pixel_size, cam_rotation, im_flip),
                                  1, 
                                  np.array([[xmin, ymin],
                                            [xmin, ymax],
                                            [xmax, ymin],
                                            [xmax, ymax]], dtype=float)
                                  )
    
    # get new min max
    min_ = np.apply_along_axis(np.min, 0, bbox_tr)
    max_ = np.apply_along_axis(np.max, 0, bbox_tr)
    
    logger.debug('new min: {}, new max: {}'.format(min_, max_))
    
    # NB: we reverse here to preserve original ymin, xmin, ymax, xmax - order
    bbox_tr_arr = np.array([list(reversed(list(min_))), list(reversed(list(max_)))], dtype=float)
    res = bbox_tr_arr.ravel()
    
    logger.debug('bbox: {}, toUnit: {}'.format(bbox, res))
    return tuple(list(res))

In [40]:
###############
# some tests
###############

logging.basicConfig(format='%(asctime)s - %(levelname)s in %(funcName)s: %(message)s', level=logging.DEBUG)

"""

cam_rot = AffineTransform(np.array([[-0.99995, -0.0096,  0],
                                    [0.0096,   -0.99995, 0],
                                    [0,        0,        1]]))
cam_rot_i = AffineTransform(inv(cam_rot.params))

im_flip = AffineTransform(scale=[1,-1])


offset = np.array([5000, 5000])
"""
cam_pixs = np.array([512, 512], dtype=float)
psz = 0.8
a = 160 / 360. * 2. * np.pi
flip = -1
r = np.pi
offset = np.array([-500,-500])

img_rot = AffineTransform(np.array([[np.cos(r), -np.sin(r),  0],
                                    [np.sin(r),   np.cos(r), 0],
                                    [0,        0,        1]]))

cam_rot = AffineTransform(np.array([[np.cos(a), -np.sin(a),  0],
                                    [np.sin(a),   np.cos(a), 0],
                                    [0,        0,        1]]))
im_flip = AffineTransform(scale=[1,flip])



for x,y in np.nditer(np.meshgrid([0,1024], [0,1024])):
    a = np.array([x,y]) 
    print('transforming point: {}'.format(a))
    a1 = (im_flip + img_rot).inverse(a-cam_pixs/2)
    print(a1)
    au = a1 * psz
    print(au)
    af = au#-cam_pixs*psz/2
    print(af)
    a2 = cam_rot.inverse(af)
    print(a2)
    a3 = a2 + offset
    print(a3)


#for a in np.nditer(np.meshgrid([-1,1], [-1,1])):
#    print((cam_rot+im_flip+AffineTransform(scale=[1/psz, 1/psz]))( a * cam_pixs*psz/2 + cam_pixs*psz/2))
    
#nv(im_flip.params)

#bbox = [a*cam_pixs for a in np.nditer(np.meshgrid([0,1], [0,1]))]

#bbox_pix2unit3((0,0,512,512), [0,0], [512,512], 0.8, )

#cr = np.array([[-0.99995, -0.0096],[0.0096,   -0.99995]])
#_pix2unit3([12000,300], [0,0], cam_pixs*psz, psz, cr)

#cam_rot_i.params


transforming point: [0 0]
[[ 256. -256.]]
[[ 204.8 -204.8]]
[[ 204.8 -204.8]]
[[-262.49477409  122.40332338]]
[[-762.49477409 -377.59667662]]
transforming point: [1024    0]
[[-768. -256.]]
[[-614.4 -204.8]]
[[-614.4 -204.8]]
[[507.30142086 402.5862248 ]]
[[  7.30142086 -97.4137752 ]]
transforming point: [   0 1024]
[[256. 768.]]
[[204.8 614.4]]
[[204.8 614.4]]
[[  17.68812732 -647.39287156]]
[[ -482.31187268 -1147.39287156]]
transforming point: [1024 1024]
[[-768.  768.]]
[[-614.4  614.4]]
[[-614.4  614.4]]
[[ 787.48432227 -367.20997015]]
[[ 287.48432227 -867.20997015]]