- grassloper now saves its results in an hdf5 file in a structured form, in addition to stdout

- grassloper now processes all frames from the input file instead of an arbitrary hardcoded value

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/)

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)