Commit b819164b authored by inhuszar's avatar inhuszar
Browse files

Version commit 2.1.3b1

parent 184d92d7
......@@ -9,3 +9,4 @@
/tirl.egg-info/
*__pycache__
.*
/pdf/
\ No newline at end of file
{% set version = "2.1.2b1" %}
{% set version = "2.1.3b1" %}
package:
name: tirl
version: {{ version }}
......
......@@ -56,7 +56,7 @@ with open("requirements.txt", "r") as fp:
#-------------------------------- TIRL Installer ------------------------------#
setup(name="tirl",
version="2.1.2b1",
version="2.1.3b1",
description="Tensor Image Registration Library",
author="Istvan N. Huszar",
ext_modules=cython_modules,
......
......@@ -10,11 +10,12 @@
"system": "linux",
"loglevel": "debug",
"logfile": null,
"paramlogfile": "/Users/inhuszar/temp/example/stage1/paramlog.log",
"paramlogfile": null,
"verbose": false,
"outputdir": "/Users/inhuszar/temp/example/stage1",
"stages": ["rotation", "rigid", "affine", "nonlinear"],
"warnings": false
"warnings": false,
"isotropic": true
},
"histology": {
"file": "/Users/inhuszar/temp/example/1_histology/histology.tif",
......@@ -35,14 +36,14 @@
"snapshot": true
},
"block": {
"file": "/Users/inhuszar/temp/example/2_tissue_block/tissue_block.tif",
"file": null,
"storage": "mem",
"dtype": "f4",
"resolution": 0.050,
"resolution": 0.05,
"mask": {
"file": "/Users/inhuszar/temp/example/2_tissue_block/tissue_block_mask.tif",
"file": null,
"normalise": true,
"function": null,
"function": "dilated_object_mask",
"automask": {
"thr": 0,
"uthr": 1
......@@ -53,15 +54,15 @@
"snapshot": false
},
"preprocessing": {
"histology": ["match_block_resolution", "histology_preprocessing"],
"histology": ["match_block_resolution", "histology_preprocessing", "pad"],
"block": ["block_preprocessing"]
},
"regparams": {
"init": {
"scale": {
"x0": [1.0, 1.0],
"lb": [0.9, 0.9],
"ub": [1.1, 1.1]
"x0": 1.0,
"lb": 0.9,
"ub": 1.1
},
"rotation": {
"x0": 0.0,
......@@ -76,8 +77,8 @@
},
"affine": {
"x0": [1.0, 0.0, 0.0, 0.0, 1.0, 0.0],
"lb": [0.8, -1.0, -5.0, -1.0, 0.8, -5.0],
"ub": [1.2, 1.0, 5.0, 1.0, 1.2, 5.0]
"lb": [0.95, -0.2, -1.0, -0.2, 0.95, -1.0],
"ub": [1.05, 0.2, 1.0, 0.2, 1.05, 1.0]
},
"nonlinear": {
"x0": 0.0,
......@@ -86,39 +87,39 @@
}
},
"rotsearch": {
"coarse": 90,
"scale": 0.1,
"coarse": 30,
"scale": 0.125,
"visualise": false,
"xtol_rel": 0.01,
"xtol_abs": [0.01, 0.001, 0.001, 0.1, 0.1],
"opt_step": [0.1, 0.01, 0.01, 1.0, 1.0]
"xtol_abs": [0.001, 0.001, 0.001, 0.001],
"opt_step": 0.5
},
"rigid": {
"scaling": [30, 20],
"smoothing": [0, 0],
"scaling": [16, 8, 4],
"smoothing": [0, 0, 0],
"xtol_rel": 0.01,
"xtol_abs": [0.01, 0.001, 0.001, 0.1, 0.1],
"opt_step": [0.1, 0.01, 0.01, 1.0, 1.0],
"xtol_abs": [0.001, 0.001, 0.001, 0.001],
"opt_step": 0.1,
"visualise": false
},
"affine": {
"scaling": [20, 10, 5, 2],
"smoothing": [0, 0, 0, 0],
"scaling": [8, 4],
"smoothing": [0, 0],
"xtol_rel": 0.01,
"xtol_abs": [0.001, 0.1, 0.1, 0.1, 0.001, 0.1],
"opt_step": [0.01, 0.01, 1, 0.01, 0.01, 1],
"xtol_abs": [0.001, 0.001, 0.001, 0.001, 0.001, 0.001],
"opt_step": 0.1,
"visualise": false
},
"nonlinear": {
"scaling": [30, 20, 10, 5, 1],
"smoothing": [0, 0, 0, 0, 0],
"scaling": [16, 8, 4],
"smoothing": [0, 0, 0],
"sigma": 1,
"truncate": 1.5,
"regweight": 0.5,
"regweight": 0.3,
"maxiter": [20, 20, 20, 20, 5],
"xtol_abs": 0.1,
"xtol_rel": 0.01,
"visualise": false
}
}
}
}
\ No newline at end of file
......@@ -10,9 +10,9 @@
"system": "linux",
"loglevel": "debug",
"logfile": null,
"paramlogfile": "/Users/inhuszar/temp/example/stage2/paramlog.log",
"paramlogfile": null,
"verbose": false,
"outputdir": "/Users/inhuszar/temp/example/stage2",
"outputdir": "",
"stages": ["rigid", "affine", "nonlinear"],
"warnings": false
},
......@@ -22,7 +22,7 @@
"dtype": "f4",
"resolution": 0.05,
"mask": {
"file": "/Users/inhuszar/temp/example/2_tissue_block/tissue_block_mask.tif",
"file": null,
"normalise": true,
"function": null,
"automask": {
......@@ -85,7 +85,7 @@
"ub": null
}
},
"sites": "/Users/inhuszar/temp/example/3_brain_slice/sites/sites.txt",
"sites": "",
"jiggle": {
"scale": 0.1,
"xrange": [10.0, 10.0],
......@@ -125,4 +125,4 @@
"visualise": false
}
}
}
}
\ No newline at end of file
......@@ -12,7 +12,7 @@
"paramlogfile": "/Users/inhuszar/temp/example/stage3/paramlogs.log",
"verbose": false,
"outputdir": "/Users/inhuszar/temp/example/stage3",
"stages": [1, 2, 3, 4, 5],
"stages": [1, 2, 3, 4, 5, 3, 4, 5],
"isotropic": true,
"cost": "MIND",
"warnings": false
......@@ -79,20 +79,20 @@
"centre": [2.2, 66.99, -6.3],
"normal": [0, 1, 0],
"offset": 0,
"thickness": 10,
"n_positions": 5,
"range": [0, 0, 0],
"n_orientations": [1, 1, 1]
"thickness": 20,
"n_positions": 11,
"range": [0, 20, 20],
"n_orientations": [1, 5, 5]
},
"iterations": 1,
"n_cpu": -1,
"n_cpu": 1,
"top": 0.3,
"slice_scaling": [4, 2, 1],
"slice_smoothing": [0, 0, 0],
"volume_scaling": [4, 2, 1],
"volume_smoothing": [0, 0, 0],
"slice_scaling": [2, 1],
"slice_smoothing": [0, 0],
"volume_scaling": [2, 1],
"volume_smoothing": [0, 0],
"constrained": true,
"try_unconstrained": true,
"try_unconstrained": false,
"opt_step": 0.1,
"stage_1a": {
"scale2d_lower_delta": 0.05,
......@@ -122,14 +122,14 @@
"source_mask": false
},
"visualise": false,
"slice_scaling": [2, 1],
"volume_scaling": [2, 1],
"slice_smoothing": [1, 0],
"volume_smoothing": [1, 0],
"slice_scaling": [1],
"volume_scaling": [1],
"slice_smoothing": [0],
"volume_smoothing": [0],
"lower_delta": [0.1, 0.1, 0.1, 1.0, 0.1, 0.1, 0.1, 1.0, 0.1, 0.1, 0.1, 1.0],
"upper_delta": [0.1, 0.1, 0.1, 1.0, 0.1, 0.1, 0.1, 1.0, 0.1, 0.1, 0.1, 1.0],
"xtol_rel": 0.01,
"opt_step": 0.5
"opt_step": 0.1
},
"stage_3": {
"export": {
......@@ -154,11 +154,14 @@
},
"optsize": 4,
"vectorder": "xy",
"slice_smoothing": [0],
"volume_smoothing": [0],
"slice_scaling": [2, 1],
"slice_smoothing": [0, 0],
"volume_scaling": [1, 1],
"volume_smoothing": [0, 0],
"lower_dxy": 5.0,
"upper_dxy": 5.0,
"regweight": 0,
"model": "multiquadric",
"opt_step": 0.1,
"xtol_abs": 0.01
},
......@@ -186,13 +189,16 @@
},
"optsize": 4,
"vectorder": "xyz",
"slice_smoothing": [0],
"volume_smoothing": [0],
"lower_dxy": 3.0,
"lower_dz": 3.0,
"upper_dxy": 3.0,
"upper_dz": 3.0,
"slice_scaling": [2, 1],
"slice_smoothing": [0, 0],
"volume_scaling": [1, 1],
"volume_smoothing": [0, 0],
"lower_dxy": 5.0,
"lower_dz": 5.0,
"upper_dxy": 5.0,
"upper_dz": 5.0,
"regweight": 0,
"model": "multiquadric",
"opt_step": 0.3,
"xtol_abs": 0.01
},
......@@ -215,4 +221,4 @@
}
}
}
}
}
\ No newline at end of file
......@@ -229,7 +229,7 @@ __author__ = "Istvan N. Huszar"
__copyright__ = "Copyright (C) 2018-2020 University of Oxford"
__credits__ = ["Mark Jenkinson"]
__license__ = "FSL"
__version__ = "2.1.2b1"
__version__ = "2.1.3b1"
__maintainer__ = "Istvan N Huszar"
__status__ = "active"
......@@ -331,5 +331,5 @@ def load(fname):
:rtype: TIRLObject
"""
from tirl import tirlobject
return tirlobject.TIRLObject.load(fname)
from tirl.tirlobject import TIRLObject
return TIRLObject.load(fname)
......@@ -24,6 +24,8 @@ cimport numpy as np
from scipy.spatial import cKDTree
from libc.stdlib cimport malloc, free
from tirl.settings import CPU_CORES as n_cpu
# 2D functions
cdef extern from "finv2d_src.h":
......@@ -118,10 +120,12 @@ def jacobian2d(field, locations, pretx=None, vectorder=(0, 1), maxiter=1000):
ic[point] = &c_input_coordinates[point, 0]
# Find nearest grid point to each location
tree = cKDTree(fwd_grid_coordinates, leafsize=10000000)
tree = cKDTree(fwd_grid_coordinates, leafsize=16, balanced_tree=False,
compact_nodes=False)
cdef np.ndarray[int, ndim=1, mode="c"] nearest_neighbour = \
np.ascontiguousarray(
tree.query(locations, k=1)[1].ravel(), dtype=ctypes.c_int)
tree.query(locations, k=1, eps=0.1, p=1, n_jobs=n_cpu)[1].ravel(),
dtype=ctypes.c_int)
cdef int* nn = <int *>malloc(n_points * sizeof(int))
nn = &nearest_neighbour[0]
......@@ -211,8 +215,10 @@ def invert(targetlocations, field, gridlocations, vectorder=None, output=None,
p_mapped_gridlocations[point] = &c_mapped_gridlocations[point, 0]
# Find the nearest gridpoint for every new location using a k-D tree
tree = cKDTree(mapped_gridlocations, leafsize=16)
nearest = tree.query(targetlocations, k=1)[1].ravel()
tree = cKDTree(mapped_gridlocations, leafsize=16, balanced_tree=False,
compact_nodes=False)
nearest = tree.query(
targetlocations, k=1, eps=0.1, p=1, n_jobs=n_cpu)[1].ravel()
cdef np.ndarray[int, ndim=1, mode="c"] c_nearest = \
np.ascontiguousarray(nearest, dtype=ctypes.c_int)
......@@ -310,8 +316,10 @@ def local_affine(field, newlocations, gridlocations, output=None,
p_gridlocations[point] = &c_gridlocations[point, 0]
# Find the nearest gridpoint for every new location using a k-D tree
tree = cKDTree(gridlocations, leafsize=10000000)
nearest = tree.query(newlocations, k=1)[1].ravel()
tree = cKDTree(gridlocations, leafsize=16, balanced_tree=False,
compact_nodes=False)
nearest = tree.query(
newlocations, k=1, eps=0.1, p=1, n_jobs=n_cpu)[1].ravel()
cdef np.ndarray[int, ndim=1, mode="c"] c_nearest = \
np.ascontiguousarray(nearest, dtype=ctypes.c_int)
......@@ -397,10 +405,12 @@ def jacobian3d(field, locations, pretx=None, vectorder=(0, 1, 2), maxiter=1000):
ic[point] = &c_input_coordinates[point, 0]
# Find nearest grid point to each location
tree = cKDTree(fwd_grid_coordinates, leafsize=10000000)
tree = cKDTree(fwd_grid_coordinates, leafsize=16, balanced_tree=False,
compact_nodes=False)
cdef np.ndarray[int, ndim=1, mode="c"] nearest_neighbour = \
np.ascontiguousarray(
tree.query(locations, k=1)[1].ravel(), dtype=ctypes.c_int)
tree.query(locations, k=1, eps=0.1, p=1, n_jobs=n_cpu)[1].ravel(),
dtype=ctypes.c_int)
cdef int* nn = <int *>malloc(n_points * sizeof(int))
nn = &nearest_neighbour[0]
......
......@@ -34,20 +34,20 @@ double Dot2D(double v[], double u[])
return Result;
}
double* Difference2D(double v[], double u[]) {
double* Result = calloc(sizeof(double), 2);
int Difference2D(double v[], double u[], double * Result) {
// double* Result = calloc(sizeof(double), 2);
for (int i = 0; i < 2; i++) {
Result[i] = v[i] - u[i];
}
return Result;
return 0;
}
int* DifferenceInt2D(int v[], int u[]) {
int* Result = calloc(sizeof(int), 2);
int DifferenceInt2D(int v[], int u[], int * Result) {
// int* Result = calloc(sizeof(int), 2);
for (int i = 0; i < 2; i++) {
Result[i] = v[i] - u[i];
}
return Result;
return 0;
}
int Unravel2DIndex(int Index, int* Shape, int* MultiIndex) {
......@@ -93,18 +93,18 @@ int IsInTriangle(double* Point, double** Simplex) {
Vertex = (FirstVertex + Run) % 3;
Vertex2 = (Vertex + 1) % 3;
Vertex3 = (Vertex + 2) % 3;
double* Edge = Difference2D(Simplex[Vertex2], Simplex[Vertex3]);
double Edge[2];
Difference2D(Simplex[Vertex2], Simplex[Vertex3], Edge);
double Normal[2] = {-Edge[1], Edge[0]};
// Vector from index vertex to point of interest
double* PDist = Difference2D(Point, Simplex[Vertex2]);
double PDist[2];
Difference2D(Point, Simplex[Vertex2], PDist);
// Vector from index vertex to secondary vertex
double* VDist = Difference2D(Simplex[Vertex], Simplex[Vertex2]);
double VDist[2];
Difference2D(Simplex[Vertex], Simplex[Vertex2], VDist);
PointInTriangle = ((Dot2D(PDist, Normal) * Dot2D(VDist, Normal)) >= 0) ? -1 : Vertex;
if (PointInTriangle >= 0)
break;
free(Edge);
free(PDist);
free(VDist);
}
return PointInTriangle;
......@@ -115,7 +115,8 @@ int StepTriangle(int** Simplex, int Vertex, int* Shape) {
// Find mirror axis: diagonal or gridline?
int Vertex2 = (Vertex + 1) % 3;
int Vertex3 = (Vertex + 2) % 3;
int* Edge = DifferenceInt2D(Simplex[Vertex2], Simplex[Vertex3]);
int Edge[2];
DifferenceInt2D(Simplex[Vertex2], Simplex[Vertex3], Edge);
// Mirror simplex around diagonal line
if ((Edge[0] != 0) && (Edge[1] != 0)) {
......
......@@ -22,8 +22,8 @@
// Basic math
double Dot2D(double v[], double u[]);
double* Difference2D(double v[], double u[]);
int* DifferenceInt2D(int v[], int u[]);
int Difference2D(double v[], double u[], double* Result);
int DifferenceInt2D(int v[], int u[], int* Result);
int Unravel2DIndex(int Index, int* Shape, int* MultiIndex);
int Ravel2DIndex(int* MultiIndex, int* Shape);
lapack_int InvertMat2D(double* A, int N);
......
......@@ -33,20 +33,18 @@ double Dot3D(double v[], double u[]) {
return Result;
}
double* Cross(double v[], double u[]) {
double * Result = calloc(3, sizeof(double));
int Cross(double v[], double u[], double* Result) {
Result[0] = v[1]*u[2] - v[2]*u[1];
Result[1] = v[2]*u[0] - v[0]*u[2];
Result[2] = v[0]*u[1] - v[1]*u[0];
return Result;
return 0;
}
double* Difference3D(double v[], double u[]) {
double* Result = calloc(3, sizeof(double));
int Difference3D(double v[], double u[], double* Result) {
for (int i = 0; i < 3; i++) {
Result[i] = v[i] - u[i];
}
return Result;
return 0;
}
int Unravel3DIndex(int Index, int* Shape, int* MultiIndex) {
......@@ -101,15 +99,18 @@ int IsInTetrahedron(double* Point, double** Simplex) {
Vertex3 = (Vertex + 2) % 4;
Vertex4 = (Vertex + 3) % 4;
// Find surface normal for the side that is opposite of Vertex1
double* Edge1 = Difference3D(Simplex[Vertex2], Simplex[Vertex3]);
double* Edge2 = Difference3D(Simplex[Vertex2], Simplex[Vertex4]);
double* Normal = Cross(Edge1, Edge2);
double Edge1[3], Edge2[3];
Difference3D(Simplex[Vertex2], Simplex[Vertex3], Edge1);
Difference3D(Simplex[Vertex2], Simplex[Vertex4], Edge2);
double Normal[3];
Cross(Edge1, Edge2, Normal);
// Vector from index vertex to point of interest
double* PDist = Difference3D(Point, Simplex[Vertex2]);
double PDist[3];
Difference3D(Point, Simplex[Vertex2], PDist);
// Vector from index vertex to secondary vertex
double* VDist = Difference3D(Simplex[Vertex], Simplex[Vertex2]);
double VDist[3];
Difference3D(Simplex[Vertex], Simplex[Vertex2], VDist);
PointInTetrahedron = ((Dot3D(PDist, Normal) * Dot3D(VDist, Normal)) >= 0) ? -1 : Vertex;
free(Edge1); free(Edge2); free(Normal); free(PDist); free(VDist);
if (PointInTetrahedron >= 0)
break;
}
......@@ -143,6 +144,9 @@ int InitialiseTetrahedron(int ** Simplex, int* Shape) {
else
{ _i4 = _i1; _j4 = _j1; _k4 = _k1 + 1; }
// Multiplex the fourth vertex if the input is practically 2D
if (_ks == 1) _k4 = _k1;
// Update the vertices of the tetrahedron
Simplex[1][0] = _i2; Simplex[1][1] = _j2; Simplex[1][2] = _k2;
Simplex[2][0] = _i3; Simplex[2][1] = _j3; Simplex[2][2] = _k3;
......
......@@ -22,9 +22,9 @@
// Basic math
double Dot3D(double v[], double u[]);
double* Cross(double v[], double u[]);
double* Difference3D(double v[], double u[]);
int* DifferenceInt3D(int v[], int u[]);
int Cross(double v[], double u[], double* Result);
int Difference3D(double v[], double u[], double* Result);
int DifferenceInt3D(int v[], int u[], int* Result);
int Unravel3DIndex(int Index, int* Shape, int* MultiIndex);
int _lin(int* Shape, int i, int j, int k);
int Ravel3DIndex(int* MultiIndex, int* Shape);
......
......@@ -180,6 +180,7 @@ class Cost(TIRLObject):
# Set target image directly if no filter was specified
if target_filter is None:
self._target = target
self._trg = target
# Attempt to apply target filter
elif callable(target_filter):
......@@ -188,6 +189,8 @@ class Cost(TIRLObject):
except Exception as exc:
raise RuntimeError("Target image filtering failed.") \
from exc
else:
self._trg = target
else:
raise TypeError(f"Target image filter is not callable: "
f"{target_filter.__class__.__name__}")
......@@ -212,6 +215,7 @@ class Cost(TIRLObject):
# Set source image directly if no prefilter was specified
if source_prefilter is None:
self._source = source
self._src = source
# Attempt to apply pre-filter
elif callable(source_prefilter):
......@@ -220,6 +224,8 @@ class Cost(TIRLObject):
except Exception as exc:
raise RuntimeError("Source image prefiltering failed.") \
from exc
else:
self._src = source
else:
raise TypeError(f"Source image prefilter is not callable: "
f"{source_prefilter.__class__.__name__}")
......
......@@ -299,7 +299,8 @@ class SpatialOperator(Operator):
dimensions.
"""
chunkshape = np.asarray(tf.tshape + (1 + 2 * self.radius,) * tf.vdim)
chunkshape = np.asarray(
tf.tshape[:tf.tdim] + (1 + 2 * self.radius,) * tf.vdim)
if np.product(chunkshape) > chunksize:
raise MemoryError(
"Not enough memory to allocate the minimal chunk.")
......
# from protocols import tirl_affine2d
# from protocols import tirl_affine3d
# from protocols import tirl_rigid2d
# from protocols import tirl_rigid3d
# from protocols import tirl_rigid2d3d
# from protocols import tirl_affine2d3d
# from protocols import tirl_nonlin2d
......@@ -29,4 +29,6 @@ The package includes the following modules:
"""
__version__ = "2020.7.1"
from tirl.scripts.mnd import general
\ No newline at end of file
......@@ -29,7 +29,7 @@ photographs by the accompanying "find_sites" script.
"""
__version__ = "2020.7.0"