Mercurial > repos > imgteam > projective_transformation_points

comparison projective_transformation_points.py @ 2:0d2707c82d29 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/BMCV/galaxy-image-analysis/tools/projective_transformation_points/ commit f1298dca5a7f5be3acbb5e3d80c98b1cd6d2795b"
author imgteam
date Fri, 08 May 2020 05:21:51 -0400
parents f1744c5654b9
children aaac58d83043
comparison
equal deleted inserted replaced
1:f1744c5654b9 2:0d2707c82d29
1 from skimage.transform import ProjectiveTransform 1 from skimage.transform import ProjectiveTransform
2 from scipy.ndimage import map_coordinates
2 import numpy as np 3 import numpy as np
3 import pandas as pd 4 import pandas as pd
4 import argparse 5 import argparse
6
7
8 def _stackcopy(a, b):
9 if a.ndim == 3:
10 a[:] = b[:, :, np.newaxis]
11 else:
12 a[:] = b
13
14
15 def warp_img_coords_batch(coord_map, shape, dtype=np.float64, batch_size=1000000):
16 rows, cols = shape[0], shape[1]
17 coords_shape = [len(shape), rows, cols]
18 if len(shape) == 3:
19 coords_shape.append(shape[2])
20 coords = np.empty(coords_shape, dtype=dtype)
21
22 tf_coords = np.indices((cols, rows), dtype=dtype).reshape(2, -1).T
23
24 for i in range(0, (tf_coords.shape[0]//batch_size+1)):
25 tf_coords[batch_size*i:batch_size*(i+1)] = coord_map(tf_coords[batch_size*i:batch_size*(i+1)])
26 tf_coords = tf_coords.T.reshape((-1, cols, rows)).swapaxes(1, 2)
27
28 _stackcopy(coords[1, ...], tf_coords[0, ...])
29 _stackcopy(coords[0, ...], tf_coords[1, ...])
30 if len(shape) == 3:
31 coords[2, ...] = range(shape[2])
32
33 return coords
5 34
6 35
7 def warp_coords_batch(coord_map, coords, dtype=np.float64, batch_size=1000000): 36 def warp_coords_batch(coord_map, coords, dtype=np.float64, batch_size=1000000):
8 tf_coords = coords.astype(np.float32)[:, ::-1] 37 tf_coords = coords.astype(np.float32)[:, ::-1]
9 38
10 for i in range(0, (tf_coords.shape[0]//batch_size)+1): 39 for i in range(0, (tf_coords.shape[0]//batch_size)+1):
11 tf_coords[batch_size*i:batch_size*(i+1)] = coord_map(tf_coords[batch_size*i:batch_size*(i+1)]) 40 tf_coords[batch_size*i:batch_size*(i+1)] = coord_map(tf_coords[batch_size*i:batch_size*(i+1)])
12 41
13 return np.unique(np.round(tf_coords).astype(coords.dtype),axis=0)[:, ::-1] 42 return tf_coords[:, ::-1]
14 43
15 44
16 def transform(coords, warp_matrix, out): 45 def transform(coords, warp_matrix, out):
17 indices = np.array(pd.read_csv(coords, delimiter="\t")) 46 roi_coords = np.array(pd.read_csv(coords, delimiter="\t"))
18 a_matrix = np.array(pd.read_csv(warp_matrix, delimiter="\t", header=None)) 47 trans_matrix = np.array(pd.read_csv(warp_matrix, delimiter="\t", header=None))
19 48
20 trans = ProjectiveTransform(matrix=a_matrix) 49 tol = 10
21 warped_coords = warp_coords_batch(trans, indices) 50 moving = np.zeros(np.max(roi_coords,axis=0)+tol, dtype=np.int8)
22 51 idx_roi_coords = (roi_coords[:,0]-1) * moving.shape[1] + roi_coords[:,1] - 1
52 moving.flat[idx_roi_coords] = 1
53
54 transP = ProjectiveTransform(matrix=trans_matrix)
55 roi_coords_warped_direct = warp_coords_batch(transP, roi_coords)
56 shape_fixed = np.round(np.max(roi_coords_warped_direct,axis=0)).astype(roi_coords.dtype)+tol
57
58 transI = ProjectiveTransform(matrix=np.linalg.inv(trans_matrix))
59 img_coords_warped = warp_img_coords_batch(transI, shape_fixed)
60
61 moving_warped = map_coordinates(moving, img_coords_warped, mode='constant', cval=0)
62 idx_roi_coords_warped = np.where(moving_warped==1)
63
23 df = pd.DataFrame() 64 df = pd.DataFrame()
24 df['x'] = warped_coords[:,0] 65 df['x'] = idx_roi_coords_warped[0] + 1
25 df['y'] = warped_coords[:,1] 66 df['y'] = idx_roi_coords_warped[1] + 1
26 df.to_csv(out, index = False, sep="\t") 67 df.to_csv(out, index = False, sep="\t")
27 68
28 69
29 if __name__ == "__main__": 70 if __name__ == "__main__":
30 parser = argparse.ArgumentParser(description="Transform coordinates") 71 parser = argparse.ArgumentParser(description="Transform coordinates")