Skip to content
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

Incorporating height profile into GrainStats #871

Merged
merged 4 commits into from
Jul 17, 2024

Conversation

ns-rse
Copy link
Collaborator

@ns-rse ns-rse commented Jul 4, 2024

Closes #748

After the protracted wrestling I had of how to get minimum feret distances out (see #755) I've decided to finish this work off now as I'm unlikely to look at it again until September and had been tinkering with it over the last week and it is more efficient to finish it off now than struggle to remember where I was in two months time.

This PR includes two commits, importantly one fixes a problem with the topostats.measure.height_profiles.interpolate_height_profiles() which was introduced in #859 (details on that below).

Height Profiles

The option extract_height_profile has been introduced to the grainstats section of default_config.yaml, when true this will trigger extraction of height profiles along the maximum feret. These are aggregated across all "grains" and images and written to a JSON file (height_profiles.json) adjacent to all_statistics.csv in the specified output directory.

JSON is used because CSV is not particularly conducive to this data structure as profiles are taken across cropped images which differ in size.

This necessitated adding a dependency for numpyencoder which correctly and safely converts Numpy dtypes to a format that json.dump() can write.

I have a suspicion that at some point it will be desirable to generate line plots of these which is possible but for now users will have to extract the data align and plot these themselves (if anyone does this and wants to share or contribute their solution to TopoStats that would be very much appreciated).

Fixing interpolate_height_profiles() #859

In #859 I had mistakenly used np.sort() on the feret co-ordinates. Two of tests developed there passed, one I opted to skip (The maximum feret was almost horizontal and this should have raised flags as it was indicative of another problem, namely that we should interpolate in such a manner as to get as many points as possible along the line of the maximum feret)

When it came to developing tests to incorporate height profile extraction into the work flow of GrainStats() it revealed that things were not as they should be as whilst the feret coordinates were correct point interpolation resulted in...

grain_with_sort01
grain_with_sort02
grain_with_sort03

Clearly not what we wanted! This arose because np.interp() expects "monotonically increasing sample points" and in some of the above the feret coordinates were returned with the highest value first and were therefore decreasing.

Solution

In order to get the interpolation correct after identifying these problems two changes were made.

  1. Assess which is the largest distance between the feret co-ordinates on the x and y axis and the largest has points linearly interpolated first.
  2. Order the feret points in ascending order (depending on whether it is the x or y axis that is the largest distance as this is the one along which interpolation will occur).

The points are now correctly interpolated regardless of the order in which the feret coordinates are returned and we get the greatest detail by doing so along the longest x/y axis first.

grain_correct01
grain_correct02
grain_correct03

As a bonus the test that was previously skipped is now looking as expected and has been reinstated.

@ns-rse ns-rse force-pushed the ns-rse/748-height-profiles-in-grainstats branch 4 times, most recently from 40156ae to aa8d5e6 Compare July 8, 2024 13:15
Saves height profiles across the maximum feret of cropped grains to JSON (`height_profiles.json`).

Controlled by the `extract_height_profile` option under `grainstats` in `default_config.yaml`.

Regression tests are updated to include this new output and new tests introduced where required.
One of the tests was deliberately skipped and looked weird as only three points were interpolated. On starting to
incorporate into the main workflow tests revealed why this is. The feret coordinates were being sorted independently
which was a mistake, it decoupled and often flipped over the points between which interpolation was occurring.

Further we were always taking the x-axis and interpolating along that but when the y-axis was the longer distance this
resulted in sparse interpolation. Instead we now perform linear interpolation on the axis with longest difference first
and then interpolate the points on the other axis to determine the points at which we want to interpolate heights.
@ns-rse ns-rse force-pushed the ns-rse/748-height-profiles-in-grainstats branch from aa8d5e6 to 3c11103 Compare July 8, 2024 13:29
@ns-rse
Copy link
Collaborator Author

ns-rse commented Jul 8, 2024

Ready for review as I figured out how to update the HDF5 (see #873 for details and future reference).

In fixing this I found there was a precision error in the heights on OSX compared to Linux (where the tests were developed) to which end I split the regression test of heights out to their own functions and skip them if the reported platform.platform()[:5] == "macOS". This adds a slight increased overhead as in essence two of the tests are run twice but assessed differently but I couldn't see an option for format float's when using json.dump().

I also noticed a missing f-string in one of the plottingfuncs.py LOGGER.info() so updated that.

@ns-rse ns-rse marked this pull request as ready for review July 8, 2024 13:59
Copy link
Collaborator

@SylviaWhittle SylviaWhittle left a comment

Choose a reason for hiding this comment

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

This works for me, and I am able to extract the height profiles from the HDF5 file and plot them.

Thank you very much for your work on this 🎉

tests/measure/test_height_profiles.py Outdated Show resolved Hide resolved
topostats/plottingfuncs.py Show resolved Hide resolved
topostats/processing.py Show resolved Hide resolved
@SylviaWhittle
Copy link
Collaborator

Here is how I plotted the heights in case anyone reads this
image

@ns-rse
Copy link
Collaborator Author

ns-rse commented Jul 17, 2024

Here is how I plotted the heights in case anyone reads this image

It should be possible to iterate over the dictionary directly using .items(). Following is untested but should achieve the same thing.

from pathlib import Path
from loguru import logger
import numpy as np
import h5py
import matplotlib.pyplot as plt

FILE_PATH = Path("./minicircle.topostats")
with h5py.File(FILE_PATH, "r") as f:
    # print keys
    logger.info(f.keys())
    height_profiles = f["height_profiles"]["above"]

fig, ax = plt.subplots(figsize=(10,5))
for key, profile in height_profiles.items()
    ax.plot(profile, label=key)
ax.set_xlabel("Pixel")
ax.set_ylabel("Height (nm)")
ax.set_title("Height Profiles")
plt.legend()
plt.show()

Copy link
Collaborator

@SylviaWhittle SylviaWhittle left a comment

Choose a reason for hiding this comment

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

🎉

@ns-rse ns-rse added this pull request to the merge queue Jul 17, 2024
Merged via the queue into main with commit 3cee358 Jul 17, 2024
14 checks passed
@ns-rse ns-rse deleted the ns-rse/748-height-profiles-in-grainstats branch July 17, 2024 13:46
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.

Add horizontal grain profiles
2 participants