added thd io / added cera to json

example/convert_projections.py

+from rq_controller.common import PyProjection
+from rq_controller.common.io.rq_json import RqJsonWriter
+from rq_controller.common.io.thd import RawLoader
+from pathlib import Path
+
+
+FOLDER = Path(r'C:\Users\swittl\Downloads\defect_detection\defect_detection\object_2')
+SAVE_FOLDER = Path(r'C:\Users\swittl\Downloads\defect_detection\defect_detection\rq_converted\object_2')
+
+def main():
+
+
+    writer = RqJsonWriter()
+    loader = RawLoader()
+
+    projections_files = FOLDER.glob('*.raw')
+    
+    for i, file in enumerate(projections_files):
+        
+        projection = loader.load_projection(file)
+        save_path = writer.get_next_projection_save_path(SAVE_FOLDER)
+        writer.write_projection(save_path, projection, optional=True)
+
+if __name__ == '__main__':
+    main()
\ No newline at end of file

rq_controller/common/io/rq_json/json_write.py

@@ -28,7 +28,7 @@ class RqJsonWriter(BaseDataWriter):
 
         self.write_json(save_path, data_dict)
 
-    def write_projection(self, save_path: Path, projection: PyProjection):
+    def write_projection(self, save_path: Path, projection: PyProjection, optional: bool = False):
         save_path_projection_geometry = save_path.parent / f'{save_path.stem}{self.porjection_geometry_suffix}'
         data_dict = dict()
 
@@ -50,6 +50,10 @@ class RqJsonWriter(BaseDataWriter):
         data_dict['current_ua'] = projection.current_ua
         data_dict['frame_id'] = projection.frame_id
 
+        # extra / converted information
+        if optional:
+            data_dict['optional'] = dict(projection_matrix_cera=projection.projection_matrix_cera_px.tolist())
+
         self.write_json(save_path_projection_geometry, data_dict)
 
         image = Image.fromarray(projection.image)

rq_controller/common/io/thd/__init__.py

+from .thd_load import RawLoader
+from .thd_write import RawWriter
\ No newline at end of file

rq_controller/common/io/thd/thd_load.py

+import numpy as np
+import json 
+from pathlib import Path
+
+try:
+    from PythonTools.raw2py import raw2py
+    from PythonTools.ezrt_header import EzrtHeader
+    from PythonTools.rek2py import rek2py
+except ModuleNotFoundError:
+    raise ModuleNotFoundError('Install PythonTools from Fraunhofer EZRT.')
+
+from ..loader import BaseDataLoader, PyProjection, PyProjectionGeometry, PyRegionOfIntrest, PyVolume
+from PIL import Image
+import pyometiff
+from scipy.spatial.transform import Rotation
+
+from ..rq_json.json_load import RqJsonLoader
+
+
+class RawLoader(BaseDataLoader):
+    def __init__(self):
+        super().__init__('.raw', '.raw', '.roi-json', '.rek')
+
+
+    def load_projection_geometry(self, load_path: Path) -> PyProjectionGeometry:
+        header = EzrtHeader.fromfile(load_path)
+
+        focal_spot_mm = np.array(header.agv_source_position) * 1000.
+        detector_postion_mm = np.array(header.agv_detector_center_position) * 1000.
+        
+        line = np.array(header.agv_detector_line_direction)
+        column = np.array(header.agv_detector_col_direction)
+        normal = np.cross(line, column)
+
+        matrix = np.eye(3)
+        matrix[:, 0] = line
+        matrix[:, 1] = column
+        matrix[:, 2] = normal
+
+        rotation = Rotation.from_matrix(matrix)
+        detector_orientation_quad = rotation.as_quat()
+
+        frame_id = header.number_of_images
+
+
+        projection_geometry = PyProjectionGeometry(
+            focal_spot_mm,
+            detector_postion_mm,
+            detector_orientation_quad,
+            frame_id,
+            focal_spot_orientation_quad=np.array([0., 0., 0., 1.]))
+        
+        return projection_geometry
+    
+    def load_projection(self, load_path: Path, switch_order: bool = True) -> PyProjection:
+        header, image = raw2py(load_path, switch_order=switch_order)
+
+        focal_spot_mm = np.array(header.agv_source_position) * 1000.
+        detector_postion_mm = np.array(header.agv_detector_center_position) * 1000.
+        
+        line = np.array(header.agv_detector_line_direction)
+        column = np.array(header.agv_detector_col_direction)
+        normal = np.cross(line, column)
+
+        matrix = np.eye(3)
+        matrix[:, 0] = line
+        matrix[:, 1] = column
+        matrix[:, 2] = normal
+
+        rotation = Rotation.from_matrix(matrix)
+        detector_orientation_quad = rotation.as_quat()
+
+        frame_id = header.number_of_images
+
+        detector_heigth_mm = header.pixel_width_in_um / 1000. * image.shape[1]
+        detector_width_mm = header.pixel_width_in_um / 1000. * image.shape[0]
+
+        voltage_kv = header.voltage_in_kv
+        current_ua = header.current_in_ua
+        exposure_time_ms = header.exposure_time_in_ms
+
+
+        projection = PyProjection(
+            focal_spot_mm,
+            detector_postion_mm,
+            detector_orientation_quad,
+            image,
+            detector_heigth_mm,
+            detector_width_mm,
+            frame_id,
+            np.array([0., 0., 0., 1.]),
+            voltage_kv,
+            current_ua,
+            exposure_time_ms)
+        
+        return projection
+    
+    def load_json(self, load_path: Path) -> dict:
+        with open(str(load_path), 'r') as f:
+            data_dict = json.load(f)
+        return data_dict
+    
+    def load_region_of_intrest(self, load_path: Path) -> PyRegionOfIntrest:
+        RqJsonLoader.load_region_of_intrest(self, load_path)
+    
+    def load_volume(self, load_path: Path) -> PyVolume:
+
+        header, volume = rek2py(load_path)
+        load_path_roi = load_path.parent / f'{load_path.stem}{self.region_of_intrest_suffix}'
+        roi = self.load_region_of_intrest(load_path_roi)
+
+        if volume.dtype == np.uint16:
+            data_type = 0
+        elif volume.dtype == np.uint8:
+            data_type = 1
+        else:
+            raise ValueError('data type not implemented.')
+        
+        return PyVolume(volume, roi, data_type)
+

rq_controller/common/io/thd/thd_write.py

+import numpy as np
+from pathlib import Path
+
+from ..writer import BaseDataWriter, PyProjection, PyProjectionGeometry, PyRegionOfIntrest, PyVolume
+
+try:
+    from PythonTools.py2raw import py2raw
+    from PythonTools.ezrt_header import EzrtHeader
+    from PythonTools.py2rek import py2rek
+except ModuleNotFoundError:
+    raise ModuleNotFoundError('Install PythonTools from Fraunhofer EZRT.')
+
+
+from ..rq_json.json_write import RqJsonWriter
+
+
+class RawWriter(BaseDataWriter):
+    def __init__(self):
+        super().__init__('.raw', '.raw', '.roi-json', '.rek')
+
+    
+
+    def write_projection_geometry(self, save_path: Path, projection_geometry: PyProjectionGeometry):
+        header = EzrtHeader()
+        header.agv_source_position = (projection_geometry.focal_spot_mm / 1000.).tolist()
+        header.agv_detector_center_position = (projection_geometry.detector_postion_mm / 1000.).tolist()
+
+        line = projection_geometry.detector_rotation_matrix[:, 0]
+        column = projection_geometry.detector_rotation_matrix[:, 0]
+
+        header.agv_detector_line_direction = line
+        header.agv_detector_col_direction = column
+
+        header.number_of_images = projection_geometry.frame_id
+        header.tofile(save_path)
+
+    def write_projection(self, save_path: Path, projection: PyProjection, swith_order: bool = False):
+        projection_geometry = projection.get_projection_geometry()
+
+        header = EzrtHeader()
+        
+        header.agv_source_position = (projection_geometry.focal_spot_mm / 1000.).tolist()
+        header.agv_detector_center_position = (projection_geometry.detector_postion_mm / 1000.).tolist()
+
+        line = projection_geometry.detector_rotation_matrix[:, 0]
+        column = projection_geometry.detector_rotation_matrix[:, 0]
+
+        header.agv_detector_line_direction = line
+        header.agv_detector_col_direction = column
+
+        header.number_of_images = projection_geometry.frame_id
+        header.pixel_width_in_um = projection.pixel_pitch_x_mm * 1000.
+        header.voltage_in_kv = projection.voltage_kv
+        header.current_in_ua = projection.current_ua
+        header.exposure_time_in_ms = projection.exposure_time_ms
+
+        py2raw(projection.image, save_path, header, swith_order)
+
+        
+    def write_region_of_intrest(self, save_path: Path, region_of_intrest: PyRegionOfIntrest):
+        RqJsonWriter.write_region_of_intrest(self, save_path, region_of_intrest)
+        
+    def write_volume(self, save_path: Path, volume: PyVolume):
+        header = EzrtHeader()
+        header.num_voxel_x = volume.shape[0]
+        header.num_voxel_y = volume.shape[1]
+        header.num_voxel_z = volume.shape[2]
+        header.pixel_width_in_um = volume.roi.resolution_mm[0]
+
+        py2rek(volume.array, save_path, header, True)
\ No newline at end of file

rq_controller/common/projection.py

@@ -254,4 +254,21 @@ class PyProjection(PyProjectionGeometry):
     def pixel_pitch_y_mm(self) -> float:
         return self.detector_heigth_mm / self.detector_heigth_px
     
+    @property
+    def to_cera_transformation_matrix(self) -> np.ndarray:
+        transformation_matrix = np.eye(3)
+        # [mm] -> [px]
+        transformation_matrix[2, 2] = 1. / self.pixel_pitch_x_mm 
+        transformation_matrix[0, 2] = -self.detector_width_px / 2.
+        transformation_matrix[1, 2] = -self.detector_heigth_px / 2.
+
+        return transformation_matrix
+
+    
+    @property
+    def projection_matrix_cera_px(self):
+        projection_geometry = self.get_projection_geometry()
+        return np.linalg.inv(self.to_cera_transformation_matrix) @ projection_geometry.projection_matrix
+
+    
     

rq_controller/common/projection_geometry.py

@@ -3,6 +3,7 @@ from __future__ import annotations
 import numpy as np
 
 from rq_interfaces.msg import ProjectionGeometry, Projection
+from scipy.spatial.transform import Rotation
 
 class PyProjectionGeometry():
     """
@@ -112,3 +113,27 @@ class PyProjectionGeometry():
         message.header.frame_id = self.frame_id
 
         return message
+    
+    @property
+    def detector_rotation_matrix(self) -> np.ndarray:
+        return Rotation.from_quat(self.detector_orientation_quad).as_matrix()
+
+    
+    @property 
+    def detector_horizontal_vector(self) -> np.ndarray:
+        return self.detector_rotation_matrix[:, 0]
+    
+    @property 
+    def detector_vertical_vector(self) -> np.ndarray:
+        return self.detector_rotation_matrix[:, 1]
+
+    
+    @property
+    def projection_matrix(self) -> np.ndarray:
+        p3x3 = np.vstack([self.detector_horizontal_vector,
+                          self.detector_vertical_vector, 
+                          self.detector_postion_mm - self.focal_spot_mm]).T
+        p3x3_inv = np.linalg.inv(p3x3)
+        p4 = (p3x3_inv @ (-self.focal_spot_mm)).reshape((3, 1))
+        matrix = np.concatenate([p3x3_inv, p4], 1).reshape((3, 4))
+        return matrix