From caf101235bafa44552278785c0bdbad54fbe1d29 Mon Sep 17 00:00:00 2001 From: ngergihun Date: Thu, 5 Jun 2025 11:38:21 +0200 Subject: [PATCH 1/5] SpikeRemoval not finished --- pySNOM/images.py | 114 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 114 insertions(+) diff --git a/pySNOM/images.py b/pySNOM/images.py index f4a5070..e0602db 100644 --- a/pySNOM/images.py +++ b/pySNOM/images.py @@ -7,6 +7,9 @@ from skimage.transform import warp from skimage.registration import optical_flow_tvl1, phase_cross_correlation from scipy.ndimage import fourier_shift +from scipy.sparse import coo_matrix +from scipy.sparse.linalg import spsolve +from scipy.ndimage import generic_filter MeasurementModes = Enum( "MeasurementModes", @@ -372,6 +375,117 @@ def calculate(self, data, mask=None): def transform(self, data, mask=None): return MaskedTransformation.transform(self, data, mask=mask) +class LaplaceFillIn(Transformation): + """ + Fill in missing (masked) regions of data using inward + interpolation via Laplace's equation. Handles edge and corner cases. + Original NATLAB code: https://github.com/EvanCzako/image-spike-removal/blob/master/remove_spikes.m + """ + + def __init__(self, mask): + self.mask = mask + + def transform(self,data): + """ + Parameters: + - data (2D np.ndarray): Input data. + - mask (2D np.ndarray): Boolean mask where True indicates missing values to fill. + + Returns: + - filled (2D np.ndarray): Image with missing values filled. + """ + + M, N = data.shape + num_pixels = M * N + + # Flattened indices + u = np.flatnonzero(self.mask) # masked (unknown) pixels + w = np.flatnonzero(~self.mask) # known pixels + + # Neighbor index offsets + u_north = u - 1 + u_north = np.where(u % M != 0, u_north, 0) # Wrap prevention for top row + u_east = u + M + u_east = np.where(u_east < num_pixels, u_east, 0) + u_south = u + 1 + u_south = np.where((u + 1) % M != 0, u_south, 0) + u_west = u - M + u_west = np.where(u_west >= 0, u_west, 0) + + a = np.stack([u_north, u_east, u_south, u_west], axis=1) + b = (a > 0).astype(float) + sum_b = b.sum(axis=1, keepdims=True) + c = -b / np.maximum(sum_b, 1e-12) + + # Sparse matrix entries + row_inds = np.concatenate([u, u, u, u, u]) + col_inds = np.concatenate([u, u_north, u_east, u_south, u_west]) + data_vals = np.concatenate([ + np.ones(len(u)), + c[:, 0], c[:, 1], c[:, 2], c[:, 3] + ]) + + # Remove invalid entries + valid = (col_inds >= 0) & (col_inds < num_pixels) + row_inds = row_inds[valid] + col_inds = col_inds[valid] + data_vals = data_vals[valid] + + # Include identity rows for known pixels + row_inds = np.concatenate([row_inds, w]) + col_inds = np.concatenate([col_inds, w]) + data_vals = np.concatenate([data_vals, np.ones(len(w))]) + + # Build sparse matrix + A = coo_matrix((data_vals, (row_inds, col_inds)), shape=(num_pixels, num_pixels)).tocsr() + + # Build RHS vector + b_vec = data.flatten() + b_vec[self.mask.flatten()] = 0 + + # Solve linear system + x = spsolve(A, b_vec) + filled = x.reshape(data.shape) + + return filled + +class RemoveSpikes(Transformation): + def __init__(self, threshold = 0.8): + self.threshold = threshold + + def transfrom(self, data): + """ + Remove impulse-like spikes from an image by replacing values that are + statistical outliers compared to their neighbors. + + Parameters: + - im_in (2D np.ndarray): Input grayscale image. + - thresh (float): Threshold in standard deviations to define spikes. + + Returns: + - z (2D np.ndarray): Image with spikes removed. + """ + z = np.array(data, dtype=float) + M, N = z.shape + + # Define a 3x3 neighborhood excluding the center + def is_spike(val): + center = val[4] + neighbors = np.concatenate([val[:4], val[5:]]) + mean = np.mean(neighbors) + std = np.std(neighbors) + return np.isnan(center) or abs(center - mean) > self.threshold * std + + # Mask of spike pixels (True where spike is detected) + spike_mask = generic_filter(z, is_spike, size=3, mode='mirror') + + # Replace spikes with NaN for region filling + z[spike_mask.astype(bool)] = np.nan + + # Use previously defined fill_region to fill spikes + filled = LaplaceFillIn(z, ).transform(np.isnan(z)) + return filled + class ScarRemoval(Transformation): def __init__(self, threshold=0.5, flip=False, datatype=DataTypes.Phase): From f23f0f735f4c0f8bc1b9be82c386bcbb3c8029b1 Mon Sep 17 00:00:00 2001 From: Gergely Nemeth Date: Thu, 5 Jun 2025 17:24:35 +0200 Subject: [PATCH 2/5] Simple thresholding for spike recognition --- pySNOM/images.py | 44 ++++++++++++-------------------------------- 1 file changed, 12 insertions(+), 32 deletions(-) diff --git a/pySNOM/images.py b/pySNOM/images.py index e0602db..228afba 100644 --- a/pySNOM/images.py +++ b/pySNOM/images.py @@ -450,41 +450,21 @@ def transform(self,data): return filled class RemoveSpikes(Transformation): - def __init__(self, threshold = 0.8): + def __init__(self, threshold = 0.8, fillin = True): self.threshold = threshold + self.fillin = fillin - def transfrom(self, data): - """ - Remove impulse-like spikes from an image by replacing values that are - statistical outliers compared to their neighbors. - - Parameters: - - im_in (2D np.ndarray): Input grayscale image. - - thresh (float): Threshold in standard deviations to define spikes. - - Returns: - - z (2D np.ndarray): Image with spikes removed. - """ - z = np.array(data, dtype=float) - M, N = z.shape - - # Define a 3x3 neighborhood excluding the center - def is_spike(val): - center = val[4] - neighbors = np.concatenate([val[:4], val[5:]]) - mean = np.mean(neighbors) - std = np.std(neighbors) - return np.isnan(center) or abs(center - mean) > self.threshold * std - - # Mask of spike pixels (True where spike is detected) - spike_mask = generic_filter(z, is_spike, size=3, mode='mirror') - - # Replace spikes with NaN for region filling - z[spike_mask.astype(bool)] = np.nan - + def transform(self, data): + norm_data = data/np.nanmedian(data) + spike_mask = mask_from_datacondition(norm_data < self.threshold) + # Use previously defined fill_region to fill spikes - filled = LaplaceFillIn(z, ).transform(np.isnan(z)) - return filled + if self.fillin: + data = LaplaceFillIn(np.isnan(spike_mask)).transform(data) + else: + data = spike_mask*data + + return data class ScarRemoval(Transformation): From 639b03043b727d0bbcce21c60292bfc13e92488e Mon Sep 17 00:00:00 2001 From: Gergely Nemeth Date: Thu, 5 Jun 2025 23:40:24 +0200 Subject: [PATCH 3/5] Test for RemoveSpike --- pySNOM/tests/test_transform.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/pySNOM/tests/test_transform.py b/pySNOM/tests/test_transform.py index 223a0a1..d04eeb3 100644 --- a/pySNOM/tests/test_transform.py +++ b/pySNOM/tests/test_transform.py @@ -9,6 +9,7 @@ SimpleNormalize, DataTypes, AlignImageStack, + RemoveSpikes, ScarRemoval, mask_from_datacondition, dict_from_imagestack, @@ -212,6 +213,16 @@ def test_min(self): out = l.transform(d, mask=mask) np.testing.assert_almost_equal(out, [-1.0, 0.0, 1.0]) +class TestRemoveSpikes(unittest.TestCase): + def test_remove(self): + d = np.ones([9, 9]) + d[4, 4] = 0.8 + + l = RemoveSpikes(threshold=0.9,fillin=True) + + out = l.transform(d) + np.testing.assert_almost_equal(out[4,4],1.0) + class TestAlignImageStack(unittest.TestCase): def test_stackalignment(self): From 13ec06f0abfd1bbc566e12e02fa188eae5bda0b6 Mon Sep 17 00:00:00 2001 From: ngergihun Date: Sat, 7 Jun 2025 15:07:58 +0200 Subject: [PATCH 4/5] Higher option for SpikeRemoval with tests --- pySNOM/images.py | 10 +++++++--- pySNOM/tests/test_transform.py | 9 +++++++++ 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/pySNOM/images.py b/pySNOM/images.py index 228afba..b9d6c79 100644 --- a/pySNOM/images.py +++ b/pySNOM/images.py @@ -450,13 +450,17 @@ def transform(self,data): return filled class RemoveSpikes(Transformation): - def __init__(self, threshold = 0.8, fillin = True): + def __init__(self, threshold = 0.8, higher=False, fillin = True): self.threshold = threshold self.fillin = fillin - + self.higher = higher + def transform(self, data): norm_data = data/np.nanmedian(data) - spike_mask = mask_from_datacondition(norm_data < self.threshold) + if self.higher: + spike_mask = mask_from_datacondition(norm_data > self.threshold) + else: + spike_mask = mask_from_datacondition(norm_data < self.threshold) # Use previously defined fill_region to fill spikes if self.fillin: diff --git a/pySNOM/tests/test_transform.py b/pySNOM/tests/test_transform.py index d04eeb3..0ca4c24 100644 --- a/pySNOM/tests/test_transform.py +++ b/pySNOM/tests/test_transform.py @@ -223,6 +223,15 @@ def test_remove(self): out = l.transform(d) np.testing.assert_almost_equal(out[4,4],1.0) + def test_remove_higher(self): + d = np.ones([9, 9]) + d[4, 4] = 1.2 + + l = RemoveSpikes(threshold=1.1,higher=True,fillin=True) + + out = l.transform(d) + np.testing.assert_almost_equal(out[4,4],1.0) + class TestAlignImageStack(unittest.TestCase): def test_stackalignment(self): From db985ddd247241fa75a32008665cb6eca88f03cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gergely=20N=C3=A9meth?= Date: Mon, 9 Jun 2025 17:56:52 +0200 Subject: [PATCH 5/5] SpikeRemoval updated, new ValueFillIn, new tests --- pySNOM/images.py | 35 ++++++++++++++++++++----- pySNOM/tests/test_transform.py | 47 +++++++++++++++++++++++++++++++--- 2 files changed, 71 insertions(+), 11 deletions(-) diff --git a/pySNOM/images.py b/pySNOM/images.py index b9d6c79..2009036 100644 --- a/pySNOM/images.py +++ b/pySNOM/images.py @@ -448,25 +448,46 @@ def transform(self,data): filled = x.reshape(data.shape) return filled + +class ValueFillIn(Transformation): + def __init__(self, mask, value): + self.value = value + self.mask = mask + + def transform(self, data): + data[np.isnan(self.mask)] = self.value + + return data class RemoveSpikes(Transformation): - def __init__(self, threshold = 0.8, higher=False, fillin = True): + def __init__(self, threshold=0.8, absolute_threshold=False, method='laplace', higher=False, value = 1.0): self.threshold = threshold - self.fillin = fillin + self.method = method self.higher = higher + self.value = value + self.absolute_threshold = absolute_threshold def transform(self, data): - norm_data = data/np.nanmedian(data) + if not self.absolute_threshold: + norm_data = data/np.nanmedian(data) + else: + norm_data = data + if self.higher: spike_mask = mask_from_datacondition(norm_data > self.threshold) else: spike_mask = mask_from_datacondition(norm_data < self.threshold) # Use previously defined fill_region to fill spikes - if self.fillin: - data = LaplaceFillIn(np.isnan(spike_mask)).transform(data) - else: - data = spike_mask*data + match self.method: + case "laplace": + data = LaplaceFillIn(np.isnan(spike_mask)).transform(data) + case "median": + data = ValueFillIn(spike_mask,np.nanmedian(data)).transform(data) + case "manual": + data = ValueFillIn(spike_mask,self.value).transform(data) + case _: + data = spike_mask*data return data diff --git a/pySNOM/tests/test_transform.py b/pySNOM/tests/test_transform.py index 0ca4c24..f532479 100644 --- a/pySNOM/tests/test_transform.py +++ b/pySNOM/tests/test_transform.py @@ -10,6 +10,7 @@ DataTypes, AlignImageStack, RemoveSpikes, + ValueFillIn, ScarRemoval, mask_from_datacondition, dict_from_imagestack, @@ -213,25 +214,63 @@ def test_min(self): out = l.transform(d, mask=mask) np.testing.assert_almost_equal(out, [-1.0, 0.0, 1.0]) +class TestFillIn(unittest.TestCase): + def test_value_fillin(self): + d = np.ones([9, 9]) + mask = np.random.rand(9,9) + mask[4,4] = np.nan + + l = ValueFillIn(mask=mask,value=np.inf) + + out = l.transform(d) + np.testing.assert_equal(out[4,4],np.inf) + class TestRemoveSpikes(unittest.TestCase): - def test_remove(self): + def test_remove_laplace(self): d = np.ones([9, 9]) d[4, 4] = 0.8 - l = RemoveSpikes(threshold=0.9,fillin=True) + l = RemoveSpikes(threshold=0.9,method='laplace') out = l.transform(d) np.testing.assert_almost_equal(out[4,4],1.0) - def test_remove_higher(self): + def test_remove_higher_laplace(self): d = np.ones([9, 9]) d[4, 4] = 1.2 - l = RemoveSpikes(threshold=1.1,higher=True,fillin=True) + l = RemoveSpikes(threshold=1.1,method='laplace',higher=True) out = l.transform(d) np.testing.assert_almost_equal(out[4,4],1.0) + def test_remove_manual(self): + d = np.ones([9, 9]) + d[4, 4] = 0.8 + + l = RemoveSpikes(threshold=0.9,method='manual',value=1.11111) + + out = l.transform(d) + np.testing.assert_almost_equal(out[4,4],1.11111) + + def test_remove_median(self): + d = np.ones([9, 9]) + d[4, 4] = 0.8 + + l = RemoveSpikes(threshold=0.9,method='median') + + out = l.transform(d) + np.testing.assert_almost_equal(out[4,4],1.0) + + def test_remove_nan(self): + d = np.ones([9, 9]) + d[4, 4] = 0.8 + + l = RemoveSpikes(threshold=0.9,method='asdfsdfhdfg') + + out = l.transform(d) + np.testing.assert_almost_equal(out[4,4],np.nan) + class TestAlignImageStack(unittest.TestCase): def test_stackalignment(self):