diff --git a/README.md b/README.md index e822939..be329d7 100644 --- a/README.md +++ b/README.md @@ -28,3 +28,5 @@ bob@stykades:~/work/grassloper$ source ./grassloper.venv/bin/activate ```sh (grassloper.venv) bob@stykades:~/work/grassloper$ grassloper --help ``` + +`grassloper` uses the open source [hdf5 file format](https://en.wikipedia.org/wiki/Hierarchical_Data_Format#HDF5) to store its data. These files can be explored by [hdfview](https://www.hdfgroup.org/downloads/hdfview/) or [HDFCompass](https://support.hdfgroup.org/projects/compass/) diff --git a/grassloper/main.py b/grassloper/main.py index 07b1760..d6180b8 100755 --- a/grassloper/main.py +++ b/grassloper/main.py @@ -1,4 +1,4 @@ -##!/usr/bin/env python3 +#!/usr/bin/env python3 import numpy as np import cv2 from pathlib import Path @@ -8,7 +8,9 @@ import time import scipy.stats import argparse from improtools.imageprocessing import ImageProcessDebugger, IMovieProcessListener, NullMovieProcessListener, MovieProcessDebugger +from typing import Dict, Tuple +GRASSLOPER_VERSION = '1.0' def create_slope_image(image_width: int, image_height: int, p1_angle: float, p1_radius: float, bg_color: float = 0.0, fg_color: float = 1.0): """ @@ -69,7 +71,7 @@ def hdf5_to_trac_data(hdf5_file_path: Path): version = f.attrs['version'] # noqa date = f.attrs['date'] # noqa - trac_data.num_frames = f.attrs['nFrames'] + trac_data.num_frames = f.attrs['nFrames'] - 2 # it seems that tractrac doesn't output particles for 2 frames trac_data.image_size = f.attrs['size'] # print(type(trac_data.image_size)) # f.attrs['size'] = I0f.shape @@ -142,12 +144,13 @@ class SlopeFinder(): x = [pt[0] for pt in surface_pts] y = [pt[1] for pt in surface_pts] lin_regress_result = scipy.stats.linregress(x, y) - print(lin_regress_result) + # print(type(lin_regress_result)) from_vertex = [0.0, lin_regress_result.intercept] x_max = float(isbead_image.shape[1]) to_vertex = [x_max, lin_regress_result.intercept + x_max * lin_regress_result.slope] mo_pro_listener.m_imageProcessListener.onLine( ((from_vertex[0], from_vertex[1]), (to_vertex[0], to_vertex[1])), 'slope') - + return lin_regress_result + # Line(lin_regress_result.slope, lin_regress_result.intercept) @staticmethod def create_isbead_image(trac_data: TracData, frame_index: int, particle_radius: float): @@ -162,9 +165,40 @@ class SlopeFinder(): return image +def save_results(results_file_path: Path, fitted_lines: Dict[int, scipy.stats._stats_mstats_common.LinregressResult], input_file_path: Path, num_frames: int, images_size: Tuple[int, int]): + frame_index_array = np.zeros(shape=(num_frames), dtype=np.uint16) + slope_array = np.zeros(shape=(num_frames), dtype=np.float32) + intercept_array = np.zeros(shape=(num_frames), dtype=np.float32) + stderr_array = np.zeros(shape=(num_frames), dtype=np.float32) + intercept_stderr_array = np.zeros(shape=(num_frames), dtype=np.float32) + + for frame_index in range(num_frames): + frame_index_array[frame_index] = frame_index + 1 + if frame_index_array[frame_index] in fitted_lines.keys(): + fitted_line = fitted_lines[frame_index_array[frame_index]] + slope_array[frame_index] = fitted_line.slope + intercept_array[frame_index] = fitted_line.intercept + stderr_array[frame_index] = fitted_line.stderr + intercept_stderr_array[frame_index] = fitted_line.intercept_stderr + + f = h5py.File(results_file_path,'w') + f.attrs['version']='HDF5 file made with grassloper v %s' % GRASSLOPER_VERSION + f.attrs['date'] = time.strftime("%d/%m/%Y") + f.attrs['input_file'] = str(input_file_path) + f.attrs['num_frames'] = num_frames + f.attrs['image_size'] = images_size + f.create_dataset("frame_index", data=frame_index_array) + f.create_dataset("slope", data=slope_array) + f.create_dataset("intercept", data=intercept_array) + f.create_dataset("stderr", data=stderr_array) + f.create_dataset("intercept_stderr", data=intercept_stderr_array) + f.close() + + def main(): - parser = argparse.ArgumentParser('grassloper', description='estimates the slope of granular surface', epilog='example: \n\tgrassloper --part-pos-file=./grassloper.git/samples/TracTrac/sample002_track.hdf5 --particle-radius=11.0 --debug-dir=\'./improdebug\'') + parser = argparse.ArgumentParser('grassloper', description='estimates the slope of granular surface', epilog='example: \n\tgrassloper --part-pos-file=./grassloper.git/samples/TracTrac/sample002_track.hdf5 --particle-radius=11.0 --debug-dir=\'./improdebug\' --results-file ./grassloper.git/samples/TracTrac/sample002_track_slopes.hdf5') parser.add_argument('--part-pos-file', type=Path, default=None, required=True, help='the input file containing the particle positions') + parser.add_argument('--results-file', type=Path, default=None, required=True, help='the output file containing for each frame one line fitting the granular surface (hdf5 file format)') parser.add_argument('--part-pos-file-format', type=str, choices=['tractrac-hdf5'], default='tractrac-hdf5', help='the file format of the input file that contains particle positions') parser.add_argument('--particle-radius', type=float, required=True, help='the radius of the particles in pixels', default=11) parser.add_argument('--debug-dir', type=Path, default=None, help='where to store image processing debug files') @@ -176,11 +210,16 @@ def main(): mo_pro_listener = NullMovieProcessListener() trac_data = hdf5_to_trac_data(args.part_pos_file) slope_finder = SlopeFinder(beads_radius=args.particle_radius) - for frame_index in range(1, 3): + print(trac_data.num_frames) + fitted_lines = {} + for frame_index in range(1, trac_data.num_frames + 1): mo_pro_listener.onImageProcessingStart() - slope_finder.find_slope(trac_data, frame_index=frame_index, mo_pro_listener=mo_pro_listener) + line = slope_finder.find_slope(trac_data, frame_index=frame_index, mo_pro_listener=mo_pro_listener) + print('frame %d: %s' % (frame_index, str(line))) + fitted_lines[frame_index] = line mo_pro_listener.onImageProcessingEnd() + save_results(args.results_file, fitted_lines, args.part_pos_file, trac_data.num_frames, trac_data.image_size) # slope_image = create_slope_image(image_width=128, image_height=256, p1_angle=0.7, p1_radius=3.0) # cv2.imwrite(str(Path('toto.tif')), slope_image)