Skip to content

Add WatershedLungspace class#418

Merged
psomhorst merged 55 commits intodevelopfrom
feature/389-watershed
Aug 25, 2025
Merged

Add WatershedLungspace class#418
psomhorst merged 55 commits intodevelopfrom
feature/389-watershed

Conversation

@psomhorst
Copy link
Copy Markdown
Contributor

This PR adds the WatershedLungspace class.

During development of this class some bugs were fixed and a lot of features were added, mainly to PixelMap and PixelMask. I also updated warnings calls in many files. I recommend going through this PR commit by commit, as I collected relevant changes in atomic commits.

Documentation for the watershed method and a notebook demonstrating its use are included.

pytest changed fixture wrapping in 8.4.0, which changes how fixutures
can be accessed outside of tests.
This makes warnings more readable in general
Update update method to properly update plot_config
A generator was necessary anyway for circular import reasons. While
writing, I remembered I wanted to write a function for non-standard
shapes anyway. This commit kills two birds with one stone.
@psomhorst psomhorst force-pushed the feature/389-watershed branch from e7c867f to a67d7fe Compare August 18, 2025 18:57
def __len__(self):
return self.pixel_impedance.shape[0]

def get_summed_impedance(self, *, return_label: str | None = None, **return_kwargs) -> ContinuousData:
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this for summing impedance within a ROI? Is it necessary to have a separate calculate_global_impedance then? Or I would remove the word "included" in the docstring for calculate_global_impedance since that suggests that there is some sort of selection as well.

Copy link
Copy Markdown
Contributor Author

@psomhorst psomhorst Aug 22, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The idea is that get_summed_impedance will replace the function of calculate_global_impedance. However, for now, calculate_global_impedance is used to calculate the global impedance of the raw signal. The idea is to mark calculate_global_impedance as deprecated in a future version, and remove it altogether in a later version.

Do you think we should already mark it deprecated in this PR? It is not a good idea to make it deprecated in this PR, because it wil result in a DeprecationWarning every time you load EIT data.



@dataclass(frozen=True, init=False)
class IntegerMap(PixelMap):
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe add the IntegerMap and PixelMask plotting examples to the test_pixelmap_imshow notebook as well?

Copy link
Copy Markdown
Contributor Author

@psomhorst psomhorst Aug 22, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added a visual test.

Also, I adapted how an IntegerMap works, since the implementation was a bit hacky. Now, a PixelMap has a dtype argument, making it possible to change how values are casted on init.

Also, there was no unit test for integer maps. I added those.

aggregated_values = aggregator(stacked, axis=0)
return cls(values=aggregated_values, **return_attrs)

def to_boolean_array(self) -> np.ndarray:
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should this method also have a keep_zeros = True option?

Copy link
Copy Markdown
Contributor Author

@psomhorst psomhorst Aug 22, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I made it pm.to_boolean_array(zero=True) and to_non_nan_array(nan=...), both comparable to np.nan_to_num(nan=...).


included_marker_indices = np.isin(peaks_loc_int, markers_inside_tiv_mask)
included_peaks = np.argwhere(included_marker_indices)
excluded_peaks = np.argwhere(~included_marker_indices)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Filtering out the background (marker == 0) might make this capture more meaningful.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are no background pixels in the watershed method. All pixels are associated with one of the peaks. I'm not sure what you are suggesting.


# Find markers that overlap with TIV mask
masked_marker_map = tiv_functional_mask.apply(marker_map)
markers_inside_tiv_mask = masked_marker_map.values[~np.isnan(masked_marker_map.values)]
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If no markers fall inside the TIV mask, should this raise a warning?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think this is mathematically possible, unless, maybe, you set the threshold to a number near (or higher than?) 1. I have no ability to test this.



@pytest.mark.parametrize("threshold", [0.1, 0.15, 0.2, 0.25, 0.3])
def test_watershed_with_simulated_data(simulated_eit_data: Callable[..., EITData], threshold: float):
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If peak_local_max returns an empty list (e.g., flat amplitude), the code may error. Maybe add a test with simulated data where amplitude is flat so no peaks are found, and assert the output mask is empty (all False) or raises a warning.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the amplitude is flat, breath detection would fail before this. I 'fixed' BreathDetection so it returns empty if no peaks are found at all. Then, WatershedLungspace will raise an exception.

assert captures, "captures should have values after runnning apply"
assert isinstance(watershed_mask, PixelMask)
assert np.nansum(watershed_mask.values) > 0
assert len(captures["local peaks"]) # There should be a non-empty list of local peaks
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some captures are checked but not all are asserted. Is there a reason or should the creation and shape of all captures be checked?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm only testing here for things relevant to the real data. I'll test the keys in capture elsewhere.

@JulietteFrancovich
Copy link
Copy Markdown
Contributor

@psomhorst going through the PR per commit was a great tip! I have added my comments, hopefully they are helpful.

@psomhorst
Copy link
Copy Markdown
Contributor Author

@JulietteFrancovich Thanks for the review and for fixing all my language issue. I addressed all your comments.

@psomhorst psomhorst force-pushed the feature/389-watershed branch from 9b8c448 to 6eacd1b Compare August 25, 2025 09:28

@dataclass(kw_only=True, frozen=True)
class TIVLungspace:
"""Create a pixel mask by thresholding the mean TIV or amplitude.
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should the AmplitudeLungspace be a class on its own?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure since the warning is very relevant that amplitude should probably not be used for FLS. However, if someone decides to anyway, it might be misleading if that results in a TIVLungspace object?

"""Create a pixel mask by thresholding the mean amplitude.

This defines the functional lung space as all pixels with an amplitude of at least the provided fractional threshold
of the maximum ampltiude.
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
of the maximum ampltiude.
of the maximum amplitude.

@psomhorst psomhorst force-pushed the feature/389-watershed branch from 61aef48 to f45fc39 Compare August 25, 2025 14:28
@psomhorst psomhorst merged commit 1637bc5 into develop Aug 25, 2025
3 checks passed
@psomhorst psomhorst deleted the feature/389-watershed branch August 25, 2025 14:40
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants