From c106bd32bb0468a514b5bb3bd547002622353ea5 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Sun, 7 Mar 2021 23:12:02 +0000 Subject: [PATCH 01/45] Issue #690: add debug prints --- deepreg/loss/deform.py | 1 + deepreg/loss/image.py | 3 +++ 2 files changed, 4 insertions(+) diff --git a/deepreg/loss/deform.py b/deepreg/loss/deform.py index 45d99a285..ea804e590 100644 --- a/deepreg/loss/deform.py +++ b/deepreg/loss/deform.py @@ -85,6 +85,7 @@ def call(self, inputs: tf.Tensor, **kwargs) -> tf.Tensor: :return: shape = () """ assert len(inputs.shape) == 5 + tf.debugging.check_numerics(inputs, "GRAIDENT ddf value NAN/INF", name=None) ddf = inputs # first order gradient # (batch, m_dim1-2, m_dim2-2, m_dim3-2, 3) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 44510bdde..4b8620900 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -231,6 +231,8 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: :return: shape = (batch,) """ # adjust + tf.debugging.check_numerics(y_true, "LNCC y_true value NAN/INF", name=None) + tf.debugging.check_numerics(y_pred, "LNCC y_pred value NAN/INF", name=None) if len(y_true.shape) == 4: y_true = tf.expand_dims(y_true, axis=4) y_pred = tf.expand_dims(y_pred, axis=4) @@ -263,6 +265,7 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] ncc = (cross * cross + EPS) / (t_var * p_var + EPS) + tf.debugging.check_numerics(ncc, "LNCC y_pred value NAN/INF", name=None) return tf.reduce_mean(ncc, axis=[1, 2, 3, 4]) def get_config(self) -> dict: From fac630b8d31fd2a18364766746e52deb68d9df57 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Mon, 22 Mar 2021 00:47:35 +0000 Subject: [PATCH 02/45] Issue #609: hardcode update_freq for tensorboard --- deepreg/train.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/deepreg/train.py b/deepreg/train.py index 5a3180376..6c0c3bb8f 100644 --- a/deepreg/train.py +++ b/deepreg/train.py @@ -142,7 +142,9 @@ def train( # build callbacks tensorboard_callback = tf.keras.callbacks.TensorBoard( - log_dir=log_dir, histogram_freq=config["train"]["save_period"] + log_dir=log_dir, + histogram_freq=config["train"]["save_period"], + update_freq=100, ) ckpt_callback, initial_epoch = build_checkpoint_callback( model=model, From a667dd64e5f458733b131a2c56aebfc26ca10017 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Mon, 22 Mar 2021 20:04:02 +0000 Subject: [PATCH 03/45] do not use zero boundary --- deepreg/model/layer_util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deepreg/model/layer_util.py b/deepreg/model/layer_util.py index 3ceb20024..ae74606de 100644 --- a/deepreg/model/layer_util.py +++ b/deepreg/model/layer_util.py @@ -218,7 +218,7 @@ def resample( vol: tf.Tensor, loc: tf.Tensor, interpolation: str = "linear", - zero_boundary: bool = True, + zero_boundary: bool = False, ) -> tf.Tensor: r""" Sample the volume at given locations. From e0f73c7a4b63b7971db552d806ed842c3ee56c75 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Tue, 23 Mar 2021 22:11:34 +0000 Subject: [PATCH 04/45] Issue #690: modify lncc implementation --- deepreg/loss/image.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 4b8620900..83f608d0b 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -258,15 +258,23 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: p_avg = p_sum / self.kernel_vol # E[p] # shape = (batch, dim1, dim2, dim3, 1) - cross = tp_sum - p_avg * t_sum # E[tp] * E[1] - E[p] * E[t] * E[1] - t_var = t2_sum - t_avg * t_sum # V[t] * E[1] - p_var = p2_sum - p_avg * p_sum # V[p] * E[1] + cross = tp_sum - p_avg * t_sum - t_avg * p_sum + p_avg * t_avg * self.kernel_vol + t_var = t2_sum - 2 * t_avg * t_sum + t_avg * t_avg * self.kernel_vol + p_var = p2_sum - 2 * p_avg * p_sum + p_avg * p_avg * self.kernel_vol + + # cross = tp_sum - p_avg * t_sum # E[tp] * E[1] - E[p] * E[t] * E[1] + # t_var = t2_sum - t_avg * t_sum # V[t] * E[1] + # p_var = p2_sum - p_avg * p_sum # V[p] * E[1] # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] - ncc = (cross * cross + EPS) / (t_var * p_var + EPS) + ncc = (cross * cross) / (t_var * p_var + EPS) + + tf.debugging.check_numerics(ncc, "LNCC ncc value NAN/INF", name=None) + + ncc = tf.reduce_mean(ncc, axis=[1, 2, 3, 4]) + tf.debugging.check_numerics(ncc, "LNCC ncc_mean value NAN/INF", name=None) - tf.debugging.check_numerics(ncc, "LNCC y_pred value NAN/INF", name=None) - return tf.reduce_mean(ncc, axis=[1, 2, 3, 4]) + return ncc def get_config(self) -> dict: """Return the config dictionary for recreating this class.""" From b141b3cc4d5123e07d827e5f0cdbd3df4eb4daa4 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Tue, 23 Mar 2021 23:01:03 +0000 Subject: [PATCH 05/45] Issue #690: add EPS to numerator for LNCC --- deepreg/loss/image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 83f608d0b..08bc397d1 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -267,7 +267,7 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # p_var = p2_sum - p_avg * p_sum # V[p] * E[1] # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] - ncc = (cross * cross) / (t_var * p_var + EPS) + ncc = (cross * cross + EPS) / (t_var * p_var + EPS) tf.debugging.check_numerics(ncc, "LNCC ncc value NAN/INF", name=None) From 952ae52db08ae4aa0f1da7832dc1dd8c57486190 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Tue, 23 Mar 2021 23:18:54 +0000 Subject: [PATCH 06/45] Issue #690: add debug print --- deepreg/loss/image.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 08bc397d1..08f687e31 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -274,6 +274,16 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: ncc = tf.reduce_mean(ncc, axis=[1, 2, 3, 4]) tf.debugging.check_numerics(ncc, "LNCC ncc_mean value NAN/INF", name=None) + tf.print( + "DEBUG", + tf.reduce_min(ncc), + tf.reduce_max(ncc), + tf.reduce_min(y_pred), + tf.reduce_max(y_pred), + tf.reduce_min(y_true), + tf.reduce_max(y_true), + ncc, + ) return ncc def get_config(self) -> dict: From 9f997d1041dd01892329bce129dff3127e9a0b87 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Tue, 23 Mar 2021 23:32:32 +0000 Subject: [PATCH 07/45] Issue #690: use tf.debugging.enable_check_numerics --- deepreg/loss/image.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 08f687e31..117ccc2e3 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -231,6 +231,7 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: :return: shape = (batch,) """ # adjust + tf.debugging.enable_check_numerics() tf.debugging.check_numerics(y_true, "LNCC y_true value NAN/INF", name=None) tf.debugging.check_numerics(y_pred, "LNCC y_pred value NAN/INF", name=None) if len(y_true.shape) == 4: @@ -274,16 +275,16 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: ncc = tf.reduce_mean(ncc, axis=[1, 2, 3, 4]) tf.debugging.check_numerics(ncc, "LNCC ncc_mean value NAN/INF", name=None) - tf.print( - "DEBUG", - tf.reduce_min(ncc), - tf.reduce_max(ncc), - tf.reduce_min(y_pred), - tf.reduce_max(y_pred), - tf.reduce_min(y_true), - tf.reduce_max(y_true), - ncc, - ) + # tf.print( + # "DEBUG", + # tf.reduce_min(ncc), + # tf.reduce_max(ncc), + # tf.reduce_min(y_pred), + # tf.reduce_max(y_pred), + # tf.reduce_min(y_true), + # tf.reduce_max(y_true), + # ncc, + # ) return ncc def get_config(self) -> dict: From 8ff1898ccbd4af0f36c37720e3594d5bb3394f77 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Wed, 24 Mar 2021 01:27:28 +0000 Subject: [PATCH 08/45] Issue #690: attempt to skip some ops --- deepreg/model/network.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/deepreg/model/network.py b/deepreg/model/network.py index 3de4e10e3..3de45eed8 100644 --- a/deepreg/model/network.py +++ b/deepreg/model/network.py @@ -60,7 +60,8 @@ def __init__( self.batch_size = batch_size self.config = config self.num_devices = num_devices - self.global_batch_size = num_devices * batch_size + self.global_batch_size = num_devices * batch_size * 1.0 + assert self.global_batch_size > 0 self._inputs = None # save inputs of self._model as dict self._outputs = None # save outputs of self._model as dict @@ -215,7 +216,7 @@ def _build_loss(self, name: str, inputs_dict: dict): config=dict_without(d=loss_config, key="weight") ) loss_value = loss_layer(**inputs_dict) / self.global_batch_size - weighted_loss = loss_value * weight + weighted_loss = loss_value * weight if weight != 1 else loss_value # add loss self._model.add_loss(weighted_loss) From 1330286cf7c3eb77ac530ca105ce434bbf17e6e9 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Wed, 24 Mar 2021 22:43:05 +0000 Subject: [PATCH 09/45] Issue #690: check input image shapes --- deepreg/dataset/loader/interface.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/deepreg/dataset/loader/interface.py b/deepreg/dataset/loader/interface.py index 1c89d5e53..bb169d047 100644 --- a/deepreg/dataset/loader/interface.py +++ b/deepreg/dataset/loader/interface.py @@ -384,6 +384,11 @@ def validate_images_and_labels( f"Sample {image_indices}'s {name}' shape should be 3D. " f"Got {arr.shape}." ) + if min(arr.shape) <= 0: + raise ValueError( + f"Sample {image_indices}'s {name}' shape should be non-empty. " + f"Got {arr.shape}." + ) # when data are labeled if moving_label is not None and fixed_label is not None: # labels should be 3D or 4D arrays From e2354ff79bbbdbe50de139db82d189d6d4aee82a Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Thu, 25 Mar 2021 01:34:50 +0000 Subject: [PATCH 10/45] Issue #690: add debug check --- deepreg/model/network.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/deepreg/model/network.py b/deepreg/model/network.py index 7800f9794..dc766b97f 100644 --- a/deepreg/model/network.py +++ b/deepreg/model/network.py @@ -224,6 +224,12 @@ def _build_loss(self, name: str, inputs_dict: dict): # add loss self._model.add_loss(weighted_loss) + tf.debugging.check_numerics( + loss_value, + f"loss {name}_{loss_layer.name} inf/nan", + name=f"op loss/{name}_{loss_layer.name}", + ) + # add metric self._model.add_metric( loss_value, name=f"loss/{name}_{loss_layer.name}", aggregation="mean" From f31f532d23228e9478a65ad83fe7d5dcecb08206 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Thu, 25 Mar 2021 01:35:54 +0000 Subject: [PATCH 11/45] Issue #690: fix op name --- deepreg/model/network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deepreg/model/network.py b/deepreg/model/network.py index dc766b97f..0f63cd958 100644 --- a/deepreg/model/network.py +++ b/deepreg/model/network.py @@ -227,7 +227,7 @@ def _build_loss(self, name: str, inputs_dict: dict): tf.debugging.check_numerics( loss_value, f"loss {name}_{loss_layer.name} inf/nan", - name=f"op loss/{name}_{loss_layer.name}", + name=f"op_loss_{name}_{loss_layer.name}", ) # add metric From 0ad5f56f55c931494c252c63f3f4b92bc6986abd Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Thu, 25 Mar 2021 01:37:38 +0000 Subject: [PATCH 12/45] Issue #690: add check into graph --- deepreg/model/network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deepreg/model/network.py b/deepreg/model/network.py index 0f63cd958..f2a577853 100644 --- a/deepreg/model/network.py +++ b/deepreg/model/network.py @@ -224,7 +224,7 @@ def _build_loss(self, name: str, inputs_dict: dict): # add loss self._model.add_loss(weighted_loss) - tf.debugging.check_numerics( + loss_value = tf.debugging.check_numerics( loss_value, f"loss {name}_{loss_layer.name} inf/nan", name=f"op_loss_{name}_{loss_layer.name}", From 40be7c3ba04af7182eec15fb4d3f36f93f8ac8f2 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Thu, 25 Mar 2021 01:41:24 +0000 Subject: [PATCH 13/45] Issue #690: add check into lncc graph --- deepreg/loss/image.py | 29 +++++++++++++---------------- 1 file changed, 13 insertions(+), 16 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 117ccc2e3..f4cb0b472 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -231,9 +231,12 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: :return: shape = (batch,) """ # adjust - tf.debugging.enable_check_numerics() - tf.debugging.check_numerics(y_true, "LNCC y_true value NAN/INF", name=None) - tf.debugging.check_numerics(y_pred, "LNCC y_pred value NAN/INF", name=None) + y_true = tf.debugging.check_numerics( + y_true, "LNCC y_true value NAN/INF", name="LNCC_y_true" + ) + y_pred = tf.debugging.check_numerics( + y_pred, "LNCC y_pred value NAN/INF", name="LNCCC_y_pred" + ) if len(y_true.shape) == 4: y_true = tf.expand_dims(y_true, axis=4) y_pred = tf.expand_dims(y_pred, axis=4) @@ -270,21 +273,15 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] ncc = (cross * cross + EPS) / (t_var * p_var + EPS) - tf.debugging.check_numerics(ncc, "LNCC ncc value NAN/INF", name=None) + ncc = tf.debugging.check_numerics( + ncc, "LNCC ncc value NAN/INF", name="LNCC_before_mean" + ) ncc = tf.reduce_mean(ncc, axis=[1, 2, 3, 4]) - tf.debugging.check_numerics(ncc, "LNCC ncc_mean value NAN/INF", name=None) - - # tf.print( - # "DEBUG", - # tf.reduce_min(ncc), - # tf.reduce_max(ncc), - # tf.reduce_min(y_pred), - # tf.reduce_max(y_pred), - # tf.reduce_min(y_true), - # tf.reduce_max(y_true), - # ncc, - # ) + ncc = tf.debugging.check_numerics( + ncc, "LNCC ncc_mean value NAN/INF", name="LNCC_after_mean" + ) + return ncc def get_config(self) -> dict: From e0345710c627a9a68f17da3edac2b59a8be8f144 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Thu, 25 Mar 2021 01:59:44 +0000 Subject: [PATCH 14/45] Issue #690: add more checks into lncc --- deepreg/loss/image.py | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index f4cb0b472..86fb706f0 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -256,11 +256,34 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: p2_sum = separable_filter(p2, kernel=self.kernel) # E[pp] * E[1] tp_sum = separable_filter(tp, kernel=self.kernel) # E[tp] * E[1] + t_sum = tf.debugging.check_numerics( + t_sum, "LNCC t_sum value NAN/INF", name="LNCC_t_sum" + ) + p_sum = tf.debugging.check_numerics( + p_sum, "LNCC p_avg value NAN/INF", name="LNCC_p_sum" + ) + t2_sum = tf.debugging.check_numerics( + t2_sum, "LNCC t2_sum value NAN/INF", name="LNCC_t2_sum" + ) + p2_sum = tf.debugging.check_numerics( + p2_sum, "LNCC p2_sum value NAN/INF", name="LNCC_p2_sum" + ) + tp_sum = tf.debugging.check_numerics( + tp_sum, "LNCC tp_sum value NAN/INF", name="LNCC_tp_sum" + ) + # average over kernel # (batch, dim1, dim2, dim3, 1) t_avg = t_sum / self.kernel_vol # E[t] p_avg = p_sum / self.kernel_vol # E[p] + t_avg = tf.debugging.check_numerics( + t_avg, "LNCC t_avg value NAN/INF", name="LNCC_t_avg" + ) + p_avg = tf.debugging.check_numerics( + p_avg, "LNCC p_avg value NAN/INF", name="LNCC_p_avg" + ) + # shape = (batch, dim1, dim2, dim3, 1) cross = tp_sum - p_avg * t_sum - t_avg * p_sum + p_avg * t_avg * self.kernel_vol t_var = t2_sum - 2 * t_avg * t_sum + t_avg * t_avg * self.kernel_vol @@ -270,6 +293,16 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # t_var = t2_sum - t_avg * t_sum # V[t] * E[1] # p_var = p2_sum - p_avg * p_sum # V[p] * E[1] + cross = tf.debugging.check_numerics( + cross, "LNCC cross value NAN/INF", name="LNCC_cross" + ) + t_var = tf.debugging.check_numerics( + t_var, "LNCC t_var value NAN/INF", name="LNCC_t_var" + ) + p_var = tf.debugging.check_numerics( + p_var, "LNCC t_var value NAN/INF", name="LNCC_t_var" + ) + # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] ncc = (cross * cross + EPS) / (t_var * p_var + EPS) From eebfbddbd47ed4690d2055e4fa609dc2415c9ca9 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Thu, 25 Mar 2021 02:48:35 +0000 Subject: [PATCH 15/45] Issue #690: assert denom >= 0 --- deepreg/loss/image.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 86fb706f0..422f7fe8b 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -306,6 +306,23 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] ncc = (cross * cross + EPS) / (t_var * p_var + EPS) + tf.debugging.Assert( + tf.greater_equal(t_var * p_var, 0), + [ + tf.reduce_min(cross), + tf.reduce_max(cross), + tf.reduce_min(t_var), + tf.reduce_max(t_var), + tf.reduce_min(p_var), + tf.reduce_max(p_var), + cross, + t_var, + p_var, + ], + summarize=None, + name="debug_div", + ) + ncc = tf.debugging.check_numerics( ncc, "LNCC ncc value NAN/INF", name="LNCC_before_mean" ) From 4dc9a08dfab097792b2e839bca2e2c76983911d1 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Thu, 25 Mar 2021 20:32:51 +0000 Subject: [PATCH 16/45] Issue #690: assert based on nan --- deepreg/loss/image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 422f7fe8b..8aa316baf 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -307,7 +307,7 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: ncc = (cross * cross + EPS) / (t_var * p_var + EPS) tf.debugging.Assert( - tf.greater_equal(t_var * p_var, 0), + tf.math.reduce_any(tf.math.is_inf(ncc)), [ tf.reduce_min(cross), tf.reduce_max(cross), From b2f63382c6f3af2bb290aecda03d9dddc3559260 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Thu, 25 Mar 2021 22:34:26 +0000 Subject: [PATCH 17/45] Issue #690: do proper debugging --- deepreg/loss/image.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 8aa316baf..9eb5ab6a8 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -306,9 +306,8 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] ncc = (cross * cross + EPS) / (t_var * p_var + EPS) - tf.debugging.Assert( - tf.math.reduce_any(tf.math.is_inf(ncc)), - [ + debug_assert = tf.debugging.Assert( + tf.math.is_finite(tf.reduce_mean(ncc))[ tf.reduce_min(cross), tf.reduce_max(cross), tf.reduce_min(t_var), @@ -322,6 +321,7 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: summarize=None, name="debug_div", ) + debug_assert.mark_used() ncc = tf.debugging.check_numerics( ncc, "LNCC ncc value NAN/INF", name="LNCC_before_mean" From e0219d85d02603ae19f8a8faf9718434fa938787 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Thu, 25 Mar 2021 22:36:27 +0000 Subject: [PATCH 18/45] Issue #690: fix typo --- deepreg/loss/image.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 9eb5ab6a8..1b554ebd9 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -307,7 +307,8 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: ncc = (cross * cross + EPS) / (t_var * p_var + EPS) debug_assert = tf.debugging.Assert( - tf.math.is_finite(tf.reduce_mean(ncc))[ + tf.math.is_finite(tf.reduce_mean(ncc)), + [ tf.reduce_min(cross), tf.reduce_max(cross), tf.reduce_min(t_var), From 90aa4344322024bd892f784421a67a1a17014f64 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Thu, 25 Mar 2021 22:39:11 +0000 Subject: [PATCH 19/45] Issue #690: add assert into graph --- deepreg/loss/image.py | 36 ++++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 1b554ebd9..e53b3250d 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -306,23 +306,27 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] ncc = (cross * cross + EPS) / (t_var * p_var + EPS) - debug_assert = tf.debugging.Assert( - tf.math.is_finite(tf.reduce_mean(ncc)), + with tf.control_dependencies( [ - tf.reduce_min(cross), - tf.reduce_max(cross), - tf.reduce_min(t_var), - tf.reduce_max(t_var), - tf.reduce_min(p_var), - tf.reduce_max(p_var), - cross, - t_var, - p_var, - ], - summarize=None, - name="debug_div", - ) - debug_assert.mark_used() + tf.debugging.Assert( + tf.math.is_finite(tf.reduce_mean(ncc)), + [ + tf.reduce_min(cross), + tf.reduce_max(cross), + tf.reduce_min(t_var), + tf.reduce_max(t_var), + tf.reduce_min(p_var), + tf.reduce_max(p_var), + cross, + t_var, + p_var, + ], + summarize=None, + name="debug_div", + ) + ] + ): + ncc = tf.identity(ncc) ncc = tf.debugging.check_numerics( ncc, "LNCC ncc value NAN/INF", name="LNCC_before_mean" From 6c27061af1ea056917e9d4b8f8fe85d2ea6ba2ce Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 01:32:33 +0000 Subject: [PATCH 20/45] Issue #690: revert change to LNCC and increase EPS --- deepreg/loss/image.py | 80 +++---------------------------------------- deepreg/loss/label.py | 2 +- deepreg/loss/util.py | 2 +- 3 files changed, 7 insertions(+), 77 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index e53b3250d..975e7eb5d 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -10,7 +10,7 @@ ) from deepreg.registry import REGISTRY -EPS = tf.keras.backend.epsilon() +EPS = 1.0e-5 @REGISTRY.register_loss(name="ssd") @@ -231,12 +231,6 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: :return: shape = (batch,) """ # adjust - y_true = tf.debugging.check_numerics( - y_true, "LNCC y_true value NAN/INF", name="LNCC_y_true" - ) - y_pred = tf.debugging.check_numerics( - y_pred, "LNCC y_pred value NAN/INF", name="LNCCC_y_pred" - ) if len(y_true.shape) == 4: y_true = tf.expand_dims(y_true, axis=4) y_pred = tf.expand_dims(y_pred, axis=4) @@ -256,88 +250,24 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: p2_sum = separable_filter(p2, kernel=self.kernel) # E[pp] * E[1] tp_sum = separable_filter(tp, kernel=self.kernel) # E[tp] * E[1] - t_sum = tf.debugging.check_numerics( - t_sum, "LNCC t_sum value NAN/INF", name="LNCC_t_sum" - ) - p_sum = tf.debugging.check_numerics( - p_sum, "LNCC p_avg value NAN/INF", name="LNCC_p_sum" - ) - t2_sum = tf.debugging.check_numerics( - t2_sum, "LNCC t2_sum value NAN/INF", name="LNCC_t2_sum" - ) - p2_sum = tf.debugging.check_numerics( - p2_sum, "LNCC p2_sum value NAN/INF", name="LNCC_p2_sum" - ) - tp_sum = tf.debugging.check_numerics( - tp_sum, "LNCC tp_sum value NAN/INF", name="LNCC_tp_sum" - ) - # average over kernel # (batch, dim1, dim2, dim3, 1) t_avg = t_sum / self.kernel_vol # E[t] p_avg = p_sum / self.kernel_vol # E[p] - t_avg = tf.debugging.check_numerics( - t_avg, "LNCC t_avg value NAN/INF", name="LNCC_t_avg" - ) - p_avg = tf.debugging.check_numerics( - p_avg, "LNCC p_avg value NAN/INF", name="LNCC_p_avg" - ) - # shape = (batch, dim1, dim2, dim3, 1) - cross = tp_sum - p_avg * t_sum - t_avg * p_sum + p_avg * t_avg * self.kernel_vol - t_var = t2_sum - 2 * t_avg * t_sum + t_avg * t_avg * self.kernel_vol - p_var = p2_sum - 2 * p_avg * p_sum + p_avg * p_avg * self.kernel_vol - - # cross = tp_sum - p_avg * t_sum # E[tp] * E[1] - E[p] * E[t] * E[1] - # t_var = t2_sum - t_avg * t_sum # V[t] * E[1] - # p_var = p2_sum - p_avg * p_sum # V[p] * E[1] - - cross = tf.debugging.check_numerics( - cross, "LNCC cross value NAN/INF", name="LNCC_cross" - ) - t_var = tf.debugging.check_numerics( - t_var, "LNCC t_var value NAN/INF", name="LNCC_t_var" - ) - p_var = tf.debugging.check_numerics( - p_var, "LNCC t_var value NAN/INF", name="LNCC_t_var" - ) + cross = tp_sum - p_avg * t_sum # E[tp] * E[1] - E[p] * E[t] * E[1] + t_var = t2_sum - t_avg * t_sum # V[t] * E[1] + p_var = p2_sum - p_avg * p_sum # V[p] * E[1] # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] ncc = (cross * cross + EPS) / (t_var * p_var + EPS) - with tf.control_dependencies( - [ - tf.debugging.Assert( - tf.math.is_finite(tf.reduce_mean(ncc)), - [ - tf.reduce_min(cross), - tf.reduce_max(cross), - tf.reduce_min(t_var), - tf.reduce_max(t_var), - tf.reduce_min(p_var), - tf.reduce_max(p_var), - cross, - t_var, - p_var, - ], - summarize=None, - name="debug_div", - ) - ] - ): - ncc = tf.identity(ncc) - ncc = tf.debugging.check_numerics( ncc, "LNCC ncc value NAN/INF", name="LNCC_before_mean" ) - ncc = tf.reduce_mean(ncc, axis=[1, 2, 3, 4]) - ncc = tf.debugging.check_numerics( - ncc, "LNCC ncc_mean value NAN/INF", name="LNCC_after_mean" - ) - - return ncc + return tf.reduce_mean(ncc, axis=[1, 2, 3, 4]) def get_config(self) -> dict: """Return the config dictionary for recreating this class.""" diff --git a/deepreg/loss/label.py b/deepreg/loss/label.py index 0e0d94bad..6b7693e30 100644 --- a/deepreg/loss/label.py +++ b/deepreg/loss/label.py @@ -9,7 +9,7 @@ from deepreg.loss.util import separable_filter from deepreg.registry import REGISTRY -EPS = tf.keras.backend.epsilon() +EPS = 1.0e-5 class MultiScaleLoss(tf.keras.losses.Loss): diff --git a/deepreg/loss/util.py b/deepreg/loss/util.py index 6f453d73c..909d9d49d 100644 --- a/deepreg/loss/util.py +++ b/deepreg/loss/util.py @@ -27,7 +27,7 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: return -super().call(y_true=y_true, y_pred=y_pred) -EPS = tf.keras.backend.epsilon() +EPS = 1.0e-5 def rectangular_kernel1d(kernel_size: int) -> tf.Tensor: From de96cad667ff14077a3539aa993d8f05f7ebaef8 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 01:47:30 +0000 Subject: [PATCH 21/45] Issue #690: add additional metrics --- deepreg/loss/image.py | 135 +++++++++++++++++++++++++++++++++++++++ deepreg/model/network.py | 5 ++ 2 files changed, 140 insertions(+) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 975e7eb5d..bbccb50e8 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -277,6 +277,141 @@ def get_config(self) -> dict: return config +class LocalNormalizedCrossCorrelationDEBUG(tf.keras.layers.Layer): + """ + Local squared zero-normalized cross-correlation. + + Denote y_true as t and y_pred as p. Consider a window having n elements. + Each position in the window corresponds a weight w_i for i=1:n. + + Define the discrete expectation in the window E[t] as + + E[t] = sum_i(w_i * t_i) / sum_i(w_i) + + Similarly, the discrete variance in the window V[t] is + + V[t] = E[t**2] - E[t] ** 2 + + The local squared zero-normalized cross-correlation is therefore + + E[ (t-E[t]) * (p-E[p]) ] ** 2 / V[t] / V[p] + + where the expectation in numerator is + + E[ (t-E[t]) * (p-E[p]) ] = E[t * p] - E[t] * E[p] + + Different kernel corresponds to different weights. + + For now, y_true and y_pred have to be at least 4d tensor, including batch axis. + + Reference: + + - Zero-normalized cross-correlation (ZNCC): + https://en.wikipedia.org/wiki/Cross-correlation + - Code: https://github.com/voxelmorph/voxelmorph/blob/legacy/src/losses.py + """ + + kernel_fn_dict = dict( + gaussian=gaussian_kernel1d, + rectangular=rectangular_kernel1d, + triangular=triangular_kernel1d, + ) + + def __init__( + self, + kernel_size: int = 9, + kernel_type: str = "rectangular", + reduction: str = tf.keras.losses.Reduction.SUM, + name: str = "LocalNormalizedCrossCorrelationDEBUG", + ): + """ + Init. + + :param kernel_size: int. Kernel size or kernel sigma for kernel_type='gauss'. + :param kernel_type: str, rectangular, triangular or gaussian + :param reduction: using SUM reduction over batch axis, + calling the loss like `loss(y_true, y_pred)` will return a scalar tensor. + :param name: name of the loss + """ + super().__init__(reduction=reduction, name=name) + if kernel_type not in self.kernel_fn_dict.keys(): + raise ValueError( + f"Wrong kernel_type {kernel_type} for LNCC loss type. " + f"Feasible values are {self.kernel_fn_dict.keys()}" + ) + self.kernel_fn = self.kernel_fn_dict[kernel_type] + self.kernel_type = kernel_type + self.kernel_size = kernel_size + + # (kernel_size, ) + self.kernel = self.kernel_fn(kernel_size=self.kernel_size) + # E[1] = sum_i(w_i), () + self.kernel_vol = tf.reduce_sum( + self.kernel[:, None, None] + * self.kernel[None, :, None] + * self.kernel[None, None, :] + ) + + def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: + """ + Return loss for a batch. + + :param y_true: shape = (batch, dim1, dim2, dim3) + or (batch, dim1, dim2, dim3, ch) + :param y_pred: shape = (batch, dim1, dim2, dim3) + or (batch, dim1, dim2, dim3, ch) + :return: shape = (batch,) + """ + # adjust + if len(y_true.shape) == 4: + y_true = tf.expand_dims(y_true, axis=4) + y_pred = tf.expand_dims(y_pred, axis=4) + assert len(y_true.shape) == len(y_pred.shape) == 5 + + # t = y_true, p = y_pred + # (batch, dim1, dim2, dim3, ch) + t2 = y_true * y_true + p2 = y_pred * y_pred + tp = y_true * y_pred + + # sum over kernel + # (batch, dim1, dim2, dim3, 1) + t_sum = separable_filter(y_true, kernel=self.kernel) # E[t] * E[1] + p_sum = separable_filter(y_pred, kernel=self.kernel) # E[p] * E[1] + t2_sum = separable_filter(t2, kernel=self.kernel) # E[tt] * E[1] + p2_sum = separable_filter(p2, kernel=self.kernel) # E[pp] * E[1] + tp_sum = separable_filter(tp, kernel=self.kernel) # E[tp] * E[1] + + # average over kernel + # (batch, dim1, dim2, dim3, 1) + t_avg = t_sum / self.kernel_vol # E[t] + p_avg = p_sum / self.kernel_vol # E[p] + + # shape = (batch, dim1, dim2, dim3, 1) + cross = tp_sum - p_avg * t_sum # E[tp] * E[1] - E[p] * E[t] * E[1] + t_var = t2_sum - t_avg * t_sum # V[t] * E[1] + p_var = p2_sum - p_avg * p_sum # V[p] * E[1] + + # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] + ncc = t_var * p_var + EPS + + ncc = tf.debugging.check_numerics( + ncc, "DEBUG LNCC ncc value NAN/INF", name="DEBUGLNCC_before_mean" + ) + cross = tf.debugging.check_numerics( + cross, "DEBUG LNCC ncc value NAN/INF", name="DEBUGLNCC_before_mean" + ) + + return tf.reduce_mean(ncc, axis=[1, 2, 3, 4]) + + def get_config(self) -> dict: + """Return the config dictionary for recreating this class.""" + config = super().get_config() + config["kernel_size"] = self.kernel_size + config["kernel_type"] = self.kernel_type + return config + + @REGISTRY.register_loss(name="lncc") class LocalNormalizedCrossCorrelationLoss( NegativeLossMixin, LocalNormalizedCrossCorrelation diff --git a/deepreg/model/network.py b/deepreg/model/network.py index f2a577853..a06a30c5e 100644 --- a/deepreg/model/network.py +++ b/deepreg/model/network.py @@ -6,6 +6,7 @@ import tensorflow as tf +from deepreg.loss.image import LocalNormalizedCrossCorrelationDEBUG from deepreg.loss.label import DiceScore, compute_centroid_distance from deepreg.model import layer, layer_util from deepreg.model.backbone import GlobalNet @@ -253,6 +254,10 @@ def build_loss(self): # image loss, conditional model does not have this if "pred_fixed_image" in self._outputs: pred_fixed_image = self._outputs["pred_fixed_image"] + debug_value = LocalNormalizedCrossCorrelationDEBUG()( + y_true=fixed_image, y_pred=pred_fixed_image + ) + self.log_tensor_stats(debug_value, name="debug-lncc") self._build_loss( name="image", inputs_dict=dict(y_true=fixed_image, y_pred=pred_fixed_image), From 64e6561bf9a60c35ed43ab63df50cbcf82d778d5 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 01:50:54 +0000 Subject: [PATCH 22/45] Issue #690: remove arg --- deepreg/loss/image.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index bbccb50e8..3467124c5 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -321,7 +321,7 @@ def __init__( self, kernel_size: int = 9, kernel_type: str = "rectangular", - reduction: str = tf.keras.losses.Reduction.SUM, + # reduction: str = tf.keras.losses.Reduction.SUM, name: str = "LocalNormalizedCrossCorrelationDEBUG", ): """ @@ -329,11 +329,9 @@ def __init__( :param kernel_size: int. Kernel size or kernel sigma for kernel_type='gauss'. :param kernel_type: str, rectangular, triangular or gaussian - :param reduction: using SUM reduction over batch axis, - calling the loss like `loss(y_true, y_pred)` will return a scalar tensor. :param name: name of the loss """ - super().__init__(reduction=reduction, name=name) + super().__init__(name=name) if kernel_type not in self.kernel_fn_dict.keys(): raise ValueError( f"Wrong kernel_type {kernel_type} for LNCC loss type. " From 15631dab898effb0c7a38f8185a6ec9740617fa5 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 01:55:40 +0000 Subject: [PATCH 23/45] Issue #690: fit eagerly --- deepreg/loss/image.py | 5 +++++ deepreg/train.py | 2 +- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 3467124c5..1ac05145f 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -263,6 +263,11 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] ncc = (cross * cross + EPS) / (t_var * p_var + EPS) + print("cross", cross) + print("t_var", t_var) + print("p_var", p_var) + print("ncc", ncc) + ncc = tf.debugging.check_numerics( ncc, "LNCC ncc value NAN/INF", name="LNCC_before_mean" ) diff --git a/deepreg/train.py b/deepreg/train.py index 60aeabac2..6de154ac0 100644 --- a/deepreg/train.py +++ b/deepreg/train.py @@ -137,7 +137,7 @@ def train( optimizer = opt.build_optimizer(optimizer_config=config["train"]["optimizer"]) # compile - model.compile(optimizer=optimizer) + model.compile(optimizer=optimizer, run_eagerly=True) model.plot_model(output_dir=log_dir) # build callbacks From 1cbf23c9242616672ccbf8fc1bc625d6fec12a14 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 02:02:31 +0000 Subject: [PATCH 24/45] Issue #690: print min max and save array --- deepreg/loss/image.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 1ac05145f..246319bef 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -1,4 +1,5 @@ """Provide different loss or metrics classes for images.""" +import numpy as np import tensorflow as tf from deepreg.loss.util import NegativeLossMixin @@ -263,10 +264,14 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] ncc = (cross * cross + EPS) / (t_var * p_var + EPS) - print("cross", cross) - print("t_var", t_var) - print("p_var", p_var) - print("ncc", ncc) + for x_name, x in [ + ("cross", cross), + ("t_var", t_var), + ("p_var", p_var), + ("ncc", ncc), + ]: + print(x_name + "_min_max", tf.reduce_min(x), tf.reduce_max(x)) + np.save(x_name + ".npy", x.numpy()) ncc = tf.debugging.check_numerics( ncc, "LNCC ncc value NAN/INF", name="LNCC_before_mean" From 5ac34c1041bbc067d97423a83f442f9f54e0da78 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 02:04:47 +0000 Subject: [PATCH 25/45] Issue #690: catch err --- deepreg/loss/image.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 246319bef..c17953a49 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -271,7 +271,10 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: ("ncc", ncc), ]: print(x_name + "_min_max", tf.reduce_min(x), tf.reduce_max(x)) - np.save(x_name + ".npy", x.numpy()) + try: + np.save(x_name + ".npy", x.numpy()) + except AttributeError: + print("save failed") ncc = tf.debugging.check_numerics( ncc, "LNCC ncc value NAN/INF", name="LNCC_before_mean" From c40cb34b5c4ea4ed07cdca7f95661c630faf7e94 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 02:36:08 +0000 Subject: [PATCH 26/45] add print into grapj --- deepreg/loss/image.py | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index c17953a49..361841324 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -264,17 +264,22 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] ncc = (cross * cross + EPS) / (t_var * p_var + EPS) - for x_name, x in [ - ("cross", cross), - ("t_var", t_var), - ("p_var", p_var), - ("ncc", ncc), - ]: - print(x_name + "_min_max", tf.reduce_min(x), tf.reduce_max(x)) - try: - np.save(x_name + ".npy", x.numpy()) - except AttributeError: - print("save failed") + @tf.function + def verify(): + for x_name, x in [ + ("cross", cross), + ("t_var", t_var), + ("p_var", p_var), + ("ncc", ncc), + ]: + tf.print(x_name + "_min_max", tf.reduce_min(x), tf.reduce_max(x)) + try: + np.save(x_name + ".npy", x.numpy()) + except AttributeError: + print("save failed") + return ncc + + ncc = verify() ncc = tf.debugging.check_numerics( ncc, "LNCC ncc value NAN/INF", name="LNCC_before_mean" From 06b76b01ad2f738d59e2cef1671fdb4a81054fad Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 02:49:13 +0000 Subject: [PATCH 27/45] Issue #690: hack ncc --- deepreg/loss/image.py | 22 +++------------------- 1 file changed, 3 insertions(+), 19 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 361841324..57e68acf2 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -1,5 +1,4 @@ """Provide different loss or metrics classes for images.""" -import numpy as np import tensorflow as tf from deepreg.loss.util import NegativeLossMixin @@ -262,24 +261,9 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: p_var = p2_sum - p_avg * p_sum # V[p] * E[1] # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] - ncc = (cross * cross + EPS) / (t_var * p_var + EPS) - - @tf.function - def verify(): - for x_name, x in [ - ("cross", cross), - ("t_var", t_var), - ("p_var", p_var), - ("ncc", ncc), - ]: - tf.print(x_name + "_min_max", tf.reduce_min(x), tf.reduce_max(x)) - try: - np.save(x_name + ".npy", x.numpy()) - except AttributeError: - print("save failed") - return ncc - - ncc = verify() + num = cross * cross + EPS + denom = t_var * p_var + EPS + ncc = tf.math.minimum(num, denom) / denom ncc = tf.debugging.check_numerics( ncc, "LNCC ncc value NAN/INF", name="LNCC_before_mean" From 259bbd9c34a4d8d562433c83a02d0c57b79fa67b Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 02:50:50 +0000 Subject: [PATCH 28/45] Issue #690: do not compile eagerly --- deepreg/train.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deepreg/train.py b/deepreg/train.py index 6de154ac0..a2fba6cdb 100644 --- a/deepreg/train.py +++ b/deepreg/train.py @@ -137,7 +137,7 @@ def train( optimizer = opt.build_optimizer(optimizer_config=config["train"]["optimizer"]) # compile - model.compile(optimizer=optimizer, run_eagerly=True) + model.compile(optimizer=optimizer) # , run_eagerly=True) model.plot_model(output_dir=log_dir) # build callbacks From 9e7203b1ffcfbe134f62a5e7b1af222b10390f56 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 03:09:00 +0000 Subject: [PATCH 29/45] Issue #690: add metrics --- deepreg/loss/image.py | 147 ++++++++++++++++++++++++++++++++++++--- deepreg/model/network.py | 13 +++- 2 files changed, 148 insertions(+), 12 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 57e68acf2..6e63d731d 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -263,7 +263,7 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] num = cross * cross + EPS denom = t_var * p_var + EPS - ncc = tf.math.minimum(num, denom) / denom + ncc = num / denom ncc = tf.debugging.check_numerics( ncc, "LNCC ncc value NAN/INF", name="LNCC_before_mean" @@ -279,7 +279,7 @@ def get_config(self) -> dict: return config -class LocalNormalizedCrossCorrelationDEBUG(tf.keras.layers.Layer): +class LocalNormalizedCrossCorrelationDEBUG1(tf.keras.layers.Layer): """ Local squared zero-normalized cross-correlation. @@ -324,7 +324,7 @@ def __init__( kernel_size: int = 9, kernel_type: str = "rectangular", # reduction: str = tf.keras.losses.Reduction.SUM, - name: str = "LocalNormalizedCrossCorrelationDEBUG", + name: str = "LocalNormalizedCrossCorrelationDEBUG1", ): """ Init. @@ -393,16 +393,145 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: p_var = p2_sum - p_avg * p_sum # V[p] * E[1] # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] - ncc = t_var * p_var + EPS + num = cross * cross + EPS + denom = t_var * p_var + EPS - ncc = tf.debugging.check_numerics( - ncc, "DEBUG LNCC ncc value NAN/INF", name="DEBUGLNCC_before_mean" + denom = tf.debugging.check_numerics( + denom, "DEBUG LNCC denom value NAN/INF", name="DEBUGLNCC_before_mean_denom" ) - cross = tf.debugging.check_numerics( - cross, "DEBUG LNCC ncc value NAN/INF", name="DEBUGLNCC_before_mean" + + return num + + def get_config(self) -> dict: + """Return the config dictionary for recreating this class.""" + config = super().get_config() + config["kernel_size"] = self.kernel_size + config["kernel_type"] = self.kernel_type + return config + + +class LocalNormalizedCrossCorrelationDEBUG2(tf.keras.layers.Layer): + """ + Local squared zero-normalized cross-correlation. + + Denote y_true as t and y_pred as p. Consider a window having n elements. + Each position in the window corresponds a weight w_i for i=1:n. + + Define the discrete expectation in the window E[t] as + + E[t] = sum_i(w_i * t_i) / sum_i(w_i) + + Similarly, the discrete variance in the window V[t] is + + V[t] = E[t**2] - E[t] ** 2 + + The local squared zero-normalized cross-correlation is therefore + + E[ (t-E[t]) * (p-E[p]) ] ** 2 / V[t] / V[p] + + where the expectation in numerator is + + E[ (t-E[t]) * (p-E[p]) ] = E[t * p] - E[t] * E[p] + + Different kernel corresponds to different weights. + + For now, y_true and y_pred have to be at least 4d tensor, including batch axis. + + Reference: + + - Zero-normalized cross-correlation (ZNCC): + https://en.wikipedia.org/wiki/Cross-correlation + - Code: https://github.com/voxelmorph/voxelmorph/blob/legacy/src/losses.py + """ + + kernel_fn_dict = dict( + gaussian=gaussian_kernel1d, + rectangular=rectangular_kernel1d, + triangular=triangular_kernel1d, + ) + + def __init__( + self, + kernel_size: int = 9, + kernel_type: str = "rectangular", + # reduction: str = tf.keras.losses.Reduction.SUM, + name: str = "LocalNormalizedCrossCorrelationDEBUG1", + ): + """ + Init. + + :param kernel_size: int. Kernel size or kernel sigma for kernel_type='gauss'. + :param kernel_type: str, rectangular, triangular or gaussian + :param name: name of the loss + """ + super().__init__(name=name) + if kernel_type not in self.kernel_fn_dict.keys(): + raise ValueError( + f"Wrong kernel_type {kernel_type} for LNCC loss type. " + f"Feasible values are {self.kernel_fn_dict.keys()}" + ) + self.kernel_fn = self.kernel_fn_dict[kernel_type] + self.kernel_type = kernel_type + self.kernel_size = kernel_size + + # (kernel_size, ) + self.kernel = self.kernel_fn(kernel_size=self.kernel_size) + # E[1] = sum_i(w_i), () + self.kernel_vol = tf.reduce_sum( + self.kernel[:, None, None] + * self.kernel[None, :, None] + * self.kernel[None, None, :] ) - return tf.reduce_mean(ncc, axis=[1, 2, 3, 4]) + def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: + """ + Return loss for a batch. + + :param y_true: shape = (batch, dim1, dim2, dim3) + or (batch, dim1, dim2, dim3, ch) + :param y_pred: shape = (batch, dim1, dim2, dim3) + or (batch, dim1, dim2, dim3, ch) + :return: shape = (batch,) + """ + # adjust + if len(y_true.shape) == 4: + y_true = tf.expand_dims(y_true, axis=4) + y_pred = tf.expand_dims(y_pred, axis=4) + assert len(y_true.shape) == len(y_pred.shape) == 5 + + # t = y_true, p = y_pred + # (batch, dim1, dim2, dim3, ch) + t2 = y_true * y_true + p2 = y_pred * y_pred + tp = y_true * y_pred + + # sum over kernel + # (batch, dim1, dim2, dim3, 1) + t_sum = separable_filter(y_true, kernel=self.kernel) # E[t] * E[1] + p_sum = separable_filter(y_pred, kernel=self.kernel) # E[p] * E[1] + t2_sum = separable_filter(t2, kernel=self.kernel) # E[tt] * E[1] + p2_sum = separable_filter(p2, kernel=self.kernel) # E[pp] * E[1] + tp_sum = separable_filter(tp, kernel=self.kernel) # E[tp] * E[1] + + # average over kernel + # (batch, dim1, dim2, dim3, 1) + t_avg = t_sum / self.kernel_vol # E[t] + p_avg = p_sum / self.kernel_vol # E[p] + + # shape = (batch, dim1, dim2, dim3, 1) + cross = tp_sum - p_avg * t_sum # E[tp] * E[1] - E[p] * E[t] * E[1] + t_var = t2_sum - t_avg * t_sum # V[t] * E[1] + p_var = p2_sum - p_avg * p_sum # V[p] * E[1] + + # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] + num = cross * cross + EPS + denom = t_var * p_var + EPS + + num = tf.debugging.check_numerics( + num, "DEBUG LNCC num value NAN/INF", name="DEBUGLNCC_before_mean_num" + ) + + return denom def get_config(self) -> dict: """Return the config dictionary for recreating this class.""" diff --git a/deepreg/model/network.py b/deepreg/model/network.py index a06a30c5e..a9765df88 100644 --- a/deepreg/model/network.py +++ b/deepreg/model/network.py @@ -6,7 +6,10 @@ import tensorflow as tf -from deepreg.loss.image import LocalNormalizedCrossCorrelationDEBUG +from deepreg.loss.image import ( + LocalNormalizedCrossCorrelationDEBUG1, + LocalNormalizedCrossCorrelationDEBUG2, +) from deepreg.loss.label import DiceScore, compute_centroid_distance from deepreg.model import layer, layer_util from deepreg.model.backbone import GlobalNet @@ -254,10 +257,14 @@ def build_loss(self): # image loss, conditional model does not have this if "pred_fixed_image" in self._outputs: pred_fixed_image = self._outputs["pred_fixed_image"] - debug_value = LocalNormalizedCrossCorrelationDEBUG()( + debug_value1 = LocalNormalizedCrossCorrelationDEBUG1()( y_true=fixed_image, y_pred=pred_fixed_image ) - self.log_tensor_stats(debug_value, name="debug-lncc") + self.log_tensor_stats(debug_value1, name="debug-lncc-num") + debug_value2 = LocalNormalizedCrossCorrelationDEBUG2()( + y_true=fixed_image, y_pred=pred_fixed_image + ) + self.log_tensor_stats(debug_value2, name="debug-lncc-num") self._build_loss( name="image", inputs_dict=dict(y_true=fixed_image, y_pred=pred_fixed_image), From b764112f4616c69e8054a68ddbbec1f6d164fa60 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 03:10:27 +0000 Subject: [PATCH 30/45] Issue #690: fix typo --- deepreg/model/network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deepreg/model/network.py b/deepreg/model/network.py index a9765df88..048fc1067 100644 --- a/deepreg/model/network.py +++ b/deepreg/model/network.py @@ -264,7 +264,7 @@ def build_loss(self): debug_value2 = LocalNormalizedCrossCorrelationDEBUG2()( y_true=fixed_image, y_pred=pred_fixed_image ) - self.log_tensor_stats(debug_value2, name="debug-lncc-num") + self.log_tensor_stats(debug_value2, name="debug-lncc-denom") self._build_loss( name="image", inputs_dict=dict(y_true=fixed_image, y_pred=pred_fixed_image), From 35beea783a7775a491275132c8133bfd898f2c47 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 03:11:26 +0000 Subject: [PATCH 31/45] Issue #690: fix typo --- deepreg/loss/image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 6e63d731d..15f11860b 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -455,7 +455,7 @@ def __init__( kernel_size: int = 9, kernel_type: str = "rectangular", # reduction: str = tf.keras.losses.Reduction.SUM, - name: str = "LocalNormalizedCrossCorrelationDEBUG1", + name: str = "LocalNormalizedCrossCorrelationDEBUG2", ): """ Init. From 4c47097dea9b69497fddf0d2d915e20217f00a5b Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 03:22:53 +0000 Subject: [PATCH 32/45] Issue #690: do not use separable filter --- deepreg/loss/image.py | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 15f11860b..45f6c959e 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -472,7 +472,7 @@ def __init__( ) self.kernel_fn = self.kernel_fn_dict[kernel_type] self.kernel_type = kernel_type - self.kernel_size = kernel_size + self.kernel_size = kernel_size ** 3 # (kernel_size, ) self.kernel = self.kernel_fn(kernel_size=self.kernel_size) @@ -499,6 +499,15 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: y_pred = tf.expand_dims(y_pred, axis=4) assert len(y_true.shape) == len(y_pred.shape) == 5 + # [dim1, dim2, dim3, d_in, d_out] + # ch must be evenly divisible by d_in + ch = y_true.shape[4] + filters = tf.ones( + shape=[self.kernel_size, self.kernel_size, self.kernel_size, ch, 1] + ) + strides = [1, 1, 1, 1, 1] + padding = "SAME" + # t = y_true, p = y_pred # (batch, dim1, dim2, dim3, ch) t2 = y_true * y_true @@ -507,11 +516,16 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # sum over kernel # (batch, dim1, dim2, dim3, 1) - t_sum = separable_filter(y_true, kernel=self.kernel) # E[t] * E[1] - p_sum = separable_filter(y_pred, kernel=self.kernel) # E[p] * E[1] - t2_sum = separable_filter(t2, kernel=self.kernel) # E[tt] * E[1] - p2_sum = separable_filter(p2, kernel=self.kernel) # E[pp] * E[1] - tp_sum = separable_filter(tp, kernel=self.kernel) # E[tp] * E[1] + # t_sum = separable_filter(y_true, kernel=self.kernel) # E[t] * E[1] + # p_sum = separable_filter(y_pred, kernel=self.kernel) # E[p] * E[1] + # t2_sum = separable_filter(t2, kernel=self.kernel) # E[tt] * E[1] + # p2_sum = separable_filter(p2, kernel=self.kernel) # E[pp] * E[1] + # tp_sum = separable_filter(tp, kernel=self.kernel) # E[tp] * E[1] + t_sum = tf.nn.conv3d(y_true, filters=filters, strides=strides, padding=padding) + p_sum = tf.nn.conv3d(y_pred, filters=filters, strides=strides, padding=padding) + t2_sum = tf.nn.conv3d(t2, filters=filters, strides=strides, padding=padding) + p2_sum = tf.nn.conv3d(p2, filters=filters, strides=strides, padding=padding) + tp_sum = tf.nn.conv3d(tp, filters=filters, strides=strides, padding=padding) # average over kernel # (batch, dim1, dim2, dim3, 1) From 357837accce7298ce797d5c9acfbf64ef5c2797d Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 03:33:36 +0000 Subject: [PATCH 33/45] Issue #690: correct lncc --- deepreg/loss/image.py | 122 ++++++++++++++++++++++-------------------- 1 file changed, 64 insertions(+), 58 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 45f6c959e..90f8202fe 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -3,11 +3,7 @@ from deepreg.loss.util import NegativeLossMixin from deepreg.loss.util import gaussian_kernel1d_size as gaussian_kernel1d -from deepreg.loss.util import ( - rectangular_kernel1d, - separable_filter, - triangular_kernel1d, -) +from deepreg.loss.util import rectangular_kernel1d, triangular_kernel1d from deepreg.registry import REGISTRY EPS = 1.0e-5 @@ -209,16 +205,14 @@ def __init__( ) self.kernel_fn = self.kernel_fn_dict[kernel_type] self.kernel_type = kernel_type - self.kernel_size = kernel_size + self.filters = tf.ones( + shape=[self.kernel_size, self.kernel_size, self.kernel_size, 1, 1] + ) + self.strides = [1, 1, 1, 1, 1] + self.padding = "SAME" - # (kernel_size, ) - self.kernel = self.kernel_fn(kernel_size=self.kernel_size) # E[1] = sum_i(w_i), () - self.kernel_vol = tf.reduce_sum( - self.kernel[:, None, None] - * self.kernel[None, :, None] - * self.kernel[None, None, :] - ) + self.kernel_vol = kernel_size ** 3 def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: """ @@ -244,11 +238,21 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # sum over kernel # (batch, dim1, dim2, dim3, 1) - t_sum = separable_filter(y_true, kernel=self.kernel) # E[t] * E[1] - p_sum = separable_filter(y_pred, kernel=self.kernel) # E[p] * E[1] - t2_sum = separable_filter(t2, kernel=self.kernel) # E[tt] * E[1] - p2_sum = separable_filter(p2, kernel=self.kernel) # E[pp] * E[1] - tp_sum = separable_filter(tp, kernel=self.kernel) # E[tp] * E[1] + t_sum = tf.nn.conv3d( + y_true, filters=self.filters, strides=self.strides, padding=self.padding + ) + p_sum = tf.nn.conv3d( + y_pred, filters=self.filters, strides=self.strides, padding=self.padding + ) + t2_sum = tf.nn.conv3d( + t2, filters=self.filters, strides=self.strides, padding=self.padding + ) + p2_sum = tf.nn.conv3d( + p2, filters=self.filters, strides=self.strides, padding=self.padding + ) + tp_sum = tf.nn.conv3d( + tp, filters=self.filters, strides=self.strides, padding=self.padding + ) # average over kernel # (batch, dim1, dim2, dim3, 1) @@ -341,16 +345,14 @@ def __init__( ) self.kernel_fn = self.kernel_fn_dict[kernel_type] self.kernel_type = kernel_type - self.kernel_size = kernel_size + self.filters = tf.ones( + shape=[self.kernel_size, self.kernel_size, self.kernel_size, 1, 1] + ) + self.strides = [1, 1, 1, 1, 1] + self.padding = "SAME" - # (kernel_size, ) - self.kernel = self.kernel_fn(kernel_size=self.kernel_size) # E[1] = sum_i(w_i), () - self.kernel_vol = tf.reduce_sum( - self.kernel[:, None, None] - * self.kernel[None, :, None] - * self.kernel[None, None, :] - ) + self.kernel_vol = kernel_size ** 3 def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: """ @@ -376,11 +378,21 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # sum over kernel # (batch, dim1, dim2, dim3, 1) - t_sum = separable_filter(y_true, kernel=self.kernel) # E[t] * E[1] - p_sum = separable_filter(y_pred, kernel=self.kernel) # E[p] * E[1] - t2_sum = separable_filter(t2, kernel=self.kernel) # E[tt] * E[1] - p2_sum = separable_filter(p2, kernel=self.kernel) # E[pp] * E[1] - tp_sum = separable_filter(tp, kernel=self.kernel) # E[tp] * E[1] + t_sum = tf.nn.conv3d( + y_true, filters=self.filters, strides=self.strides, padding=self.padding + ) + p_sum = tf.nn.conv3d( + y_pred, filters=self.filters, strides=self.strides, padding=self.padding + ) + t2_sum = tf.nn.conv3d( + t2, filters=self.filters, strides=self.strides, padding=self.padding + ) + p2_sum = tf.nn.conv3d( + p2, filters=self.filters, strides=self.strides, padding=self.padding + ) + tp_sum = tf.nn.conv3d( + tp, filters=self.filters, strides=self.strides, padding=self.padding + ) # average over kernel # (batch, dim1, dim2, dim3, 1) @@ -472,16 +484,14 @@ def __init__( ) self.kernel_fn = self.kernel_fn_dict[kernel_type] self.kernel_type = kernel_type - self.kernel_size = kernel_size ** 3 + self.filters = tf.ones( + shape=[self.kernel_size, self.kernel_size, self.kernel_size, 1, 1] + ) + self.strides = [1, 1, 1, 1, 1] + self.padding = "SAME" - # (kernel_size, ) - self.kernel = self.kernel_fn(kernel_size=self.kernel_size) # E[1] = sum_i(w_i), () - self.kernel_vol = tf.reduce_sum( - self.kernel[:, None, None] - * self.kernel[None, :, None] - * self.kernel[None, None, :] - ) + self.kernel_vol = kernel_size ** 3 def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: """ @@ -499,15 +509,6 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: y_pred = tf.expand_dims(y_pred, axis=4) assert len(y_true.shape) == len(y_pred.shape) == 5 - # [dim1, dim2, dim3, d_in, d_out] - # ch must be evenly divisible by d_in - ch = y_true.shape[4] - filters = tf.ones( - shape=[self.kernel_size, self.kernel_size, self.kernel_size, ch, 1] - ) - strides = [1, 1, 1, 1, 1] - padding = "SAME" - # t = y_true, p = y_pred # (batch, dim1, dim2, dim3, ch) t2 = y_true * y_true @@ -516,16 +517,21 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # sum over kernel # (batch, dim1, dim2, dim3, 1) - # t_sum = separable_filter(y_true, kernel=self.kernel) # E[t] * E[1] - # p_sum = separable_filter(y_pred, kernel=self.kernel) # E[p] * E[1] - # t2_sum = separable_filter(t2, kernel=self.kernel) # E[tt] * E[1] - # p2_sum = separable_filter(p2, kernel=self.kernel) # E[pp] * E[1] - # tp_sum = separable_filter(tp, kernel=self.kernel) # E[tp] * E[1] - t_sum = tf.nn.conv3d(y_true, filters=filters, strides=strides, padding=padding) - p_sum = tf.nn.conv3d(y_pred, filters=filters, strides=strides, padding=padding) - t2_sum = tf.nn.conv3d(t2, filters=filters, strides=strides, padding=padding) - p2_sum = tf.nn.conv3d(p2, filters=filters, strides=strides, padding=padding) - tp_sum = tf.nn.conv3d(tp, filters=filters, strides=strides, padding=padding) + t_sum = tf.nn.conv3d( + y_true, filters=self.filters, strides=self.strides, padding=self.padding + ) + p_sum = tf.nn.conv3d( + y_pred, filters=self.filters, strides=self.strides, padding=self.padding + ) + t2_sum = tf.nn.conv3d( + t2, filters=self.filters, strides=self.strides, padding=self.padding + ) + p2_sum = tf.nn.conv3d( + p2, filters=self.filters, strides=self.strides, padding=self.padding + ) + tp_sum = tf.nn.conv3d( + tp, filters=self.filters, strides=self.strides, padding=self.padding + ) # average over kernel # (batch, dim1, dim2, dim3, 1) From 23645a08751b653b8501b29e3f47fed60d872d5e Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 03:36:35 +0000 Subject: [PATCH 34/45] Issue #690: fix bug --- deepreg/loss/image.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 90f8202fe..f48a9c028 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -205,6 +205,7 @@ def __init__( ) self.kernel_fn = self.kernel_fn_dict[kernel_type] self.kernel_type = kernel_type + self.kernel_size = kernel_size self.filters = tf.ones( shape=[self.kernel_size, self.kernel_size, self.kernel_size, 1, 1] ) @@ -345,6 +346,7 @@ def __init__( ) self.kernel_fn = self.kernel_fn_dict[kernel_type] self.kernel_type = kernel_type + self.kernel_size = kernel_size self.filters = tf.ones( shape=[self.kernel_size, self.kernel_size, self.kernel_size, 1, 1] ) @@ -484,6 +486,7 @@ def __init__( ) self.kernel_fn = self.kernel_fn_dict[kernel_type] self.kernel_type = kernel_type + self.kernel_size = kernel_size self.filters = tf.ones( shape=[self.kernel_size, self.kernel_size, self.kernel_size, 1, 1] ) From 01ba18a45c1ad96b338ec8b1eb0a9e54624971f5 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 19:06:43 +0000 Subject: [PATCH 35/45] Issue #690: divide filters --- deepreg/loss/image.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index f48a9c028..9780fdde0 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -487,8 +487,9 @@ def __init__( self.kernel_fn = self.kernel_fn_dict[kernel_type] self.kernel_type = kernel_type self.kernel_size = kernel_size - self.filters = tf.ones( - shape=[self.kernel_size, self.kernel_size, self.kernel_size, 1, 1] + self.filters = ( + tf.ones(shape=[self.kernel_size, self.kernel_size, self.kernel_size, 1, 1]) + / kernel_size ** 3 ) self.strides = [1, 1, 1, 1, 1] self.padding = "SAME" @@ -538,13 +539,16 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: # average over kernel # (batch, dim1, dim2, dim3, 1) - t_avg = t_sum / self.kernel_vol # E[t] - p_avg = p_sum / self.kernel_vol # E[p] + # t_avg = t_sum / self.kernel_vol # E[t] + # p_avg = p_sum / self.kernel_vol # E[p] # shape = (batch, dim1, dim2, dim3, 1) - cross = tp_sum - p_avg * t_sum # E[tp] * E[1] - E[p] * E[t] * E[1] - t_var = t2_sum - t_avg * t_sum # V[t] * E[1] - p_var = p2_sum - p_avg * p_sum # V[p] * E[1] + # cross = tp_sum - p_avg * t_sum # E[tp] * E[1] - E[p] * E[t] * E[1] + # t_var = t2_sum - t_avg * t_sum # V[t] * E[1] + # p_var = p2_sum - p_avg * p_sum # V[p] * E[1] + cross = tp_sum - p_sum * t_sum # E[tp] * E[1] - E[p] * E[t] * E[1] + t_var = t2_sum - t_sum * t_sum # V[t] * E[1] + p_var = p2_sum - p_sum * p_sum # V[p] * E[1] # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] num = cross * cross + EPS From d1c63bf935fc9373d097ef7103860a771f5d5048 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 19:18:09 +0000 Subject: [PATCH 36/45] Issue #690: simplify code --- deepreg/loss/image.py | 286 ++------------------------------------- deepreg/model/network.py | 14 +- 2 files changed, 15 insertions(+), 285 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 9780fdde0..de97325b7 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -1,4 +1,6 @@ """Provide different loss or metrics classes for images.""" +from typing import Tuple + import tensorflow as tf from deepreg.loss.util import NegativeLossMixin @@ -225,49 +227,9 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: or (batch, dim1, dim2, dim3, ch) :return: shape = (batch,) """ - # adjust - if len(y_true.shape) == 4: - y_true = tf.expand_dims(y_true, axis=4) - y_pred = tf.expand_dims(y_pred, axis=4) - assert len(y_true.shape) == len(y_pred.shape) == 5 - - # t = y_true, p = y_pred - # (batch, dim1, dim2, dim3, ch) - t2 = y_true * y_true - p2 = y_pred * y_pred - tp = y_true * y_pred - - # sum over kernel - # (batch, dim1, dim2, dim3, 1) - t_sum = tf.nn.conv3d( - y_true, filters=self.filters, strides=self.strides, padding=self.padding - ) - p_sum = tf.nn.conv3d( - y_pred, filters=self.filters, strides=self.strides, padding=self.padding - ) - t2_sum = tf.nn.conv3d( - t2, filters=self.filters, strides=self.strides, padding=self.padding - ) - p2_sum = tf.nn.conv3d( - p2, filters=self.filters, strides=self.strides, padding=self.padding - ) - tp_sum = tf.nn.conv3d( - tp, filters=self.filters, strides=self.strides, padding=self.padding - ) - - # average over kernel - # (batch, dim1, dim2, dim3, 1) - t_avg = t_sum / self.kernel_vol # E[t] - p_avg = p_sum / self.kernel_vol # E[p] - - # shape = (batch, dim1, dim2, dim3, 1) - cross = tp_sum - p_avg * t_sum # E[tp] * E[1] - E[p] * E[t] * E[1] - t_var = t2_sum - t_avg * t_sum # V[t] * E[1] - p_var = p2_sum - p_avg * p_sum # V[p] * E[1] - - # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] - num = cross * cross + EPS - denom = t_var * p_var + EPS + num, denom = self._call(y_true, y_pred) + num = num + EPS + denom = denom + EPS ncc = num / denom ncc = tf.debugging.check_numerics( @@ -276,87 +238,9 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: return tf.reduce_mean(ncc, axis=[1, 2, 3, 4]) - def get_config(self) -> dict: - """Return the config dictionary for recreating this class.""" - config = super().get_config() - config["kernel_size"] = self.kernel_size - config["kernel_type"] = self.kernel_type - return config - - -class LocalNormalizedCrossCorrelationDEBUG1(tf.keras.layers.Layer): - """ - Local squared zero-normalized cross-correlation. - - Denote y_true as t and y_pred as p. Consider a window having n elements. - Each position in the window corresponds a weight w_i for i=1:n. - - Define the discrete expectation in the window E[t] as - - E[t] = sum_i(w_i * t_i) / sum_i(w_i) - - Similarly, the discrete variance in the window V[t] is - - V[t] = E[t**2] - E[t] ** 2 - - The local squared zero-normalized cross-correlation is therefore - - E[ (t-E[t]) * (p-E[p]) ] ** 2 / V[t] / V[p] - - where the expectation in numerator is - - E[ (t-E[t]) * (p-E[p]) ] = E[t * p] - E[t] * E[p] - - Different kernel corresponds to different weights. - - For now, y_true and y_pred have to be at least 4d tensor, including batch axis. - - Reference: - - - Zero-normalized cross-correlation (ZNCC): - https://en.wikipedia.org/wiki/Cross-correlation - - Code: https://github.com/voxelmorph/voxelmorph/blob/legacy/src/losses.py - """ - - kernel_fn_dict = dict( - gaussian=gaussian_kernel1d, - rectangular=rectangular_kernel1d, - triangular=triangular_kernel1d, - ) - - def __init__( - self, - kernel_size: int = 9, - kernel_type: str = "rectangular", - # reduction: str = tf.keras.losses.Reduction.SUM, - name: str = "LocalNormalizedCrossCorrelationDEBUG1", - ): - """ - Init. - - :param kernel_size: int. Kernel size or kernel sigma for kernel_type='gauss'. - :param kernel_type: str, rectangular, triangular or gaussian - :param name: name of the loss - """ - super().__init__(name=name) - if kernel_type not in self.kernel_fn_dict.keys(): - raise ValueError( - f"Wrong kernel_type {kernel_type} for LNCC loss type. " - f"Feasible values are {self.kernel_fn_dict.keys()}" - ) - self.kernel_fn = self.kernel_fn_dict[kernel_type] - self.kernel_type = kernel_type - self.kernel_size = kernel_size - self.filters = tf.ones( - shape=[self.kernel_size, self.kernel_size, self.kernel_size, 1, 1] - ) - self.strides = [1, 1, 1, 1, 1] - self.padding = "SAME" - - # E[1] = sum_i(w_i), () - self.kernel_vol = kernel_size ** 3 - - def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: + def _call( + self, y_true: tf.Tensor, y_pred: tf.Tensor + ) -> Tuple[tf.Tensor, tf.Tensor]: """ Return loss for a batch. @@ -407,158 +291,10 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: p_var = p2_sum - p_avg * p_sum # V[p] * E[1] # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] - num = cross * cross + EPS - denom = t_var * p_var + EPS - - denom = tf.debugging.check_numerics( - denom, "DEBUG LNCC denom value NAN/INF", name="DEBUGLNCC_before_mean_denom" - ) - - return num - - def get_config(self) -> dict: - """Return the config dictionary for recreating this class.""" - config = super().get_config() - config["kernel_size"] = self.kernel_size - config["kernel_type"] = self.kernel_type - return config - - -class LocalNormalizedCrossCorrelationDEBUG2(tf.keras.layers.Layer): - """ - Local squared zero-normalized cross-correlation. - - Denote y_true as t and y_pred as p. Consider a window having n elements. - Each position in the window corresponds a weight w_i for i=1:n. - - Define the discrete expectation in the window E[t] as - - E[t] = sum_i(w_i * t_i) / sum_i(w_i) - - Similarly, the discrete variance in the window V[t] is - - V[t] = E[t**2] - E[t] ** 2 - - The local squared zero-normalized cross-correlation is therefore - - E[ (t-E[t]) * (p-E[p]) ] ** 2 / V[t] / V[p] - - where the expectation in numerator is - - E[ (t-E[t]) * (p-E[p]) ] = E[t * p] - E[t] * E[p] - - Different kernel corresponds to different weights. - - For now, y_true and y_pred have to be at least 4d tensor, including batch axis. - - Reference: - - - Zero-normalized cross-correlation (ZNCC): - https://en.wikipedia.org/wiki/Cross-correlation - - Code: https://github.com/voxelmorph/voxelmorph/blob/legacy/src/losses.py - """ - - kernel_fn_dict = dict( - gaussian=gaussian_kernel1d, - rectangular=rectangular_kernel1d, - triangular=triangular_kernel1d, - ) - - def __init__( - self, - kernel_size: int = 9, - kernel_type: str = "rectangular", - # reduction: str = tf.keras.losses.Reduction.SUM, - name: str = "LocalNormalizedCrossCorrelationDEBUG2", - ): - """ - Init. - - :param kernel_size: int. Kernel size or kernel sigma for kernel_type='gauss'. - :param kernel_type: str, rectangular, triangular or gaussian - :param name: name of the loss - """ - super().__init__(name=name) - if kernel_type not in self.kernel_fn_dict.keys(): - raise ValueError( - f"Wrong kernel_type {kernel_type} for LNCC loss type. " - f"Feasible values are {self.kernel_fn_dict.keys()}" - ) - self.kernel_fn = self.kernel_fn_dict[kernel_type] - self.kernel_type = kernel_type - self.kernel_size = kernel_size - self.filters = ( - tf.ones(shape=[self.kernel_size, self.kernel_size, self.kernel_size, 1, 1]) - / kernel_size ** 3 - ) - self.strides = [1, 1, 1, 1, 1] - self.padding = "SAME" - - # E[1] = sum_i(w_i), () - self.kernel_vol = kernel_size ** 3 - - def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: - """ - Return loss for a batch. - - :param y_true: shape = (batch, dim1, dim2, dim3) - or (batch, dim1, dim2, dim3, ch) - :param y_pred: shape = (batch, dim1, dim2, dim3) - or (batch, dim1, dim2, dim3, ch) - :return: shape = (batch,) - """ - # adjust - if len(y_true.shape) == 4: - y_true = tf.expand_dims(y_true, axis=4) - y_pred = tf.expand_dims(y_pred, axis=4) - assert len(y_true.shape) == len(y_pred.shape) == 5 - - # t = y_true, p = y_pred - # (batch, dim1, dim2, dim3, ch) - t2 = y_true * y_true - p2 = y_pred * y_pred - tp = y_true * y_pred - - # sum over kernel - # (batch, dim1, dim2, dim3, 1) - t_sum = tf.nn.conv3d( - y_true, filters=self.filters, strides=self.strides, padding=self.padding - ) - p_sum = tf.nn.conv3d( - y_pred, filters=self.filters, strides=self.strides, padding=self.padding - ) - t2_sum = tf.nn.conv3d( - t2, filters=self.filters, strides=self.strides, padding=self.padding - ) - p2_sum = tf.nn.conv3d( - p2, filters=self.filters, strides=self.strides, padding=self.padding - ) - tp_sum = tf.nn.conv3d( - tp, filters=self.filters, strides=self.strides, padding=self.padding - ) - - # average over kernel - # (batch, dim1, dim2, dim3, 1) - # t_avg = t_sum / self.kernel_vol # E[t] - # p_avg = p_sum / self.kernel_vol # E[p] - - # shape = (batch, dim1, dim2, dim3, 1) - # cross = tp_sum - p_avg * t_sum # E[tp] * E[1] - E[p] * E[t] * E[1] - # t_var = t2_sum - t_avg * t_sum # V[t] * E[1] - # p_var = p2_sum - p_avg * p_sum # V[p] * E[1] - cross = tp_sum - p_sum * t_sum # E[tp] * E[1] - E[p] * E[t] * E[1] - t_var = t2_sum - t_sum * t_sum # V[t] * E[1] - p_var = p2_sum - p_sum * p_sum # V[p] * E[1] - - # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] - num = cross * cross + EPS - denom = t_var * p_var + EPS - - num = tf.debugging.check_numerics( - num, "DEBUG LNCC num value NAN/INF", name="DEBUGLNCC_before_mean_num" - ) + num = cross * cross + denom = t_var * p_var - return denom + return num, denom def get_config(self) -> dict: """Return the config dictionary for recreating this class.""" diff --git a/deepreg/model/network.py b/deepreg/model/network.py index 048fc1067..271e64b77 100644 --- a/deepreg/model/network.py +++ b/deepreg/model/network.py @@ -6,10 +6,7 @@ import tensorflow as tf -from deepreg.loss.image import ( - LocalNormalizedCrossCorrelationDEBUG1, - LocalNormalizedCrossCorrelationDEBUG2, -) +from deepreg.loss.image import LocalNormalizedCrossCorrelation from deepreg.loss.label import DiceScore, compute_centroid_distance from deepreg.model import layer, layer_util from deepreg.model.backbone import GlobalNet @@ -257,14 +254,11 @@ def build_loss(self): # image loss, conditional model does not have this if "pred_fixed_image" in self._outputs: pred_fixed_image = self._outputs["pred_fixed_image"] - debug_value1 = LocalNormalizedCrossCorrelationDEBUG1()( + num, denom = LocalNormalizedCrossCorrelation()._call( y_true=fixed_image, y_pred=pred_fixed_image ) - self.log_tensor_stats(debug_value1, name="debug-lncc-num") - debug_value2 = LocalNormalizedCrossCorrelationDEBUG2()( - y_true=fixed_image, y_pred=pred_fixed_image - ) - self.log_tensor_stats(debug_value2, name="debug-lncc-denom") + self.log_tensor_stats(num, name="debug-lncc-num") + self.log_tensor_stats(denom, name="debug-lncc-denom") self._build_loss( name="image", inputs_dict=dict(y_true=fixed_image, y_pred=pred_fixed_image), From fc98187a5592dba7a95479a3570962fe57971ec0 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 19:20:54 +0000 Subject: [PATCH 37/45] Issue #690: modify filters --- deepreg/loss/image.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index de97325b7..e9b217f5c 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -210,12 +210,12 @@ def __init__( self.kernel_size = kernel_size self.filters = tf.ones( shape=[self.kernel_size, self.kernel_size, self.kernel_size, 1, 1] - ) + ) / (kernel_size ** 3) self.strides = [1, 1, 1, 1, 1] self.padding = "SAME" # E[1] = sum_i(w_i), () - self.kernel_vol = kernel_size ** 3 + # self.kernel_vol = kernel_size ** 3 def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: """ @@ -282,13 +282,13 @@ def _call( # average over kernel # (batch, dim1, dim2, dim3, 1) - t_avg = t_sum / self.kernel_vol # E[t] - p_avg = p_sum / self.kernel_vol # E[p] + # t_avg = t_sum / self.kernel_vol # E[t] + # p_avg = p_sum / self.kernel_vol # E[p] # shape = (batch, dim1, dim2, dim3, 1) - cross = tp_sum - p_avg * t_sum # E[tp] * E[1] - E[p] * E[t] * E[1] - t_var = t2_sum - t_avg * t_sum # V[t] * E[1] - p_var = p2_sum - p_avg * p_sum # V[p] * E[1] + cross = tp_sum - p_sum * t_sum # E[tp] * E[1] - E[p] * E[t] * E[1] + t_var = t2_sum - t_sum * t_sum # V[t] * E[1] + p_var = p2_sum - p_sum * p_sum # V[p] * E[1] # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] num = cross * cross From d696cb44b4faa17f4ba681093470b00d3bd9dfb6 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 22:34:43 +0000 Subject: [PATCH 38/45] Issue #690: change back to 1d kernel --- deepreg/loss/image.py | 46 ++++++++++++++----------------------------- deepreg/loss/util.py | 20 +++++++++++++------ 2 files changed, 29 insertions(+), 37 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index e9b217f5c..7cfd751af 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -5,7 +5,11 @@ from deepreg.loss.util import NegativeLossMixin from deepreg.loss.util import gaussian_kernel1d_size as gaussian_kernel1d -from deepreg.loss.util import rectangular_kernel1d, triangular_kernel1d +from deepreg.loss.util import ( + rectangular_kernel1d, + separable_filter, + triangular_kernel1d, +) from deepreg.registry import REGISTRY EPS = 1.0e-5 @@ -208,14 +212,9 @@ def __init__( self.kernel_fn = self.kernel_fn_dict[kernel_type] self.kernel_type = kernel_type self.kernel_size = kernel_size - self.filters = tf.ones( - shape=[self.kernel_size, self.kernel_size, self.kernel_size, 1, 1] - ) / (kernel_size ** 3) - self.strides = [1, 1, 1, 1, 1] - self.padding = "SAME" - # E[1] = sum_i(w_i), () - # self.kernel_vol = kernel_size ** 3 + # (kernel_size, ) + self.kernel = self.kernel_fn(kernel_size=self.kernel_size) def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: """ @@ -264,31 +263,16 @@ def _call( # sum over kernel # (batch, dim1, dim2, dim3, 1) - t_sum = tf.nn.conv3d( - y_true, filters=self.filters, strides=self.strides, padding=self.padding - ) - p_sum = tf.nn.conv3d( - y_pred, filters=self.filters, strides=self.strides, padding=self.padding - ) - t2_sum = tf.nn.conv3d( - t2, filters=self.filters, strides=self.strides, padding=self.padding - ) - p2_sum = tf.nn.conv3d( - p2, filters=self.filters, strides=self.strides, padding=self.padding - ) - tp_sum = tf.nn.conv3d( - tp, filters=self.filters, strides=self.strides, padding=self.padding - ) - - # average over kernel - # (batch, dim1, dim2, dim3, 1) - # t_avg = t_sum / self.kernel_vol # E[t] - # p_avg = p_sum / self.kernel_vol # E[p] + t_sum = separable_filter(y_true, kernel=self.kernel) + p_sum = separable_filter(y_pred, kernel=self.kernel) + t2_sum = separable_filter(t2, kernel=self.kernel) + p2_sum = separable_filter(p2, kernel=self.kernel) + tp_sum = separable_filter(tp, kernel=self.kernel) # shape = (batch, dim1, dim2, dim3, 1) - cross = tp_sum - p_sum * t_sum # E[tp] * E[1] - E[p] * E[t] * E[1] - t_var = t2_sum - t_sum * t_sum # V[t] * E[1] - p_var = p2_sum - p_sum * p_sum # V[p] * E[1] + cross = tp_sum - p_sum * t_sum + t_var = t2_sum - t_sum * t_sum + p_var = p2_sum - p_sum * p_sum # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] num = cross * cross diff --git a/deepreg/loss/util.py b/deepreg/loss/util.py index 909d9d49d..7c7522402 100644 --- a/deepreg/loss/util.py +++ b/deepreg/loss/util.py @@ -32,20 +32,23 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: def rectangular_kernel1d(kernel_size: int) -> tf.Tensor: """ - Return a the 1D filter for separable convolution equivalent to a 3-D rectangular - kernel for LocalNormalizedCrossCorrelation. + Return a the 1D rectangular kernel for LocalNormalizedCrossCorrelation. + + Sum of the weights is 1. :param kernel_size: scalar, size of the 1-D kernel :return: kernel_weights, of shape (kernel_size, ) """ - kernel = tf.ones(shape=(kernel_size,), dtype=tf.float32) + kernel = tf.ones(shape=(kernel_size,), dtype=tf.float32) / float(kernel_size) return kernel def triangular_kernel1d(kernel_size: int) -> tf.Tensor: """ - 1D triangular kernel. + Return a the 1D triangular kernel for LocalNormalizedCrossCorrelation. + + Sum of the weights is 1. Assume kernel_size is odd, it will be a smoothed from a kernel which center part is zero. @@ -73,13 +76,17 @@ def triangular_kernel1d(kernel_size: int) -> tf.Tensor: kernel = tf.nn.conv1d( kernel[None, :, None], filters=filters, stride=[1, 1, 1], padding="SAME" ) + kernel = kernel / tf.reduce_sum(kernel) + return kernel[0, :, 0] def gaussian_kernel1d_size(kernel_size: int) -> tf.Tensor: """ - Return a the 1D filter for separable convolution equivalent to a 3-D Gaussian - kernel for LocalNormalizedCrossCorrelation. + Return a the 1D Gaussian kernel for LocalNormalizedCrossCorrelation. + + Sum of the weights is 1. + :param kernel_size: scalar, size of the 1-D kernel :return: filters, of shape (kernel_size, ) """ @@ -88,6 +95,7 @@ def gaussian_kernel1d_size(kernel_size: int) -> tf.Tensor: grid = tf.range(0, kernel_size, dtype=tf.float32) filters = tf.exp(-tf.square(grid - mean) / (2 * sigma ** 2)) + filters = filters / tf.reduce_sum(filters) return filters From e910bdabcbeec94f8281d98cb8b67f099ce11c9a Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 22:51:31 +0000 Subject: [PATCH 39/45] Issue #690: clip variance and remove debug code --- deepreg/loss/image.py | 32 +++++--------------------------- deepreg/model/network.py | 10 ++-------- 2 files changed, 7 insertions(+), 35 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 7cfd751af..132707d9d 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -1,6 +1,4 @@ """Provide different loss or metrics classes for images.""" -from typing import Tuple - import tensorflow as tf from deepreg.loss.util import NegativeLossMixin @@ -226,30 +224,7 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: or (batch, dim1, dim2, dim3, ch) :return: shape = (batch,) """ - num, denom = self._call(y_true, y_pred) - num = num + EPS - denom = denom + EPS - ncc = num / denom - - ncc = tf.debugging.check_numerics( - ncc, "LNCC ncc value NAN/INF", name="LNCC_before_mean" - ) - - return tf.reduce_mean(ncc, axis=[1, 2, 3, 4]) - - def _call( - self, y_true: tf.Tensor, y_pred: tf.Tensor - ) -> Tuple[tf.Tensor, tf.Tensor]: - """ - Return loss for a batch. - - :param y_true: shape = (batch, dim1, dim2, dim3) - or (batch, dim1, dim2, dim3, ch) - :param y_pred: shape = (batch, dim1, dim2, dim3) - or (batch, dim1, dim2, dim3, ch) - :return: shape = (batch,) - """ - # adjust + # adjust shape if len(y_true.shape) == 4: y_true = tf.expand_dims(y_true, axis=4) y_pred = tf.expand_dims(y_pred, axis=4) @@ -277,8 +252,11 @@ def _call( # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] num = cross * cross denom = t_var * p_var + denom = tf.maximum(denom, 0) # make sure variance >=0 - return num, denom + ncc = (num + EPS) / (denom + EPS) + + return tf.reduce_mean(ncc, axis=[1, 2, 3, 4]) def get_config(self) -> dict: """Return the config dictionary for recreating this class.""" diff --git a/deepreg/model/network.py b/deepreg/model/network.py index 271e64b77..770511a93 100644 --- a/deepreg/model/network.py +++ b/deepreg/model/network.py @@ -6,7 +6,6 @@ import tensorflow as tf -from deepreg.loss.image import LocalNormalizedCrossCorrelation from deepreg.loss.label import DiceScore, compute_centroid_distance from deepreg.model import layer, layer_util from deepreg.model.backbone import GlobalNet @@ -31,8 +30,8 @@ class RegistrationModel(tf.keras.Model): def __init__( self, - moving_image_size: tuple, - fixed_image_size: tuple, + moving_image_size: Tuple, + fixed_image_size: Tuple, index_size: int, labeled: bool, batch_size: int, @@ -254,11 +253,6 @@ def build_loss(self): # image loss, conditional model does not have this if "pred_fixed_image" in self._outputs: pred_fixed_image = self._outputs["pred_fixed_image"] - num, denom = LocalNormalizedCrossCorrelation()._call( - y_true=fixed_image, y_pred=pred_fixed_image - ) - self.log_tensor_stats(num, name="debug-lncc-num") - self.log_tensor_stats(denom, name="debug-lncc-denom") self._build_loss( name="image", inputs_dict=dict(y_true=fixed_image, y_pred=pred_fixed_image), From 3a09c82be567eebada350884e7b8e7b872b958be Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 23:13:08 +0000 Subject: [PATCH 40/45] Issue #690: fix tests and remove debug changes --- deepreg/constant.py | 3 +++ deepreg/loss/deform.py | 1 - deepreg/loss/image.py | 3 +-- deepreg/loss/label.py | 3 +-- deepreg/loss/util.py | 3 --- deepreg/model/layer_util.py | 2 +- deepreg/model/network.py | 11 ++--------- deepreg/train.py | 2 +- test/unit/test_loss_label.py | 13 +++++++------ test/unit/test_loss_util.py | 3 +++ test/unit/util.py | 6 ++++-- 11 files changed, 23 insertions(+), 27 deletions(-) create mode 100644 deepreg/constant.py diff --git a/deepreg/constant.py b/deepreg/constant.py new file mode 100644 index 000000000..72746ee3b --- /dev/null +++ b/deepreg/constant.py @@ -0,0 +1,3 @@ +"""Module defining global constants.""" + +EPS = 1.0e-5 diff --git a/deepreg/loss/deform.py b/deepreg/loss/deform.py index ea804e590..45d99a285 100644 --- a/deepreg/loss/deform.py +++ b/deepreg/loss/deform.py @@ -85,7 +85,6 @@ def call(self, inputs: tf.Tensor, **kwargs) -> tf.Tensor: :return: shape = () """ assert len(inputs.shape) == 5 - tf.debugging.check_numerics(inputs, "GRAIDENT ddf value NAN/INF", name=None) ddf = inputs # first order gradient # (batch, m_dim1-2, m_dim2-2, m_dim3-2, 3) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index 132707d9d..fa0f33727 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -1,6 +1,7 @@ """Provide different loss or metrics classes for images.""" import tensorflow as tf +from deepreg.constant import EPS from deepreg.loss.util import NegativeLossMixin from deepreg.loss.util import gaussian_kernel1d_size as gaussian_kernel1d from deepreg.loss.util import ( @@ -10,8 +11,6 @@ ) from deepreg.registry import REGISTRY -EPS = 1.0e-5 - @REGISTRY.register_loss(name="ssd") class SumSquaredDifference(tf.keras.losses.Loss): diff --git a/deepreg/loss/label.py b/deepreg/loss/label.py index 6b7693e30..37429f51a 100644 --- a/deepreg/loss/label.py +++ b/deepreg/loss/label.py @@ -4,13 +4,12 @@ import tensorflow as tf +from deepreg.constant import EPS from deepreg.loss.util import NegativeLossMixin, cauchy_kernel1d from deepreg.loss.util import gaussian_kernel1d_sigma as gaussian_kernel1d from deepreg.loss.util import separable_filter from deepreg.registry import REGISTRY -EPS = 1.0e-5 - class MultiScaleLoss(tf.keras.losses.Loss): """ diff --git a/deepreg/loss/util.py b/deepreg/loss/util.py index 7c7522402..908d54b7e 100644 --- a/deepreg/loss/util.py +++ b/deepreg/loss/util.py @@ -27,9 +27,6 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: return -super().call(y_true=y_true, y_pred=y_pred) -EPS = 1.0e-5 - - def rectangular_kernel1d(kernel_size: int) -> tf.Tensor: """ Return a the 1D rectangular kernel for LocalNormalizedCrossCorrelation. diff --git a/deepreg/model/layer_util.py b/deepreg/model/layer_util.py index ae74606de..3ceb20024 100644 --- a/deepreg/model/layer_util.py +++ b/deepreg/model/layer_util.py @@ -218,7 +218,7 @@ def resample( vol: tf.Tensor, loc: tf.Tensor, interpolation: str = "linear", - zero_boundary: bool = False, + zero_boundary: bool = True, ) -> tf.Tensor: r""" Sample the volume at given locations. diff --git a/deepreg/model/network.py b/deepreg/model/network.py index 770511a93..a51808ec8 100644 --- a/deepreg/model/network.py +++ b/deepreg/model/network.py @@ -60,7 +60,7 @@ def __init__( self.batch_size = batch_size self.config = config self.num_devices = num_devices - self.global_batch_size = num_devices * batch_size * 1.0 + self.global_batch_size = num_devices * batch_size assert self.global_batch_size > 0 self._inputs = None # save inputs of self._model as dict @@ -219,17 +219,10 @@ def _build_loss(self, name: str, inputs_dict: dict): config=dict_without(d=loss_config, key="weight") ) loss_value = loss_layer(**inputs_dict) / self.global_batch_size - weighted_loss = loss_value * weight if weight != 1 else loss_value + weighted_loss = loss_value * weight # add loss self._model.add_loss(weighted_loss) - - loss_value = tf.debugging.check_numerics( - loss_value, - f"loss {name}_{loss_layer.name} inf/nan", - name=f"op_loss_{name}_{loss_layer.name}", - ) - # add metric self._model.add_metric( loss_value, name=f"loss/{name}_{loss_layer.name}", aggregation="mean" diff --git a/deepreg/train.py b/deepreg/train.py index a2fba6cdb..60aeabac2 100644 --- a/deepreg/train.py +++ b/deepreg/train.py @@ -137,7 +137,7 @@ def train( optimizer = opt.build_optimizer(optimizer_config=config["train"]["optimizer"]) # compile - model.compile(optimizer=optimizer) # , run_eagerly=True) + model.compile(optimizer=optimizer) model.plot_model(output_dir=log_dir) # build callbacks diff --git a/test/unit/test_loss_label.py b/test/unit/test_loss_label.py index ba43f2eca..60007d685 100644 --- a/test/unit/test_loss_label.py +++ b/test/unit/test_loss_label.py @@ -12,6 +12,7 @@ import tensorflow as tf import deepreg.loss.label as label +from deepreg.constant import EPS class TestMultiScaleLoss: @@ -91,11 +92,11 @@ def y_pred(self): @pytest.mark.parametrize( "binary,background_weight,scales,expected", [ - (True, 0.0, None, -np.log(1.0e-7)), - (False, 0.0, None, -0.6 * np.log(0.3)), - (False, 0.2, None, -0.48 * np.log(0.3) - 0.08 * np.log(0.7)), - (False, 0.2, [0, 0], -0.48 * np.log(0.3) - 0.08 * np.log(0.7)), - (False, 0.2, [0, 1], 0.5239637), + (True, 0.0, None, -np.log(EPS)), + (False, 0.0, None, -0.6 * np.log(0.3 + EPS)), + (False, 0.2, None, -0.48 * np.log(0.3 + EPS) - 0.08 * np.log(0.7 + EPS)), + (False, 0.2, [0, 0], -0.48 * np.log(0.3 + EPS) - 0.08 * np.log(0.7 + EPS)), + (False, 0.2, [0, 1], 0.5239465), ], ) def test_call(self, y_true, y_pred, binary, background_weight, scales, expected): @@ -135,7 +136,7 @@ def y_pred(self): (True, None, 0), (False, None, 0.25), (False, [0, 0], 0.25), - (False, [0, 1], 0.17484076), + (False, [0, 1], 0.17485845), ], ) def test_call(self, y_true, y_pred, binary, scales, expected): diff --git a/test/unit/test_loss_util.py b/test/unit/test_loss_util.py index 187760555..20e456ee2 100644 --- a/test/unit/test_loss_util.py +++ b/test/unit/test_loss_util.py @@ -62,6 +62,7 @@ def test_gaussian_kernel1d_size(kernel_size): grid = tf.range(0, kernel_size, dtype=tf.float32) expected = tf.exp(-tf.square(grid - mean) / (2 * sigma ** 2)) + expected = expected / tf.reduce_sum(expected) got = gaussian_kernel1d_size(kernel_size) assert is_equal_tf(got, expected) @@ -75,6 +76,7 @@ def test_rectangular_kernel1d(kernel_size): :return: """ expected = tf.ones(shape=(kernel_size,), dtype=tf.float32) + expected = expected / tf.reduce_sum(expected) got = rectangular_kernel1d(kernel_size) assert is_equal_tf(got, expected) @@ -91,6 +93,7 @@ def test_triangular_kernel1d(kernel_size): for it_k in range(kernel_size // 2): expected[it_k] = it_k + 1 expected[-it_k - 1] = it_k + 1 + expected = expected / tf.reduce_sum(expected) got = triangular_kernel1d(kernel_size) assert is_equal_tf(got, expected) diff --git a/test/unit/util.py b/test/unit/util.py index 02866779a..c40128d69 100644 --- a/test/unit/util.py +++ b/test/unit/util.py @@ -3,9 +3,11 @@ import numpy as np import tensorflow as tf +from deepreg.constant import EPS + def is_equal_np( - x: Union[np.ndarray, List], y: Union[np.ndarray, List], atol: float = 1.0e-7 + x: Union[np.ndarray, List], y: Union[np.ndarray, List], atol: float = EPS ) -> bool: """ Check if two numpy arrays are identical. @@ -23,7 +25,7 @@ def is_equal_np( def is_equal_tf( x: Union[tf.Tensor, np.ndarray, List], y: Union[tf.Tensor, np.ndarray, List], - atol: float = 1.0e-7, + atol: float = EPS, ) -> bool: """ Check if two tf tensors are identical. From c8d62ae11a8d83d56a5c9df51a6dbf504a42ae5b Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 23:18:27 +0000 Subject: [PATCH 41/45] Issue #690: change code so that no need to add test --- deepreg/dataset/loader/interface.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/deepreg/dataset/loader/interface.py b/deepreg/dataset/loader/interface.py index bb169d047..bf5bfb87a 100644 --- a/deepreg/dataset/loader/interface.py +++ b/deepreg/dataset/loader/interface.py @@ -379,14 +379,9 @@ def validate_images_and_labels( for arr, name in zip( [moving_image, fixed_image], ["moving_image", "fixed_image"] ): - if len(arr.shape) != 3: + if len(arr.shape) != 3 or min(arr.shape) <= 0: raise ValueError( - f"Sample {image_indices}'s {name}' shape should be 3D. " - f"Got {arr.shape}." - ) - if min(arr.shape) <= 0: - raise ValueError( - f"Sample {image_indices}'s {name}' shape should be non-empty. " + f"Sample {image_indices}'s {name}' shape should be 3D and non-empty. " f"Got {arr.shape}." ) # when data are labeled From 493300bac2706a9a1345b14cc8f98cc72e5d6e74 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 23:44:47 +0000 Subject: [PATCH 42/45] Issue #690: fix pylint --- deepreg/dataset/loader/interface.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/deepreg/dataset/loader/interface.py b/deepreg/dataset/loader/interface.py index bf5bfb87a..997448195 100644 --- a/deepreg/dataset/loader/interface.py +++ b/deepreg/dataset/loader/interface.py @@ -381,8 +381,8 @@ def validate_images_and_labels( ): if len(arr.shape) != 3 or min(arr.shape) <= 0: raise ValueError( - f"Sample {image_indices}'s {name}' shape should be 3D and non-empty. " - f"Got {arr.shape}." + f"Sample {image_indices}'s {name}' shape should be 3D" + f" and non-empty, got {arr.shape}." ) # when data are labeled if moving_label is not None and fixed_label is not None: From c129e588fe3ec959423b62c9b007c4de515c9206 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Fri, 26 Mar 2021 23:58:35 +0000 Subject: [PATCH 43/45] Issue #690: update changelog and fix test --- CHANGELOG.md | 2 ++ test/unit/test_interface.py | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d0089f17f..720ca54f1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,7 @@ compatible with the updates. ### Changed +- Increased all EPS to 1e-5. - Clarify the suggestion in doc to use all-zero masks for missing labels. - Moved contributor list to a separate page. - Changed `no-test` flag to `full` for demo scripts. @@ -37,6 +38,7 @@ compatible with the updates. ### Fixed +- Fixed LNCC loss regarding INF values. - Removed loss weight checks to be more robust. - Fixed import error under python 3.6. - Fixed the residual module in local net architecture, compatible for previous diff --git a/test/unit/test_interface.py b/test/unit/test_interface.py index 0315938d5..22778f610 100644 --- a/test/unit/test_interface.py +++ b/test/unit/test_interface.py @@ -356,7 +356,7 @@ def mock_sample_index_generator(): fixed_label=None, image_indices=[1], ) - assert "Sample [1]'s moving_image' shape should be 3D. " in str(err_info.value) + assert "Sample [1]'s moving_image' shape should be 3D" in str(err_info.value) with pytest.raises(ValueError) as err_info: generator.validate_images_and_labels( fixed_image=dummy_array, From c2269339d19f91dc6ba718768206b18c98cbe473 Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Sat, 27 Mar 2021 17:11:48 +0000 Subject: [PATCH 44/45] Issue #690: modify LNCC implementation to be more stable --- deepreg/loss/image.py | 44 ++++++++++++++++--------------------------- 1 file changed, 16 insertions(+), 28 deletions(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index fa0f33727..fea29f11f 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -157,16 +157,12 @@ class LocalNormalizedCrossCorrelation(tf.keras.losses.Loss): Similarly, the discrete variance in the window V[t] is - V[t] = E[t**2] - E[t] ** 2 + V[t] = E[(t - E[t])**2] The local squared zero-normalized cross-correlation is therefore E[ (t-E[t]) * (p-E[p]) ] ** 2 / V[t] / V[p] - where the expectation in numerator is - - E[ (t-E[t]) * (p-E[p]) ] = E[t * p] - E[t] * E[p] - Different kernel corresponds to different weights. For now, y_true and y_pred have to be at least 4d tensor, including batch axis. @@ -175,7 +171,7 @@ class LocalNormalizedCrossCorrelation(tf.keras.losses.Loss): - Zero-normalized cross-correlation (ZNCC): https://en.wikipedia.org/wiki/Cross-correlation - - Code: https://github.com/voxelmorph/voxelmorph/blob/legacy/src/losses.py + - https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Reliability_weights """ kernel_fn_dict = dict( @@ -223,36 +219,28 @@ def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: or (batch, dim1, dim2, dim3, ch) :return: shape = (batch,) """ - # adjust shape + # adjust shape to be (batch, dim1, dim2, dim3, ch) if len(y_true.shape) == 4: y_true = tf.expand_dims(y_true, axis=4) y_pred = tf.expand_dims(y_pred, axis=4) assert len(y_true.shape) == len(y_pred.shape) == 5 # t = y_true, p = y_pred - # (batch, dim1, dim2, dim3, ch) - t2 = y_true * y_true - p2 = y_pred * y_pred - tp = y_true * y_pred - - # sum over kernel - # (batch, dim1, dim2, dim3, 1) - t_sum = separable_filter(y_true, kernel=self.kernel) - p_sum = separable_filter(y_pred, kernel=self.kernel) - t2_sum = separable_filter(t2, kernel=self.kernel) - p2_sum = separable_filter(p2, kernel=self.kernel) - tp_sum = separable_filter(tp, kernel=self.kernel) - - # shape = (batch, dim1, dim2, dim3, 1) - cross = tp_sum - p_sum * t_sum - t_var = t2_sum - t_sum * t_sum - p_var = p2_sum - p_sum * p_sum - - # (E[tp] - E[p] * E[t]) ** 2 / V[t] / V[p] + t_mean = separable_filter(y_true, kernel=self.kernel) + p_mean = separable_filter(y_pred, kernel=self.kernel) + + t = y_true - t_mean + p = y_pred - p_mean + + # the variance can be biased but as both num and denom are biased + # it got cancelled + # https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Reliability_weights + cross = separable_filter(t * p, kernel=self.kernel) + t_var = separable_filter(t * t, kernel=self.kernel) + p_var = separable_filter(p * p, kernel=self.kernel) + num = cross * cross denom = t_var * p_var - denom = tf.maximum(denom, 0) # make sure variance >=0 - ncc = (num + EPS) / (denom + EPS) return tf.reduce_mean(ncc, axis=[1, 2, 3, 4]) From ac3a46daad27aa664543e977431381f960dbbcab Mon Sep 17 00:00:00 2001 From: Yunguan Fu Date: Sat, 27 Mar 2021 19:52:59 +0000 Subject: [PATCH 45/45] Issue #690: add reference for LNCC --- deepreg/loss/image.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/deepreg/loss/image.py b/deepreg/loss/image.py index fea29f11f..5300a33e4 100644 --- a/deepreg/loss/image.py +++ b/deepreg/loss/image.py @@ -155,6 +155,8 @@ class LocalNormalizedCrossCorrelation(tf.keras.losses.Loss): E[t] = sum_i(w_i * t_i) / sum_i(w_i) + Here, we assume sum_i(w_i) == 1, means the weights have been normalized. + Similarly, the discrete variance in the window V[t] is V[t] = E[(t - E[t])**2] @@ -163,7 +165,12 @@ class LocalNormalizedCrossCorrelation(tf.keras.losses.Loss): E[ (t-E[t]) * (p-E[p]) ] ** 2 / V[t] / V[p] - Different kernel corresponds to different weights. + When calculating variance, we choose to subtract the mean first then calculte + variance instead of calculating E[t**2] - E[t] ** 2, the reason is that when + E[t**2] and E[t] ** 2 are both very large or very small, the subtraction may + have large rounding error and makes the result inaccurate. Also, it is not + guaranteed that the result >= 0. For more details, please read "Algorithms for + computing the sample variance: Analysis and recommendations." page 1. For now, y_true and y_pred have to be at least 4d tensor, including batch axis. @@ -172,6 +179,9 @@ class LocalNormalizedCrossCorrelation(tf.keras.losses.Loss): - Zero-normalized cross-correlation (ZNCC): https://en.wikipedia.org/wiki/Cross-correlation - https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Reliability_weights + - Chan, Tony F., Gene H. Golub, and Randall J. LeVeque. + "Algorithms for computing the sample variance: Analysis and recommendations." + The American Statistician 37.3 (1983): 242-247. """ kernel_fn_dict = dict( @@ -207,6 +217,7 @@ def __init__( self.kernel_size = kernel_size # (kernel_size, ) + # sum of the kernel weights would be one self.kernel = self.kernel_fn(kernel_size=self.kernel_size) def call(self, y_true: tf.Tensor, y_pred: tf.Tensor) -> tf.Tensor: