diff --git a/lumicks/pylake/kymotracker/detail/denoising.py b/lumicks/pylake/kymotracker/detail/denoising.py index c6a32db8e..b2bee917b 100644 --- a/lumicks/pylake/kymotracker/detail/denoising.py +++ b/lumicks/pylake/kymotracker/detail/denoising.py @@ -289,6 +289,7 @@ def filter_image(self, image, false_detection_rate=0.1, verbose=False): verbose : bool Show extra output """ + image = np.array(image, dtype=float) # integer arrays lose quality when convolved detail_coefficients, remainder = self._calculate_wavelets(image, stabilize=True) if verbose: @@ -314,3 +315,105 @@ def filter_image(self, image, false_detection_rate=0.1, verbose=False): output_image[output_image < 0] = 0 return output_image, significant + + def filter_regularized( + self, + image, + false_detection_rate=0.1, + num_iter=10, + remove_background=True, + verbose=False, + ): + """Filter the image using a regularized wavelet reconstruction. + + This reconstructs the image, but with the additional constraint that the + resulting image has to be positive and sparse (L1 regularization). This regularization + procedure helps, since the IUWT is overcomplete. + + The regularization (shrinkage) is performed by shrinking coefficients towards zero. This + is done via a soft-thresholding procedure where coefficients below a threshold are set to + zero, while values above the threshold are shrunk towards zero. The threshold is then + decreased each step (this is a requirement for the algorithm). Note that at each step, the + significant structures are forced to be maintained. + + As a result, we will get a sparse representation that fulfills positivity, while + keeping all the coefficients that were significant. + + Parameters + ---------- + image : np.ndarray + Image array + false_detection_rate : float + Represents the False Detection Rate when it comes to pixels with significant + signal. The probability of erroneously detecting spots in a spot-free homogeneous + noise is upper bounded by this value. + num_iter : int + Number of iterations to run. + remove_background : bool + Remove the background after approximating the image? This amounts to not adding + the final approximation layer. + verbose : bool + Show some output while it is running. + """ + image = np.array(image, dtype=float) # integer arrays lose quality when convolved + + # Grab the significant wavelet coefficients first. + _, significant = self.filter_image(image, false_detection_rate=false_detection_rate) + + detailed_coeffs, remainder = self._calculate_wavelets(image, stabilize=False) + + def positivity_projector(coefficients): + """Calculate the coefficients that will lead to a fully positive reconstruction + + This corresponds to the operator Qs2 in the paper [1]. + + Parameters + ---------- + coefficients : list of np.ndarray + Wavelet detail coefficients + """ + img = sum(coefficients) + remainder + positive_img = np.clip(img, 0, np.inf) + coefficients, rem = self._calculate_wavelets(positive_img, stabilize=False) + return coefficients + + def significance_enforcer(coefficients): + """Force significant structures to remain. + + This corresponds to the operator Ps in the paper [1]. + + Parameters + ---------- + coefficients : list of np.ndarray + Wavelet detail coefficients + """ + for c, d, sig in zip(coefficients, detailed_coeffs, significant): + c[sig] = d[sig] + return coefficients + + current_coeffs = [np.copy(d) for d in detailed_coeffs] + beta = 1.0 + + for k in range(num_iter): + positive = positivity_projector(current_coeffs) + next_solution = significance_enforcer(positive) + + if verbose: + print(f"Iter: {k}: Beta={beta}") + + def soft_threshold(detail_coeffs, treshold): + """Soft-thresholding in the context of L1 optimization involves stepping towards + zero, and truncating any values below the stepsize to zero""" + detail_coeffs[abs(detail_coeffs) < treshold] = 0.0 + detail_coeffs[detail_coeffs > treshold] -= treshold + detail_coeffs[detail_coeffs < -treshold] += treshold + return detail_coeffs + + current_coeffs = [soft_threshold(d.copy(), beta) for d in next_solution] + beta -= 1.0 / num_iter + + reconstructed_img = sum(current_coeffs) + if not remove_background: + reconstructed_img += remainder + + return np.clip(reconstructed_img, 0.0, np.inf) diff --git a/lumicks/pylake/kymotracker/tests/reference_data/denoising_basic/output_image_1d.npz b/lumicks/pylake/kymotracker/tests/reference_data/denoising_basic/output_image_1d.npz index 844d3b67b..e299aebbc 100644 --- a/lumicks/pylake/kymotracker/tests/reference_data/denoising_basic/output_image_1d.npz +++ b/lumicks/pylake/kymotracker/tests/reference_data/denoising_basic/output_image_1d.npz @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:c9374d8e03b7508f3fb927cfd594fc1bf72f86f6bb3810d1e61537b9410e5aac +oid sha256:c2fe0f933f9dbbd615124a34f4740a81b5049755ab7623b06f28867ef375a00d size 13712 diff --git a/lumicks/pylake/kymotracker/tests/reference_data/denoising_basic/output_image_2d.npz b/lumicks/pylake/kymotracker/tests/reference_data/denoising_basic/output_image_2d.npz index bcf85b162..54a64033d 100644 --- a/lumicks/pylake/kymotracker/tests/reference_data/denoising_basic/output_image_2d.npz +++ b/lumicks/pylake/kymotracker/tests/reference_data/denoising_basic/output_image_2d.npz @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:38838017aa036aef256767230b5603494b049ad1eec378931ffed43e834e0a85 +oid sha256:d94c8b51348c9a459fb961c2544a2bf9760e057c320c6fe0a27d87491d138fb2 size 13712 diff --git a/lumicks/pylake/kymotracker/tests/reference_data/denoising_regularized/1d_background.npz b/lumicks/pylake/kymotracker/tests/reference_data/denoising_regularized/1d_background.npz new file mode 100644 index 000000000..2e76b7fe7 --- /dev/null +++ b/lumicks/pylake/kymotracker/tests/reference_data/denoising_regularized/1d_background.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:858ffedfdf118e389e2e6c95f11e388b0a0c1f5edd965224610d112959e549ac +size 13712 diff --git a/lumicks/pylake/kymotracker/tests/reference_data/denoising_regularized/1d_no_background.npz b/lumicks/pylake/kymotracker/tests/reference_data/denoising_regularized/1d_no_background.npz new file mode 100644 index 000000000..9da01ccaf --- /dev/null +++ b/lumicks/pylake/kymotracker/tests/reference_data/denoising_regularized/1d_no_background.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c11d2c5e6001330e6a4146462930d778e8a4694f5745038d340ef7c946eda452 +size 13712 diff --git a/lumicks/pylake/kymotracker/tests/reference_data/denoising_regularized/2d_background.npz b/lumicks/pylake/kymotracker/tests/reference_data/denoising_regularized/2d_background.npz new file mode 100644 index 000000000..78985f825 --- /dev/null +++ b/lumicks/pylake/kymotracker/tests/reference_data/denoising_regularized/2d_background.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e5b479df67e9f1c0aa5092e5dc3dceae7219c15f6ac526a35c202d14256db133 +size 13712 diff --git a/lumicks/pylake/kymotracker/tests/reference_data/denoising_regularized/2d_no_background.npz b/lumicks/pylake/kymotracker/tests/reference_data/denoising_regularized/2d_no_background.npz new file mode 100644 index 000000000..e105611e6 --- /dev/null +++ b/lumicks/pylake/kymotracker/tests/reference_data/denoising_regularized/2d_no_background.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:50cf6835700d853e4778bcd0c29bfab1678b6d48b34b0f392344a7537a1962f0 +size 13712 diff --git a/lumicks/pylake/kymotracker/tests/test_denoise.py b/lumicks/pylake/kymotracker/tests/test_denoise.py index db3f3448f..16c0a288b 100644 --- a/lumicks/pylake/kymotracker/tests/test_denoise.py +++ b/lumicks/pylake/kymotracker/tests/test_denoise.py @@ -161,3 +161,23 @@ def test_denoising_basic(two_dimensional, reference_data): output_image, file_name="output_image_2d" if two_dimensional else "output_image_1d" ), ) + + +@pytest.mark.parametrize( + "two_dimensional, remove_background", + [(False, False), (False, True), (True, False), (True, True)], +) +def test_denoising_regularized(two_dimensional, remove_background, reference_data): + data, truth = denoising_source_image() + kernels = generate_bspline_kernels(3) + msvst = MultiScaleVarianceStabilizingTransform(kernels, two_dimensional=two_dimensional) + image = msvst.filter_regularized( + data, false_detection_rate=0.1, remove_background=remove_background + ) + np.testing.assert_allclose( + image, + reference_data( + image, + file_name=f"{1 + two_dimensional}d_{'no_' if remove_background else ''}background", + ), + )