# HG changeset patch
# User imgteam
# Date 1658169679 0
# Node ID aee73493bf539f1bbf885f7a5591c8d781a5aade
# Parent  9ccd642e7ae2f929177dff0e91c6bc55a470bcac
planemo upload for repository https://github.com/BMCV/galaxy-image-analysis/tools/landmark_registration/ commit ba383a4f617791d0c84737a179e5b4b66cecc334

diff -r 9ccd642e7ae2 -r aee73493bf53 landmark_registration.py
--- a/landmark_registration.py	Sat Feb 26 17:14:05 2022 +0000
+++ b/landmark_registration.py	Mon Jul 18 18:41:19 2022 +0000
@@ -10,12 +10,39 @@
 
 import numpy as np
 import pandas as pd
+from scipy import spatial
 from scipy.linalg import lstsq
 from skimage.measure import ransac
 from skimage.transform import AffineTransform
 
 
-def landmark_registration(pts_f1, pts_f2, out_f, res_th=None, max_ite=None, delimiter="\t"):
+class pwlTransform(object):
+
+    def __init__(self):
+        self.triangulation = None
+        self.affines = None
+
+    def estimate(self, src, dst):
+        self.triangulation = spatial.Delaunay(src)
+        success = True
+        self.affines = []
+        for tri in self.triangulation.simplices:
+            affine = AffineTransform()
+            success &= affine.estimate(src[tri, :], dst[tri, :])
+            self.affines.append(affine)
+        return success
+
+    def __call__(self, coords):
+        simplex = self.triangulation.find_simplex(coords)
+        simplex[simplex == -1] = 0    # todo: dealing with points outside the triangulation
+        out = np.empty_like(coords, np.float64)
+        for i in range(len(self.triangulation.simplices)):
+            idx = simplex == i
+            out[idx, :] = self.affines[i](coords[idx, :])
+        return out
+
+
+def landmark_registration(pts_f1, pts_f2, out_f, pts_f=None, res_th=None, max_ite=None, delimiter="\t"):
 
     points1 = pd.read_csv(pts_f1, delimiter=delimiter)
     points2 = pd.read_csv(pts_f2, delimiter=delimiter)
@@ -31,6 +58,21 @@
         model_robust, inliers = ransac((src, dst), AffineTransform, min_samples=3, residual_threshold=res_th, max_trials=max_ite)
         pd.DataFrame(model_robust.params).to_csv(out_f, header=None, index=False, sep="\t")
 
+    elif pts_f is not None:
+        pwlt = pwlTransform()
+        pwlt.estimate(src, dst)
+
+        pts_df = pd.read_csv(pts_f, delimiter=delimiter)
+        pts = np.concatenate([np.array(pts_df['x']).reshape([-1, 1]),
+                              np.array(pts_df['y']).reshape([-1, 1])],
+                             axis=-1)
+        pts_pwlt = pwlt(pts)
+
+        df = pd.DataFrame()
+        df['x'] = pts_pwlt[:, 0]
+        df['y'] = pts_pwlt[:, 1]
+        df.to_csv(out_f, index=False, sep="\t", float_format='%.1f')
+
     else:
         A = np.zeros((src.size, 6))
         A[0:src.shape[0], [0, 1]] = src
@@ -47,14 +89,18 @@
 
 
 if __name__ == "__main__":
-    parser = argparse.ArgumentParser(description="Estimate affine transformation matrix based on landmark coordinates")
-    parser.add_argument("fn_pts1", help="Coordinates of SRC landmarks (tsv file)")
-    parser.add_argument("fn_pts2", help="Coordinates of DST landmarks (tsv file)")
-    parser.add_argument("fn_tmat", help="Path the output (transformation matrix)")
+    parser = argparse.ArgumentParser(description="estimates the affine transformation matrix or performs piecewiese affine transformation based on landmarks")
+    parser.add_argument("fn_lmkmov", help="Coordinates of moving landmarks (tsv file)")
+    parser.add_argument("fn_lmkfix", help="Coordinates of fixed landmarks (tsv file)")
+    parser.add_argument("fn_out", help="Path to the output")
+    parser.add_argument("--pwlt", dest="fn_ptsmov", help="Coordinates of points to be transformed (tsv file)")
     parser.add_argument("--res_th", dest="res_th", type=float, help="Maximum distance for a data point to be classified as an inlier")
     parser.add_argument("--max_ite", dest="max_ite", type=int, help="Maximum number of iterations for random sample selection")
     args = parser.parse_args()
 
+    fn_ptsmov = None
+    if args.fn_ptsmov:
+        fn_ptsmov = args.fn_ptsmov
     res_th = None
     if args.res_th:
         res_th = args.res_th
@@ -62,4 +108,4 @@
     if args.max_ite:
         max_ite = args.max_ite
 
-    landmark_registration(args.fn_pts1, args.fn_pts2, args.fn_tmat, res_th=res_th, max_ite=max_ite)
+    landmark_registration(args.fn_lmkmov, args.fn_lmkfix, args.fn_out, pts_f=fn_ptsmov, res_th=res_th, max_ite=max_ite)
diff -r 9ccd642e7ae2 -r aee73493bf53 landmark_registration.xml
--- a/landmark_registration.xml	Sat Feb 26 17:14:05 2022 +0000
+++ b/landmark_registration.xml	Mon Jul 18 18:41:19 2022 +0000
@@ -1,5 +1,5 @@
-<tool id="ip_landmark_registration" name="Landmark Registration" version="0.0.4" profile="20.05">
-    <description>estimates the affine transformation matrix</description>
+<tool id="ip_landmark_registration" name="Landmark Registration" version="0.0.5" profile="20.05">
+    <description>estimates the affine transformation matrix or performs piecewise affine transformation</description>
     <requirements>
         <requirement type="package" version="0.18.1">scikit-image</requirement>
         <requirement type="package" version="1.6.2">scipy</requirement>
@@ -9,58 +9,69 @@
     <command detect_errors="aggressive">
     <![CDATA[ 
         python '$__tool_directory__/landmark_registration.py'
-        '$fn_pts1'
-        '$fn_pts2'
-        '$fn_tmat'
-        #if '$algo_option.algo' == 'ransac':
-            --res_th $algo_option.res_thr
-            --max_ite $algo_option.max_iter
+        '$fn_lmkmov'
+        '$fn_lmkfix'
+        '$fn_out'
+        #if $algo_option.algo == "pwat"
+        --pwlt '$algo_option.fn_ptsmov'
+        #elif $algo_option.algo == "ransac"
+        --res_th $algo_option.res_th
+        --max_ite $algo_option.max_ite
         #end if
-]]></command>
+    ]]>
+    </command>
     <inputs>
-        <param name="fn_pts1" type="data" format="tabular" label="Coordinates of SRC landmarks (tsv file)" />
-        <param name="fn_pts2" type="data" format="tabular" label="Coordinates of DST landmarks (tsv file)" />
+        <param name="fn_lmkmov" type="data" format="tabular" label="Coordinates of moving landmarks (tsv file)" />
+        <param name="fn_lmkfix" type="data" format="tabular" label="Coordinates of fixed landmarks (tsv file)" />
         <conditional name="algo_option">
             <param name="algo" type="select" label="Select the algorithm">
-                <option value="ransac" selected="True">RANSAC</option>
-                <option value="ls">Least Squares</option>
+                <option value="ransac">Affine Transformation (based on RANSAC)</option>
+                <option value="ls">Affine Transformation (based on Least Squares)</option>
+                <option value="pwat" selected="True">Piecewise Affine Transformation</option>
             </param>
             <when value="ransac">
-                <param name="res_thr" type="float" value="2.0" label="Maximum distance for a data point to be classified as an inlier." />
-                <param name="max_iter" type="integer" value="100" label="Maximum number of iterations for random sample selection." />
+                <param name="res_th" type="float" value="2.0" label="Maximum distance for a data point to be classified as an inlier." />
+                <param name="max_ite" type="integer" value="100" label="Maximum number of iterations for random sample selection." />
             </when>
             <when value="ls"></when>
+            <when value="pwat">
+                <param name="fn_ptsmov" type="data" format="tabular" label="Coordinates of points to be transformed (tsv file)" />
+            </when>
         </conditional>
     </inputs>
     <outputs>
-        <data format="tabular" name="fn_tmat" />
+        <data format="tabular" name="fn_out" />
     </outputs>
     <tests>
         <test>
