-
Notifications
You must be signed in to change notification settings - Fork 10
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add horizontal grain profiles #748
Comments
Ultimately the profile of any horizontal line is available as the Numpy array for the given grain so its "just" a case of finding the centroid of a grain (and skimage has functions for doing this) and extracting that row and the plotting. Some random questions off the top of my head (not necessarily directed at you @SylviaWhittle ) as whilst its relatively straight-forward to implement a solution understanding the context and purpose of a feature can heavily influence the solution that is implemented (e.g. if other profiles based on the widest or thinnest part of a molecule are required then slices need taking from orientations other than horizontal). At what point in the processing should we extract this information, is it from the untraced or traced grain? I suspect the former but would like to check What should be done for linear molecules? In the future is any other profile likely to be required other than that which passes through the centroid? If so will these always be horizontal or are other angles/paths required such as the profile through the widest part of a circular molecule1? Footnotes
|
This would be done before tracing, as an optional add on and would be used on "blob" grains (which are not usually traced) to give a profile spanning the grain and the background from outside to give a trace of the blob height. Useful add ons would be to make this average over a few lines either side or not, depending on the users preference. Horizontally is used to avoid masking errors. |
Currently the classification of grains into A simple solution would be to just plot the required profile regardless of type, even if its meaningless, and that can be done, I was just curious at which stage the cross-section should be taken during the processing pipeline, for example after applying a Gaussian blur the values may be changed. What is the intended use of the data? I ask as I'm curious about why only the horizontal axis is required. These are the same grain with the same heights, but orientated differently.
...and you'd get very different profiles from the horizontal only.
|
I think this would be an option - like “crop” in the config file, Where if the suer wants to they would plot a line of “user defined length” or through the centre of the blob from the centre only. This is to use as single examples alongside the numbers obtained from topostats e.g. height, are volume, partly for sense checking, partly for visualisation.
If you did this over all the molecules, you would get an average. Which would show you heterogeneity you display below.
… On 5 Dec 2023, at 09:52, Neil Shephard ***@***.***> wrote:
Currently the classification of grains into circle or linear is done during tracing and there is no classification into blob.
A simple solution would be to just plot the required profile regardless of type, even if its meaningless, and that can be done, I was just curious at which stage the cross-section should be taken during the processing pipeline, for example after applying a Gaussian blur the values may be changed.
What is the intended use of the data? I ask as I'm curious about why only the horizontal axis is required.
These are the same grain with the same heights, but orientated differently.
0 0 0 0 0
0 1 1 1 0
0 1 2 1 0
0 1 2 1 0
0 1 2 1 0
0 1 2 1 0
0 1 2 1 0
0 1 2 1 0
0 1 1 1 0
0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 1 1 1 1 1 1 1 1 0
0 1 2 2 2 2 2 2 1 0
0 1 1 1 1 1 1 1 1 0
0 0 0 0 0 0 0 0 0 0
...and you'd get very different profiles from the horizontal only.
What information is sought/trying to be extracted?
Is there some sort of step prior to or during scanning that ensures grains are orientated in a specific way?
minicircle.spm looks fairly random to me but I'm very wary of generalising based on a single sample.
—
Reply to this email directly, view it on GitHub <#748 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADJSFMJLTT42LOEU6K6RUYTYH3VG3AVCNFSM6AAAAABAG4DU2WVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNBQGQYDQMRSG4>.
You are receiving this because you commented.
|
Wouldn't you have heterogeneous measurements across a single molecule depending on the orientation though? Would it perhaps be useful to use a slice along the line parallel to the minimum feret diameter which passes through the centroid as that would then pick up the height profile along (roughly1) the longest slice? I'm probably over thinking this though so please do say but I'm trying to anticipate other scenarios and think about general solutions at the outset. Footnotes
|
You could do that, open to measuring the line along the max feret length
through the centroid, or through the max and min feret as a sort of cross.
…On Tue, 5 Dec 2023 at 10:04, Neil Shephard ***@***.***> wrote:
Wouldn't you have heterogeneous measurements across a single molecule
depending on the orientation though?
Would it perhaps be useful to use a slice along the line parallel to the
minimum feret diameter which passes through the centroid as that would then
pick up the height profile along (roughly1
<#m_3666661370039346476_user-content-fn-1-5f08325f9d97e3fadf3048e740bccab3>)
the longest slice?
I'm probably over thinking this though so please do say but I'm trying to
anticipate other scenarios and think about general solutions at the outset.
Footnotes
1.
I say roughly because there is no guarantee the extremes will always
be through the centroid ↩
<#m_3666661370039346476_user-content-fnref-1-5f08325f9d97e3fadf3048e740bccab3>
—
Reply to this email directly, view it on GitHub
<#748 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/ADJSFMJGVFMMDAJ7ZP5GERDYH3WS5AVCNFSM6AAAAABAG4DU2WVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNBQGQZDOMZYGI>
.
You are receiving this because you commented.Message ID:
***@***.***>
|
Further notes after todays meeting with @SylviaWhittle , @llwiggins and @MaxGamill-Sheffield
This feature request should take priority over writing tests for the maxgamill-sheffield/cats branch which underpins the forthcoming paper which is due to be finalised 2024-01. |
I've started work on this and have found something that is potentially useful, its written up in a separate issue, see #750. |
Writing my thoughts down here before coding... Algorithm
How to do thisScikit-image has methods/code for calculating feret distance (see from typing import Tuple
import numpy as np
import numpy.typing as npt
def calculate_feret(image: npt.NDArray) -> Tuple[npt.NDArray, npt.NDArray]:
# Import necessary functions from skimage
from skimage.morphology.convex_hull import convex_hull_image
from skimage.measure import find_contours, marching_cubes
image_convex = convex_hull_image(image)
spacing =
identity_convex_hull = np.pad(image_convex, 2, mode="constant", constant_values=0)
if self._ndim == 2:
coordinates = np.vstack(find_contours(identity_convex_hull, 0.5, fully_connected="high"))
elif self._ndim == 3:
coordinates, _, _, _ = marching_cubes(identity_convex_hull, level=0.5)
distances = pdist(coordinates * self._spacing, "sqeuclidean")
return distances, coordinates We could then easily extract the co-ordinates for the min_feret_coordinates = coordinates[np.argwhere(np.min(distances) == distances],]
max_feret_coordinates = coordinates[np.argwhere(np.min(distances) == distances],] These can be applied to the Questions (mostly to myself)
|
Not so straight-forward as the following shows, the from skimage import draw
from skimage.morphology.convex_hull import convex_hull_image
from skimage.measure import find_contours, marching_cubes
from skimage.measure._regionprops_utils import _normalize_spacing
from scipy.spatial.distance import pdist, squareform
holo_circle = np.zeros((13, 13), dtype=np.uint8)
rr, cc = draw.circle_perimeter(6, 6, 5)
holo_circle[rr, cc] = 1
holo_circle
holo_ellipse_vertical = np.zeros((16, 10), dtype=np.uint8)
rr, cc = draw.ellipse_perimeter(8, 5, 3, 6, orientation=np.deg2rad(90))
holo_ellipse_vertical[rr, cc] = 1
holo_ellipse_vertical
holo_ellipse_horizontal = np.zeros((10, 16), dtype=np.uint8)
rr, cc = draw.ellipse_perimeter(5, 8, 6, 3, orientation=np.deg2rad(90))
holo_ellipse_horizontal[rr, cc] = 1
holo_ellipse_horizontal
holo_ellipse_angled = np.zeros((12, 14), dtype=np.uint8)
rr, cc = draw.ellipse_perimeter(6, 7, 3, 5, orientation=np.deg2rad(30))
holo_ellipse_angled[rr, cc] = 1
holo_ellipse_angled
curved_line = np.zeros((10, 10), dtype=np.uint8)
rr, cc = draw.bezier_curve(1, 5, 5, -2, 8, 8, 2)
curved_line[rr, cc] = 1
curved_line
def calculate_feret(image: npt.NDArray, spacing = None) -> Tuple[int, int]:
"""Calculate the feret distances.
NB - This method is taken from that implemented in 'skimage.measure.regionprops()' which calculates the maximum
feret distance.
Parameters
----------
image : npt.NDArray
Labelled image.
Returns
-------
Tuple
Tuple of the distances and the corresponding co-ordinates.
"""
image_convex = convex_hull_image(image)
if spacing is None:
spacing = np.full(image.ndim, 1.0)
_spacing = _normalize_spacing(spacing, image.ndim)
identity_convex_hull = np.pad(image_convex, 2, mode="constant", constant_values=0)
if image.ndim == 2:
coordinates = np.vstack(find_contours(identity_convex_hull, 0.5, fully_connected="high"))
elif image.ndim == 3:
coordinates, _, _, _ = marching_cubes(identity_convex_hull, level=0.5)
print(f"{coordinates}")
distances = pdist(coordinates * spacing, "sqeuclidean")
return distances, coordinates
holo_circle_distances, holo_circle_coordinates = calculate_feret(holo_circle) Fortunately I found in this issue reference to calculating the minimum feret distance for 2-D objects written up in a gist by @VolkerH so will give that a whirl. |
New submodule `measures.feret` with functions for calculation of min and max feret and the associated co-ordinates. Utilises code from Skan developer see [gist](https://gist.github.com/VolkerH/0d07d05d5cb189b56362e8ee41882abf) and as [suggested](scikit-image/scikit-image#4817 (comment)) it adds tests. In doing so I found sorting points prior to calculation of upper and lower convex hulls missed some points. I'm also not sure about the `rotating_calipers()` function as it doesn't seem to return all pairs formed across the points from the upper and lower convex hulls and so have included but not used the `all_pairs()` function which does, although if used in place this doesn't return the minimum feret correctly, rather just the smallest distance. Currently some of the tests for the `curved_line` fail too. The intention is to use this without instantiating the `GrainStats` class to be used in profiles (#748) and could serve as a concise stand-alone set of functions outside of `GrainStats` class which currently has static methods for the calculations (#750). Benchmarking In light of #750 and re-implementing feret calculations as a new sub-module I wanted to see how this compares in terms of performance to those in `GrainStats()` class so have run some very basic benchmarking. Note this is not the full pipeline of taking labelled images and finding the outlines/edge points which are required as inputs to the calculations, its purely on the calculation of min-max feret from the edge points. ``` import timeit import numpy as np from skimage import draw from topostats.measure import feret from topostats.grainstats import GrainStats holo_circle = np.zeros((14, 14), dtype=np.uint8) rr, cc = draw.circle_perimeter(6, 6, 5) holo_circle[rr, cc] = 1 holo_circle_edge_points = np.argwhere(holo_circle == 1) %timeit feret.min_max_feret(holo_circle_edge_points) 83 µs ± 686 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each) %timeit GrainStats.get_max_min_ferets(holo_circle_edge_points) 1.06 ms ± 10.4 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each) ``` So this new implementation is faster.
New submodule `measures.feret` with functions for calculation of min and max feret and the associated co-ordinates. Utilises code from Skan developer see [gist](https://gist.github.com/VolkerH/0d07d05d5cb189b56362e8ee41882abf) and as [suggested](scikit-image/scikit-image#4817 (comment)) it adds tests. In doing so I found sorting points prior to calculation of upper and lower convex hulls missed some points. I'm also not sure about the `rotating_calipers()` function as it doesn't seem to return all pairs formed across the points from the upper and lower convex hulls and so have included but not used the `all_pairs()` function which does, although if used in place this doesn't return the minimum feret correctly, rather just the smallest distance. Currently some of the tests for the `curved_line` fail too. The intention is to use this without instantiating the `GrainStats` class to be used in profiles (#748) and could serve as a concise stand-alone set of functions outside of `GrainStats` class which currently has static methods for the calculations (#750). Benchmarking In light of #750 and re-implementing feret calculations as a new sub-module I wanted to see how this compares in terms of performance to those in `GrainStats()` class so have run some very basic benchmarking. Note this is not the full pipeline of taking labelled images and finding the outlines/edge points which are required as inputs to the calculations, its purely on the calculation of min-max feret from the edge points. ``` import timeit import numpy as np from skimage import draw from topostats.measure import feret from topostats.grainstats import GrainStats holo_circle = np.zeros((14, 14), dtype=np.uint8) rr, cc = draw.circle_perimeter(6, 6, 5) holo_circle[rr, cc] = 1 holo_circle_edge_points = np.argwhere(holo_circle == 1) %timeit feret.min_max_feret(holo_circle_edge_points) 83 µs ± 686 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each) %timeit GrainStats.get_max_min_ferets(holo_circle_edge_points) 1.06 ms ± 10.4 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each) ``` So this new implementation is faster.
See comments in #755 (in particular here). The algorithm proposed in the Mathworks article works for everything except triangles and it took me a while (most of yesterday) to work through why, carefully comparing the method @SylviaWhittle implemented in
The clause is because of the condition... next_index = 1 if index == n_pairs - 1 else index + 1 ...which on the final iteration means the upper hull is using the last element ([1, 1]) and the first element ([2, 1]) which are the same points on the lower hull which have already been seen as the base. It also means that the base formed by [1, 2] and [2, 1] has not been observed which is why I get the minimum height returned as I can't work out a way of adapting this algorithm so have gone with @SylviaWhittle solution of calculating the minimum feret based on the triangle formed by a rotating caliper pair (of which one point is the apex) and the next adjacent point on one of the hulls (albeit with some tweaks to the code layout). Minimum Feret Co-OrdinatesAfter some revision the stated desire is to have...
This is probably most easily obtained by determining the mid-point between the maximum feret co-ordinates and taking a profile that is prependicular to the maximum feret at this point. But I can envisage two issues based on the limited sample of images I've seen.
For 1) I think perhaps returning both the suggested smaller profiles might be helpful, and a clear description of what these represents should be included in the documentations data dictionary that describes the grain statistics that are returned. For 2) I propose taking the nearest co-ordinate to the mid-point of the base of the triangle formed in calculating the minimum feret distance and returning the profile between this and the apex of the triangle. A very crude schematic to maybe aid with visualisation (would take too long for me to plot this out).
|
Further to this adds the `topostats.measure.height_profiles` sub-module which for a given pair of feret coordinates... 1. Determines the orientation of of the feret co-ordinates relative to horizontal. 2. Rotates the image so that the feret is horizontal. 3. Recalculates the co-ordinates of the feret after rotation. 4. Extracts the height profile along the feret. Includes a function in `topostats.plotting.plot_height_profiles()` which produces a line-plot of multiple height profiles. ToDo... Still a fair few steps to integrate this into the processing. + Add configuration options to `topostats/default_config.yaml` of whether to calculate `height_profiles`. + Update `GrainStats()` to calculate the height profile for the image being processed if required. + Return the `height_profile` (1-D Numpy array). + Collect `height_profile` acrss grains into a dictionary (may require current code as written to be adapted to work with dictionaries, currently works with lists in `plot_height_profiles()`). + Add functionality to write the profiles to JSON for subsequent use/plotting (e.g. customising style/axis labels/etc. of plot) Related : #748 #755
Further to this adds the `topostats.measure.height_profiles` sub-module which for a given pair of feret coordinates... 1. Determines the orientation of of the feret co-ordinates relative to horizontal. 2. Rotates the image so that the feret is horizontal. 3. Recalculates the co-ordinates of the feret after rotation. 4. Extracts the height profile along the feret. Includes a function in `topostats.plotting.plot_height_profiles()` which produces a line-plot of multiple height profiles. ToDo... Still a fair few steps to integrate this into the processing. + Add configuration options to `topostats/default_config.yaml` of whether to calculate `height_profiles`. + Update `GrainStats()` to calculate the height profile for the image being processed if required. + Return the `height_profile` (1-D Numpy array). + Collect `height_profile` acrss grains into a dictionary (may require current code as written to be adapted to work with dictionaries, currently works with lists in `plot_height_profiles()`). + Add functionality to write the profiles to JSON for subsequent use/plotting (e.g. customising style/axis labels/etc. of plot) Related : #748 #755
Further to #755 this adds the `topostats.measure.height_profiles` sub-module which interpolates the heights between the maximum feret co-ordinates. Includes a function in `topostats.plotting.plot_height_profiles()` which produces a line-plot of multiple height profiles. ToDo... Still a fair few steps to integrate this into the processing. + Add configuration options to `topostats/default_config.yaml` of whether to calculate `height_profiles`. + Add configuration option for the `scipy.interpolate.RegularGridInterpolator()` options which are passed via `**kwargs`. + Update `GrainStats()` to calculate the height profile for the image being processed if required. + Return the `height_profile` (1-D Numpy array). + Collect `height_profile` acrss grains into a dictionary (may require current code as written to be adapted to work with dictionaries, currently works with lists in `plot_height_profiles()`). + Add functionality to write the profiles to JSON for subsequent use/plotting (e.g. customising style/axis labels/etc. of plot) Related : #748 #755
Further to #755 this adds the `topostats.measure.height_profiles` sub-module which interpolates the heights between the maximum feret co-ordinates. Includes a function in `topostats.plotting.plot_height_profiles()` which produces a line-plot of multiple height profiles. ToDo... Still a fair few steps to integrate this into the processing. + Add configuration options to `topostats/default_config.yaml` of whether to calculate `height_profiles`. + Add configuration option for the `scipy.interpolate.RegularGridInterpolator()` options which are passed via `**kwargs`. + Update `GrainStats()` to calculate the height profile for the image being processed if required. + Return the `height_profile` (1-D Numpy array). + Collect `height_profile` acrss grains into a dictionary (may require current code as written to be adapted to work with dictionaries, currently works with lists in `plot_height_profiles()`). + Add functionality to write the profiles to JSON for subsequent use/plotting (e.g. customising style/axis labels/etc. of plot) Related : #748 #755
Further to #755 this adds the `topostats.measure.height_profiles` sub-module which interpolates the heights between the maximum feret co-ordinates. Includes a function in `topostats.plotting.plot_height_profiles()` which produces a line-plot of multiple height profiles. ToDo... Still a fair few steps to integrate this into the processing. + Add configuration options to `topostats/default_config.yaml` of whether to calculate `height_profiles`. + Add configuration option for the `scipy.interpolate.RegularGridInterpolator()` options which are passed via `**kwargs`. + Update `GrainStats()` to calculate the height profile for the image being processed if required. + Return the `height_profile` (1-D Numpy array). + Collect `height_profile` acrss grains into a dictionary (may require current code as written to be adapted to work with dictionaries, currently works with lists in `plot_height_profiles()`). + Add functionality to write the profiles to JSON for subsequent use/plotting (e.g. customising style/axis labels/etc. of plot) Related : #748 #755
Further to #755 this adds the `topostats.measure.height_profiles` sub-module which interpolates the heights between the maximum feret co-ordinates. Includes a function in `topostats.plotting.plot_height_profiles()` which produces a line-plot of multiple height profiles. ToDo... Still a fair few steps to integrate this into the processing. + Add configuration options to `topostats/default_config.yaml` of whether to calculate `height_profiles`. + Add configuration option for the `scipy.interpolate.RegularGridInterpolator()` options which are passed via `**kwargs`. + Update `GrainStats()` to calculate the height profile for the image being processed if required. + Return the `height_profile` (1-D Numpy array). + Collect `height_profile` acrss grains into a dictionary (may require current code as written to be adapted to work with dictionaries, currently works with lists in `plot_height_profiles()`). + Add functionality to write the profiles to JSON for subsequent use/plotting (e.g. customising style/axis labels/etc. of plot) Related : #748 #755
Further to #755 this adds the `topostats.measure.height_profiles` sub-module which interpolates the heights between the maximum feret co-ordinates. Includes a function in `topostats.plotting.plot_height_profiles()` which produces a line-plot of multiple height profiles. ToDo... Still a fair few steps to integrate this into the processing. + Add configuration options to `topostats/default_config.yaml` of whether to calculate `height_profiles`. + Add configuration option for the `scipy.interpolate.RegularGridInterpolator()` options which are passed via `**kwargs`. + Update `GrainStats()` to calculate the height profile for the image being processed if required. + Return the `height_profile` (1-D Numpy array). + Collect `height_profile` acrss grains into a dictionary (may require current code as written to be adapted to work with dictionaries, currently works with lists in `plot_height_profiles()`). + Add functionality to write the profiles to JSON for subsequent use/plotting (e.g. customising style/axis labels/etc. of plot) Related : #748 #755
Further to #755 this adds the `topostats.measure.height_profiles` sub-module which interpolates the heights between the maximum feret co-ordinates. Includes a function in `topostats.plotting.plot_height_profiles()` which produces a line-plot of multiple height profiles. ToDo... Still a fair few steps to integrate this into the processing. + Add configuration options to `topostats/default_config.yaml` of whether to calculate `height_profiles`. + Add configuration option for the `scipy.interpolate.RegularGridInterpolator()` options which are passed via `**kwargs`. + Update `GrainStats()` to calculate the height profile for the image being processed if required. + Return the `height_profile` (1-D Numpy array). + Collect `height_profile` acrss grains into a dictionary (may require current code as written to be adapted to work with dictionaries, currently works with lists in `plot_height_profiles()`). + Add functionality to write the profiles to JSON for subsequent use/plotting (e.g. customising style/axis labels/etc. of plot) Related : #748 #755
Is your feature request related to a problem? Please describe.
It has been requested that topostats have a feature similar to Gwyddion where the horizontal (in the image) height profile of a grain would be taken. In other words, a horizontal line is placed over the molecule, and the pixels that form the line have their heights measured to form a profile for that molecule.
Describe the solution you'd like
During grainstats, a horizontal, one pixel thick line of pixels should be formed, with the length sufficient to cover the entire molecule from left to right, and the line intersecting the centroid of the molecule. The heights for those pixels should be obtained and saved as a line plot. Probably only saved with image_set all.
While it does seem arbitrary to just take one line through the centre, and we will miss a lot of information, this is is not an issue, as this line profile is just to get an idea of the height profile of the molecule.
This is a heavily requested feature and should be quite high priority.
It should be relatively trivial to implement. (I think?)
The text was updated successfully, but these errors were encountered: