# HG changeset patch
# User imgteam
# Date 1588929711 14400
# Node ID 0d2707c82d29dc14188381f6bb90f78bd29da763
# Parent  f1744c5654b96ae1987738c8b50b5b3f3a3d73db
"planemo upload for repository https://github.com/BMCV/galaxy-image-analysis/tools/projective_transformation_points/ commit f1298dca5a7f5be3acbb5e3d80c98b1cd6d2795b"

diff -r f1744c5654b9 -r 0d2707c82d29 projective_transformation_points.py
--- a/projective_transformation_points.py	Wed Mar 27 14:54:32 2019 -0400
+++ b/projective_transformation_points.py	Fri May 08 05:21:51 2020 -0400
@@ -1,28 +1,69 @@
 from skimage.transform import ProjectiveTransform
+from scipy.ndimage import map_coordinates
 import numpy as np
 import pandas as pd
 import argparse
 
 
+def _stackcopy(a, b):
+    if a.ndim == 3:
+        a[:] = b[:, :, np.newaxis]
+    else:
+        a[:] = b
+
+
+def warp_img_coords_batch(coord_map, shape, dtype=np.float64, batch_size=1000000):
+    rows, cols = shape[0], shape[1]
+    coords_shape = [len(shape), rows, cols]
+    if len(shape) == 3:
+        coords_shape.append(shape[2])
+    coords = np.empty(coords_shape, dtype=dtype)
+
+    tf_coords = np.indices((cols, rows), dtype=dtype).reshape(2, -1).T
+
+    for i in range(0, (tf_coords.shape[0]//batch_size+1)):
+        tf_coords[batch_size*i:batch_size*(i+1)] = coord_map(tf_coords[batch_size*i:batch_size*(i+1)])
+    tf_coords = tf_coords.T.reshape((-1, cols, rows)).swapaxes(1, 2)
+
+    _stackcopy(coords[1, ...], tf_coords[0, ...])
+    _stackcopy(coords[0, ...], tf_coords[1, ...])
+    if len(shape) == 3:
+        coords[2, ...] = range(shape[2])
+
+    return coords
+
+
 def warp_coords_batch(coord_map, coords, dtype=np.float64, batch_size=1000000):
     tf_coords = coords.astype(np.float32)[:, ::-1]
 
     for i in range(0, (tf_coords.shape[0]//batch_size)+1):
         tf_coords[batch_size*i:batch_size*(i+1)] = coord_map(tf_coords[batch_size*i:batch_size*(i+1)])
 
-    return np.unique(np.round(tf_coords).astype(coords.dtype),axis=0)[:, ::-1]
-
+    return tf_coords[:, ::-1]
+ 
 
 def transform(coords, warp_matrix, out):
-    indices = np.array(pd.read_csv(coords, delimiter="\t"))
-    a_matrix = np.array(pd.read_csv(warp_matrix, delimiter="\t", header=None))
+    roi_coords = np.array(pd.read_csv(coords, delimiter="\t"))
+    trans_matrix = np.array(pd.read_csv(warp_matrix, delimiter="\t", header=None))
+    
+    tol = 10
+    moving = np.zeros(np.max(roi_coords,axis=0)+tol, dtype=np.int8)
+    idx_roi_coords = (roi_coords[:,0]-1) * moving.shape[1] + roi_coords[:,1] - 1
+    moving.flat[idx_roi_coords] = 1
     
-    trans = ProjectiveTransform(matrix=a_matrix)
-    warped_coords = warp_coords_batch(trans, indices)
-
+    transP = ProjectiveTransform(matrix=trans_matrix)
+    roi_coords_warped_direct = warp_coords_batch(transP, roi_coords)
+    shape_fixed = np.round(np.max(roi_coords_warped_direct,axis=0)).astype(roi_coords.dtype)+tol
+    
+    transI = ProjectiveTransform(matrix=np.linalg.inv(trans_matrix))
+    img_coords_warped = warp_img_coords_batch(transI, shape_fixed)
+    
+    moving_warped = map_coordinates(moving, img_coords_warped, mode='constant', cval=0)
+    idx_roi_coords_warped = np.where(moving_warped==1)
+    
     df = pd.DataFrame()
-    df['x'] = warped_coords[:,0]
-    df['y'] = warped_coords[:,1]
+    df['x'] = idx_roi_coords_warped[0] + 1
+    df['y'] = idx_roi_coords_warped[1] + 1
     df.to_csv(out, index = False, sep="\t")
 
 
diff -r f1744c5654b9 -r 0d2707c82d29 projective_transformation_points.xml
--- a/projective_transformation_points.xml	Wed Mar 27 14:54:32 2019 -0400
+++ b/projective_transformation_points.xml	Fri May 08 05:21:51 2020 -0400
@@ -1,9 +1,11 @@
-<tool id="ip_projective_transformation_points" name="Projective Transformation" version="0.0.3">
+<tool id="ip_projective_transformation_points" name="Projective Transformation" version="0.0.4">
     <description>of input points</description>
     <requirements>
         <requirement type="package" version="0.14.2">scikit-image</requirement>
+        <requirement type="package" version="1.2.1">scipy</requirement>
         <requirement type="package" version="0.23.4">pandas</requirement>
         <requirement type="package" version="1.15.2">numpy</requirement><!--conflict of 1.15.4 with mkl_fft!?-->
+        <requirement type="package" version="0.15.1">tifffile</requirement>
     </requirements>
     <command>
     <![CDATA[
@@ -19,7 +21,7 @@
     </inputs>
     <outputs>
        <data format="tabular" name="out" />
-    </outputs>
+    </outputs> 
     <tests>
       <test>
         <param name="moving_points" value="points_tsv.tsv" />
diff -r f1744c5654b9 -r 0d2707c82d29 test-data/out.tsv
--- a/test-data/out.tsv	Wed Mar 27 14:54:32 2019 -0400
+++ b/test-data/out.tsv	Fri May 08 05:21:51 2020 -0400
@@ -1,7 +1,7 @@
 x	y
-30	12
 28	14
 28	15
 28	16
 28	17
+30	12
 35	25