-            <param name="fn_pts1" value="points_moving.tsv"/>
-            <param name="fn_pts2" value="points_fixed.tsv"/>
-            <conditional name="algo_option">
-                <param name="algo" value="ls"/>
-            </conditional>
-            <output name="fn_tmat" value="tmat.tsv" ftype="tabular" compare="diff" lines_diff="6"/>
+            <param name="fn_lmkmov" value="points_moving.tsv"/>
+            <param name="fn_lmkfix" value="points_fixed.tsv"/>
+            <param name="algo" value="ls"/>
+            <output name="fn_out" value="tmat.tsv" ftype="tabular" compare="diff" lines_diff="6"/>
         </test>
         <test>
-            <param name="fn_pts1" value="points_moving.tsv"/>
-            <param name="fn_pts2" value="points_fixed.tsv"/>
-            <conditional name="algo_option">
-                <param name="algo" value="ransac"/>
-                <param name="res_thr" value="2"/>
-                <param name="max_iter" value="100"/>
-            </conditional>
-            <output name="fn_tmat" value="tmat.tsv" ftype="tabular" compare="diff" lines_diff="6"/>
+            <param name="fn_lmkmov" value="points_moving.tsv"/>
+            <param name="fn_lmkfix" value="points_fixed.tsv"/>
+            <param name="algo" value="ransac"/>
+            <param name="res_th" value="2.0"/>
+            <param name="max_ite" value="100"/>
+            <output name="fn_out" value="tmat.tsv" ftype="tabular" compare="diff" lines_diff="6"/>
+        </test>
+        <test>
+            <param name="fn_lmkmov" value="lmk_mov.tsv"/>
+            <param name="fn_lmkfix" value="lmk_fix.tsv"/>
+            <param name="algo" value="pwat"/>
+            <param name="fn_ptsmov" value="points.tsv"/>
+            <output name="fn_out" value="points_pwlt.tsv" ftype="tabular"/>
         </test>
     </tests>
     <help>
     **What it does**
 
-    This tool estimates the transformation matrix between two sets of 2d points.
+    1) Estimation of the affine transformation matrix between two sets of 2d points; 
+    2) Piecewise affine transformation of points based on landmark pairs.
 
-    About the format of landmark coordinates in the input TSV table: Columns with header "x" and "y" are for x- and y-coordinate, respectively.
+    About the format of landmark/point coordinates in the input TSV table: Columns with header "x" and "y" are for x- and y-coordinate, respectively.
     </help>
     <citations>
         <citation type="doi">10.1016/j.jbiotec.2017.07.019</citation>
diff -r 9ccd642e7ae2 -r aee73493bf53 test-data/lmk_fix.tsv
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/lmk_fix.tsv	Mon Jul 18 18:41:19 2022 +0000
@@ -0,0 +1,7 @@
+x	y
+6	1710
+546	474
+909	150
+1443	51
+1983	225
+2352	678
\ No newline at end of file
diff -r 9ccd642e7ae2 -r aee73493bf53 test-data/lmk_mov.tsv
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/lmk_mov.tsv	Mon Jul 18 18:41:19 2022 +0000
@@ -0,0 +1,7 @@
+x	y
+51.253	1661.5
+292.57	536.03
+640.67	232.78
+1106.2	106.78
+1576	192.2
+1969	533.89
\ No newline at end of file
diff -r 9ccd642e7ae2 -r aee73493bf53 test-data/points.tsv
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/points.tsv	Mon Jul 18 18:41:19 2022 +0000
@@ -0,0 +1,20 @@
+x	y
+563.79	386.54
+429.25	602.23
+448.47	681.24
+841.67	405.33
+1088.3	350.67
+767	674.67
+951	506.67
+1140.3	533.33
+1350	568
+1368	612
+1296	576
+1834	406
+1608	318
+1110	1194
+366	1770
+1322	1282
+1912	662
+867	322.67
+764.33	424
\ No newline at end of file
diff -r 9ccd642e7ae2 -r aee73493bf53 test-data/points_pwlt.tsv
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/points_pwlt.tsv	Mon Jul 18 18:41:19 2022 +0000
@@ -0,0 +1,20 @@
+x	y
+832.3	326.6
+676.7	565.3
+677.8	656.5
+1132.4	383.5
+1395.9	347.8
+1022.4	688.3
+1254.2	519.6
+1459.2	575.2
+1676.9	640.1
+1685.3	691.8
+1616.7	642.5
+2238.4	517.6
+1983.1	373.4
+1262.6	1314.6
+318.0	1870.8
+1469.0	1439.7
+2258.7	815.1
+1156.4	286.6
+1049.8	396.6