Bugfix in OptNL, added isotropic scaling to stage 1.

requirements.txt

@@ -12,7 +12,7 @@ nibabel
 nlopt
 numba
 numpy
-opencv
+opencv-python
 openslide-wrapper
 pillow
 psutil

src/tirl/optimisation/optnl.py

@@ -223,7 +223,7 @@ class OptNL(Optimiser):
             self.log(self.revert_parameters(x))
             return self.measure_cost_normalised(x, grad)
         else:
-            self.log(self.revert_parameters(x))
+            self.log(x)
             return self.measure_cost_absolute(x, grad)
 
     def measure_cost_absolute(self, x, grad=None):

src/tirl/scripts/mnd/histology_to_block.py

@@ -83,7 +83,7 @@ from tirl.optimisation.optgroup import OptimisationGroup
 from tirl.optimisation.gnoptdiff import GNOptimiserDiffusion
 from tirl.regularisers.diffusion import DiffusionRegulariser
 
-from tirl.transformations.linear.scale import TxScale
+from tirl.transformations.linear.scale import TxScale, TxIsoScale
 from tirl.transformations.linear.rotation import TxRotation2D
 from tirl.transformations.linear.translation import TxTranslation
 from tirl.transformations.linear.affine import TxAffine
@@ -214,15 +214,22 @@ def labkmeans(histo, **kwargs):
     return mask
 
 
-def initialise_transformations(fixed, q):
+def initialise_transformations(fixed, p):
     """
     Create transformation chain that will be optimised.
 
     """
+    q = p.regparams
+
     # Scale
     lb = np.asarray(q.init.scale.lb)
     ub = np.asarray(q.init.scale.ub)
-    tx_scale = TxScale(*q.init.scale.x0, bounds=(lb, ub), name="scale")
+    if p.general.isotropic:
+        bounds = (float(lb), float(ub))
+        tx_scale = TxIsoScale(
+            float(q.init.scale.x0), dim=2, bounds=bounds, name="scale")
+    else:
+        tx_scale = TxScale(*q.init.scale.x0, bounds=(lb, ub), name="scale")
 
     # Rotation
     if str(q.init.rotation.mode).lower() == "deg":
@@ -295,7 +302,7 @@ def register(fixed, moving, cnf):
     # Create transformation chain (does not change the fixed image)
     # rotation -> scale -> translation -> affine -> nonlinear
     logger.info("Initialising transformation chain...")
-    chain = initialise_transformations(fixed, q)
+    chain = initialise_transformations(fixed, p)
     logger.info("Transformation chain has been initialised.")
 
     # Generate output: initial alignment
@@ -402,7 +409,8 @@ def rotation_search2d(fixed, moving, cnf):
         tx_rotation.parameters.parameters[0] = angle
         tx_rotation.parameters.set_lower_bounds(angle - step / 2)
         tx_rotation.parameters.set_upper_bounds(angle + step / 2)
-        cost = CostMI(moving, fixed, normalise=True, bins=32)()
+        cost = CostMIND(moving, fixed, normalise=True, kernel=MK_FULL)()
+        # cost = CostMI(moving, fixed, normalise=True, bins=32)()
         logger.debug(f"{degrees(angle)} deg: {cost}")
         costvals.append([cost, angle])
     else:
@@ -421,27 +429,28 @@ def rotation_search2d(fixed, moving, cnf):
     tx_translation = fixed.domain.chain["translation"]
     og = OptimisationGroup(tx_rotation, tx_scale, tx_translation)
     scale_orig = tx_scale.parameters.parameters.copy()
-    translation_orig = tx_scale.parameters.parameters.copy()
+    translation_orig = tx_translation.parameters.parameters.copy()
     costvals = []
     for angle in best_angles:
         tx_rotation.parameters.parameters[0] = angle
         tx_rotation.parameters.set_lower_bounds(angle - step / 2)
         tx_rotation.parameters.set_upper_bounds(angle + step / 2)
         tx_rotation.parameters.unlock()
