now the globules area estimator's results are stored as hdf5 files

- saving results as hdf5 allows easy usage or exploration of results using hdfview or hdfcompass, or by any that supports hdf5 files.
- I had to add support for hdf5 file export, which was more complicate than using h5py library (which is not available in fiji's  jython). So instead of using h5py, I wrote a generic hdf5 data structure (not restricted to jython), that can be saved using the jython jhdf5 library.
This commit is contained in:
Guillaume Raffy 2020-03-20 14:36:21 +01:00
parent f4aa0423ac
commit 22c112a18d
5 changed files with 263 additions and 24 deletions

View File

@ -0,0 +1 @@
from hdf5_data import Group, DataSet, ElementType

View File

@ -0,0 +1,103 @@
"""structures and classes for storing hdf5 data
see https://portal.hdfgroup.org/display/HDF5/Introduction+to+HDF5
note: these structures are intentionally made independent of h5py library as this h5py library is not available in all python environments. More specifically, h5py is not available in jython. So the purpose of this library is to allow storage of hdf5 data in all environments.
"""
class ElementType(object):
U1=0
U8=1
U32=2
S32=3
F32=4
F64=5
class DataSet(object):
""" an hdf5 dataset
"""
def __init__(self, name, size, element_type=ElementType.U32):
"""
:param str name: the name of the dataset
:param tuple(int) size: the size of each dimension of the dataset
"""
self.name = name
self.size = size
self._element_type = element_type
self._elements = []
for i in range(self.num_elements): # pylint: disable=unused-variable
self._elements.append(0)
@property
def element_type(self):
return self._element_type
@property
def dimension(self):
return len(self.size)
@property
def elements(self):
return self._elements
@property
def num_elements(self):
num_elements = 1
for i in range(self.dimension):
num_elements *= self.size[i]
return num_elements
def get_element(self, pos):
assert len(pos) == len(self.size)
return self._elements[self._pos_to_index(pos)]
def __getitem__(self, pos):
assert len(pos) == len(self.size)
return self._elements[self._pos_to_index(pos)]
def set_element(self, pos, value):
assert len(pos) == len(self.size)
self._elements[self._pos_to_index(pos)] = value
def __setitem__(self, pos, value):
assert len(pos) == len(self.size)
self._elements[self._pos_to_index(pos)] = value
def _pos_to_index(self, pos):
index = 0
for i in range(self.dimension):
index += pos[i]
if i < self.dimension - 1:
index *= self.size[i]
return index
class Group(object):
def __init__(self, name):
"""
:param str name: the name of the group
"""
self.name = name
self._datasets = {}
self._groups = {}
def add_dataset(self, dataset):
assert isinstance(dataset, DataSet)
self._datasets[dataset.name] = dataset
def add_group(self, group):
assert isinstance(group, Group)
self._datasets[group.name] = group
@property
def datasets(self):
return self._datasets
@property
def groups(self):
return self._groups
def __getitem__(self, dataset_name):
return self._datasets[dataset_name]

View File

@ -0,0 +1,137 @@
# h5py is not available in fiji's jython
from ncsa.hdf.hdf5lib import HDFArray # found in FIJI_HOME/jars/jhdf5-14.12.6.jar pylint: disable=import-error
# samples in https://support.hdfgroup.org/HDF5/examples/hdf-java.html
from ncsa.hdf.hdf5lib.H5 import H5Fcreate, H5Fclose # pylint: disable=import-error
from ncsa.hdf.hdf5lib.H5 import H5Gcreate, H5Gclose # pylint: disable=import-error
from ncsa.hdf.hdf5lib.H5 import H5Screate_simple, H5Sclose # pylint: disable=import-error
from ncsa.hdf.hdf5lib.H5 import H5Dcreate, H5Dclose, H5Dwrite, H5Dget_type # pylint: disable=import-error
from ncsa.hdf.hdf5lib.HDF5Constants import H5F_ACC_TRUNC # pylint: disable=import-error
from ncsa.hdf.hdf5lib.HDF5Constants import H5S_ALL # pylint: disable=import-error
from ncsa.hdf.hdf5lib.HDF5Constants import H5P_DEFAULT # pylint: disable=import-error
from ncsa.hdf.hdf5lib.HDF5Constants import H5T_STD_U8LE, H5T_STD_U8LE, H5T_STD_U32LE, H5T_STD_I32LE, H5T_IEEE_F32LE, H5T_IEEE_F64LE, H5T_NATIVE_INT32 # pylint: disable=import-error
from jarray import zeros, array # pylint: disable=import-error
from lipase.hdf5 import Group, ElementType, DataSet
def create_hdf5_test_file1(filepath):
hdf5_file = H5Fcreate(filepath, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)
group_id1 = H5Gcreate(hdf5_file, 'g1', H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)
dataspace_id1 = H5Screate_simple(2, [20, 10], None)
dataset_id = H5Dcreate(group_id1, "2D 32-bit integer 20x10", H5T_NATIVE_INT32, dataspace_id1, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)
H5Dclose(dataset_id)
H5Sclose(dataspace_id1)
H5Gclose(group_id1)
H5Fclose(hdf5_file)
def create_hdf5_test_file2(filepath):
hdf5_file = H5Fcreate(filepath, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)
group_id1 = H5Gcreate(hdf5_file, 'g1', H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)
dataspace_id1 = H5Screate_simple(1, [5], None)
dataset_id = H5Dcreate(group_id1, "1D 32-bit integer 5", H5T_NATIVE_INT32, dataspace_id1, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)
refs=array([1,2,3,4,5], 'i')
tid = H5Dget_type(dataset_id)
# sid = H5.H5Dget_space(did);
H5Dwrite(dataset_id, tid, H5S_ALL, H5S_ALL, H5P_DEFAULT, refs)
# int[][] dataRead = new int[(int) dims2D[0]][(int) (dims2D[1])];
# try {
# if (dataset_id >= 0)
# H5.H5Dread(dataset_id, HDF5Constants.H5T_NATIVE_INT,
# HDF5Constants.H5S_ALL, HDF5Constants.H5S_ALL,
# HDF5Constants.H5P_DEFAULT, dataRead);
H5Dclose(dataset_id)
H5Sclose(dataspace_id1)
H5Gclose(group_id1)
H5Fclose(hdf5_file)
# https://www.jython.org/jython-old-sites/archive/21/docs/jarray.html
ELEMENT_TYPE_TO_JARRAY_TYPE={
ElementType.U1: 'z',
ElementType.U8: 'b',
ElementType.U32: 'i',
ElementType.S32: 'i',
ElementType.F32: 'f',
ElementType.F64: 'd',
}
# https://support.hdfgroup.org/HDF5/doc/UG/HDF5_Users_Guide-Responsive%20HTML5/index.html#t=HDF5_Users_Guide%2FDatatypes%2FHDF5_Datatypes.htm
# https://support.hdfgroup.org/HDF5/doc/RM/PredefDTypes.html
ELEMENT_TYPE_TO_HDF5_TYPE={
ElementType.U1: H5T_STD_U8LE,
ElementType.U8: H5T_STD_U8LE,
ElementType.U32: H5T_STD_U32LE,
ElementType.S32: H5T_STD_I32LE,
ElementType.F32: H5T_IEEE_F32LE,
ElementType.F64: H5T_IEEE_F64LE,
}
def _save_hdf5_data_set(hdf5_container, data_set):
"""
:param hid_t hdf5_file: hdf5 file handle open in write mode or hdf5 group handle
:param DataSet data_set:
"""
dataspace_id = H5Screate_simple(data_set.dimension, list(data_set.size), None)
element_type = data_set.element_type
dataset_id = H5Dcreate(hdf5_container, data_set.name, ELEMENT_TYPE_TO_HDF5_TYPE[element_type], dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)
data = array(data_set.elements, ELEMENT_TYPE_TO_JARRAY_TYPE[element_type])
type_id = H5Dget_type(dataset_id)
H5Dwrite(dataset_id, type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, data)
H5Dclose(dataset_id)
H5Sclose(dataspace_id)
def _save_hdf5_group_contents(hdf5_container, group):
for data_set in group.datasets.itervalues():
_save_hdf5_data_set(hdf5_container, data_set)
for data_set in group.groups.itervalues():
_save_hdf5_group(hdf5_container, group)
def _save_hdf5_group(hdf5_container, group):
"""
:param hid_t hdf5_container: hdf5 file handle open in write mode or hdf5 group handle
:param Group group:
"""
group_id1 = H5Gcreate(hdf5_container, group.name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)
_save_hdf5_group_contents(hdf5_container, group)
H5Gclose(group_id1)
def save_hdf5_file(filepath, group):
"""
:param str filepath:
:param Group group:
"""
hdf5_file = H5Fcreate(filepath, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)
# _save_hdf5_group(hdf5_file, group)
_save_hdf5_group_contents(hdf5_file, group)
H5Fclose(hdf5_file)

View File

@ -15,10 +15,10 @@ from catalog import Sequence, ImageCatalog
import telemos
from imageengine import IImageEngine, IImageProcessingDebugger, NullDebugger, StackImageFeeder, IHyperStack
# greeting = "Hello, " + name + "!"
from hdf5.hdf5_data import Group, DataSet, ElementType
from imagej.hdf5serializer import save_hdf5_file
import abc
# h5py is not available in fiji's jython
from ncsa.hdf.hdf5lib import HDFArray # found in FIJI_HOME/jars/jhdf5-14.12.6.jar
ABC = abc.ABCMeta('ABC', (object,), {})
@ -166,20 +166,6 @@ class GlobulesAreaEstimator(object):
""" An image processing suite targetted to visible images of traps with moving particles (globules)
"""
class FrameResult(object):
""" The results related to a frame
"""
def __init__(self, globules_area_ratio):
"""
:param int globules_area_ratio: area of globules in this frame, expressed as a ratio of the frame area
"""
self.globules_area_ratio = globules_area_ratio
class Results(object):
def __init__(self):
self.frames = {}
def __init__(self, background_estimator, particle_threshold):
"""
:param IBackgroundEstimator background_estimator: the method that is used to estimate the background image.
@ -190,13 +176,19 @@ class GlobulesAreaEstimator(object):
def detect_particles(self, visible_traps_sequence):
"""
:param IHyperStack visible_traps_sequence:
:param IHyperStack visible_traps_sequence:
:rtype hdf5_data.Group:
"""
background_image = self.background_estimator.estimate_background(visible_traps_sequence)
ie = IImageEngine.get_instance()
# particles_stack = ie.create_hyperstack(width=background_image.width, height=background_image.height, num_slices=1, num_frames=visible_traps_sequence.num_frames(), num_channels=1, pixel_type=PixelType.F32)
results = Group('')
globules_area = DataSet('globuls_area_ratio', (visible_traps_sequence.num_frames(), ), ElementType.F32)
frame_indices = DataSet('frame index', (visible_traps_sequence.num_frames(), ), ElementType.U32)
results.add_dataset(globules_area)
results.add_dataset(frame_indices)
results = GlobulesAreaEstimator.Results()
frame_area = background_image.get_width() * background_image.get_height()
for frame_index in range(visible_traps_sequence.num_frames()):
particle_image = ie.subtract(visible_traps_sequence.get_image(frame_index=frame_index), background_image)
@ -208,9 +200,11 @@ class GlobulesAreaEstimator(object):
num_particle_pixels = int(measured_mean_value * num_pixels / particle_pixel_value)
print("num_particle_pixels: %d " % num_particle_pixels)
globules_area_ratio = float(num_particle_pixels)/frame_area
results.frames[frame_index] = GlobulesAreaEstimator.FrameResult(globules_area_ratio=globules_area_ratio)
globules_area[(frame_index, )] = globules_area_ratio
frame_indices[(frame_index, )] = frame_index
# results.frames[frame_index] = GlobulesAreaEstimator.FrameResult(globules_area_ratio=globules_area_ratio)
# particles_stack.set_image(frame_index=frame_index, image=particle_image)
return results
def test_find_white():
@ -229,7 +223,7 @@ def test_find_white():
src_image_file_path = sequence.get_image_file_path(channel_index=channel_index, frame_index=0, z_index=0)
src_image = IJ.openImage(src_image_file_path)
# src_image = IJ.openImage(src_image_file_path)
depth_index = lipase.find_depth_index(src_image, white_seq)
depth_index = lipase.find_depth_index(src_image, white_seq) # pylint: disable=unused-variable
def run_script():

View File

@ -16,7 +16,7 @@ from lipase.traps_detector import TrapsDetector # pylint: disable=import-error
from lipase.catalog import ImageCatalog, Sequence # pylint: disable=import-error
from lipase.lipase import Lipase, ImageLogger # pylint: disable=import-error
from lipase.lipase import GlobulesAreaEstimator, EmptyFrameBackgroundEstimator # pylint: disable=import-error
from lipase.imagej.hdf5serializer import save_hdf5_file # pylint: disable=import-error
class TestLipase(unittest.TestCase):
@ -96,12 +96,16 @@ class TestLipase(unittest.TestCase):
self.assertAlmostEqual(measured_mean_value, expected_mean_value, delta=0.01)
def test_visible_traps_sequence_processing(self):
traps_sequence = self.catalog.sequences['res_soleil2018/GGH/GGH_2018_cin2_phiG_I_327_vis_-40_1/Pos0']
visible_traps_sequence = traps_sequence.as_hyperstack(['DM300_nofilter_vis'])
background_estimator = EmptyFrameBackgroundEstimator(empty_frame_index=39)
processor = GlobulesAreaEstimator(background_estimator=background_estimator, particle_threshold=2000.0)
processor.detect_particles(visible_traps_sequence)
results = processor.detect_particles(visible_traps_sequence)
save_hdf5_file('results.h5', results)
# results file could be checked with "h5dump --xml ./lipase.git/results.h5"
first_frame_measured_ratio = results['globuls_area_ratio'][(0,)]
first_frame_expected_ratio = 0.008
self.assertAlmostEqual(first_frame_measured_ratio, first_frame_expected_ratio, delta=0.01)
# def test_lipase_process(self):
# lipase = Lipase(self.catalog, debugger=NullDebugger())