-        tx_scale.parameters.parameters[:] = scale_orig
-        tx_translation.parameters.parameters[:] = translation_orig
+        tx_scale.parameters.parameters[:] = scale_orig.copy()
+        tx_translation.parameters.parameters[:] = translation_orig.copy()
         # Rescale to the same predefined resolution as above
         fixed.resample(float(q.scale), copy=False)
         moving.resample(float(q.scale), copy=False)
         # Set cost function
-        # cost = CostMIND(moving, fixed, maskmode="and", normalise=False)
-        cost = CostMI(moving, fixed, maskmode="and", normalise=False)
+        cost = CostMIND(moving, fixed, maskmode="and", normalise=True)
+        # cost = CostMI(moving, fixed, maskmode="and", normalise=True)
         # Start optimisation
         logger.info(f"Co-optimising scale and translation at "
                     f"{degrees(angle)} deg...")
+
         OptNL(og, cost, method="LN_BOBYQA", visualise=q.visualise,
               xtol_abs=q.xtol_abs, xtol_rel=q.xtol_rel, step=q.opt_step,
-              logger=logger)()
+              logger=logger, normalised=True)()
         costvals.append((cost(), og.parameters[:]))
 
     # Find best initialisation based on cost
@@ -479,13 +488,17 @@ def rigid2d(fixed, moving, cnf):
         tx_scale = fixed_smooth.domain.chain["scale"]
         tx_translation = fixed_smooth.domain.chain["translation"]
         og = OptimisationGroup(tx_rotation, tx_scale, tx_translation)
+        lb, ub = og.get_bounds().T
+        lb = lb - np.finfo(lb.dtype).eps
+        ub = ub + np.finfo(lb.dtype).eps
+        og.set_bounds(lb, ub)
         # Set cost function
-        cost = CostMI(moving_smooth, fixed_smooth, normalise=True)
-        # cost = CostMIND(moving_smooth, fixed_smooth, normalise=False)
+        # cost = CostMI(moving_smooth, fixed_smooth, normalise=True)
+        cost = CostMIND(moving_smooth, fixed_smooth, normalise=True)
         # Start optimisation
         OptNL(og, cost, method="LN_BOBYQA", visualise=q.visualise,
               xtol_abs=q.xtol_abs, xtol_rel=q.xtol_rel, step=q.opt_step,
-              logger=logger)()
+              logger=logger, normalised=True)()
         # Transfer optimised transformations to the non-smoothed images
         fixed.domain = fixed_smooth.domain
         moving.domain = moving_smooth.domain
@@ -516,12 +529,13 @@ def affine2d(fixed, moving, cnf):
         # Prepare transformation to optimise
         tx_affine = fixed_smooth.domain.chain["affine"]
         # Set cost function
-        cost = CostMI(moving_smooth, fixed_smooth, normalise=True)
-        # cost = CostMIND(moving_smooth, fixed_smooth, normalise=False)
+        # cost = CostMI(moving_smooth, fixed_smooth, normalise=True)
+        cost = CostMIND(moving_smooth, fixed_smooth, normalise=True)
         # Start optimisation
         OptNL(tx_affine, cost, method="LN_BOBYQA",
               xtol_rel=q.xtol_rel, xtol_abs=q.xtol_abs,
-              visualise=q.visualise, step=q.opt_step, logger=logger)()
+              visualise=q.visualise, step=q.opt_step, logger=logger,
+              normalised=True)()
         # Transfer optimised transformations to the non-smoothed images
         fixed.domain = fixed_smooth.domain
         moving.domain = moving_smooth.domain
@@ -654,6 +668,13 @@ def block_preprocessing(block, **kwargs):
     return tirl.scripts.mnd.image.rgb2yiq(block.tensors[:3]).tensors[0]
 
 
+def hsv_sat(timg, **kwargs):
+    from skimage.color import rgb2hsv
+    hsv = rgb2hsv(timg.data)
+    return hsv[..., 1]
+    # return TImage.fromarray(hsv[..., 1], tensor_axes=())
+
+
 @beta_function
 def mask_roi_defects(histo, block, p):