.gitignore

@@ -3,4 +3,5 @@
 /dist/
 __pycache__
 *.egg-info
-*temp.*
\ No newline at end of file
+*temp.*
+/dir2025/
\ No newline at end of file

pyproject.toml

 [project]
 name = "cyxtrax"
-version = "0.1.0"
+version = "0.1.1"
 description = "Add your description here"
 readme = "README.md"
 authors = [
@@ -14,8 +14,7 @@ dependencies = [
     "jax>=0.4.35 ; sys_platform == 'win32'",
     "jax[cuda12]>=0.4.35 ; sys_platform == 'linux'",
     "kiwisolver>=1.4.7",
-    "matplotlib>=3.9.2",
-    "numpy>=2.1.3",
+    "numpy>=2.0.2",
     "ruff>=0.7.3",
     "scipy>=1.14.1",
     "xraylib>=4.1.5",

scripts/02_generate_atlas.py

 from cyxtrax.io import record_atlas
 from cyxtrax.simulation import CylindricalProjection, load_mesh
-from cyxtrax.data import olaf_mesh
+from cyxtrax.data import olaf_mesh, MeshObject
 from pathlib import Path
+import numpy as np
 
 
 # Just get the temp folder to store the atlas
@@ -13,16 +14,18 @@ def main():
     # code_embedder:A start
     # Create the cylindrical projection interface
     cyl_proj = CylindricalProjection()
-    load_mesh(cyl_proj, [olaf_mesh])
+    obj = MeshObject(Path(r'C:\Program Files\BAM\aRTist 2.12\Data\Library\ExampleParts\Fun\Smurf.stl'),
+                     np.array([0., 0., 0.]), np.array([0., 0., 0., 1.]))
+    load_mesh(cyl_proj, [obj])
     # make an atlas
 
     save_path = record_atlas(
         cyl_proj,
-        [olaf_mesh],
+        [obj],
         "test",
         TEMP_FOLDER,
         map_bounding_box=[-10, 20],
-        number_of_maps=2,
+        number_of_maps=1,
     )
     # code_embedder:A end
     print(f"Atlas saved to the path: {save_path}")

scripts/03_load_atlas.py

@@ -5,12 +5,12 @@ from pathlib import Path
 
 # Just get the temp folder to store the atlas
 FOLDER = Path(__file__)
-TEMP_FOLDER = FOLDER.parent.parent / "temp"
+TEMP_FOLDER = FOLDER.parent.parent / "dir2025"
 
 
 def main():
     # code_embedder:A start
-    files = list(TEMP_FOLDER.glob("*.h5"))
+    files = list(TEMP_FOLDER.glob("*1_smurf.h5"))
     if not files:
         raise FileNotFoundError("No .h5 files in temp folder. Run Example 02!")
 

scripts/04_mapping.py

@@ -10,6 +10,7 @@ from matplotlib import pyplot as plt, figure
 # Just get the temp folder to store the atlas
 FOLDER = Path("/mnt/data/xray_studies_data/")
 TEMP_FOLDER = FOLDER / "dataset_64"
+TEMP_FOLDER = Path('temp')
 
 
 def main():
@@ -24,9 +25,9 @@ def main():
     print(f"Cone Beam Mode set: {cyl_proj.cone_mode}")
 
     rot = 30
-    x, y = 2000 * jnp.cos(rot * jnp.pi / 180), 500 * jnp.sin(rot * jnp.pi / 180)
+    x, y = 500 * jnp.cos(rot * jnp.pi / 180), 500 * jnp.sin(rot * jnp.pi / 180)
     source_position = jnp.array([x, y, 0.0]).reshape((1, 3))
-    detector_position = jnp.array([-x, -y, -10.0]).reshape((1, 3))
+    detector_position = jnp.array([-x, -y, 0.0]).reshape((1, 3)) / 2.
     detctor_orientation_scipy = Rotation.from_euler("xyz", [0, -90, rot], degrees=True)
 
     temp_projection = Path("temp.tif")
@@ -39,7 +40,7 @@ def main():
         "D", detector_position[0, 0], detector_position[0, 1], detector_position[0, 2]
     )
     cyl_proj.api.rotate_from_quat("D", detctor_orientation[0])
-    cyl_proj.api.save_image(temp_projection, save_mode=SAVEMODES.FLOAT_TIFF)
+    cyl_proj.api.save_image(temp_projection, save_mode=SAVEMODES.FLOAT_TIFF, save_projection_geometry=1)
 
     projection, geometry = utility.load_projection(temp_projection)
 
@@ -88,11 +89,12 @@ def main():
     ax2 = subfig2.add_subplot(111)
     ax2.set_title("Cylindrical Mapping")
     ax2.axis(False)
-    for _ in range(3):
-        for _ in range(3):
+    nums_subplots = 1
+    for _ in range(nums_subplots):
+        for _ in range(nums_subplots):
             if counter >= number_of_maps:
                 continue
-            ax = subfig2.add_subplot(3, 3, counter + 1)
+            ax = subfig2.add_subplot(nums_subplots, nums_subplots, counter + 1)
             ax.imshow(maps[:, :, counter], vmin=maps.min(), vmax=maps.max())
             ax.scatter(
                 angles_calc[0, counter], angles_calc[1, counter], c=colors[counter]

scripts/07_presentation.py

@@ -22,11 +22,11 @@ MAP_FOLDER = Path(r"C:\data\XrayTransformLow")
 
 
 def main():
-    files = list(MAP_FOLDER.glob("*.h5"))
+    files = list(TEMP_FOLDER.glob("*.h5"))
     if not files:
         raise FileNotFoundError("No .h5 files in temp folder. Run Example 02!")
 
-    maps, points, mesh_object_list = load_atlas(files[3])
+    maps, points, mesh_object_list = load_atlas(files[-1])
     cyl_proj = CylindricalProjection()
 
     load_mesh(cyl_proj, mesh_object_list, False)

scripts/09_compute_smurf_map.py

+from cyxtrax.io import record_map
+from cyxtrax.simulation import CylindricalProjection, load_mesh
+from cyxtrax.data import olaf_mesh, MeshObject
+from pathlib import Path
+import numpy as np
+
+
+# Just get the temp folder to store the atlas
+FOLDER = Path(__file__)
+TEMP_FOLDER = FOLDER.parent.parent / "temp"
+
+
+def main():
+    # code_embedder:A start
+    # Create the cylindrical projection interface
+    cyl_proj = CylindricalProjection()
+    obj = MeshObject(Path(r'C:\Program Files\BAM\aRTist 2.12\Data\Library\ExampleParts\Fun\Smurf.stl'),
+                     np.array([0., 0., 0.]), np.array([0., 0., 0., 1.]))
+    load_mesh(cyl_proj, [obj])
+    # make an atlas
+
+    save_path = record_map(
+        cyl_proj,
+        [obj],
+        "test",
+        TEMP_FOLDER,
+        position=np.array([56., 10.5, 2.])
+    )
+    # code_embedder:A end
+    print(f"Atlas saved to the path: {save_path}")
+
+
+if __name__ == "__main__":
+    main()

scripts/10_dir2025.py

+from cyxtrax.io import record_atlas
+from cyxtrax.simulation import CylindricalProjection, load_mesh
+from cyxtrax.data import olaf_mesh, MeshObject
+from pathlib import Path
+import numpy as np
+
+
+# Just get the temp folder to store the atlas
+FOLDER = Path(__file__)
+TEMP_FOLDER = FOLDER.parent.parent / "dir2025"
+
+
+def main():
+    # code_embedder:A start
+    # Create the cylindrical projection interface
+    cyl_proj = CylindricalProjection()
+    obj = MeshObject(Path(r'C:\Program Files\BAM\aRTist 2.12\Data\Library\ExampleParts\Fun\Smurf.stl'),
+                     np.array([0., 0., 0.]), np.array([0., 0., 0., 1.]))
+    load_mesh(cyl_proj, [obj])
+    # make an atlas
+
+    save_path = record_atlas(
+        cyl_proj,
+        [obj],
+        "smurf",
+        TEMP_FOLDER,
+        map_bounding_box=[-50, 50],
+        number_of_maps=4,
+    )
+    # code_embedder:A end
+    print(f"Atlas saved to the path: {save_path}")
+
+
+if __name__ == "__main__":
+    main()

scripts/10_dir2025_2.py

+from cyxtrax.io import record_atlas
+from cyxtrax.simulation import CylindricalProjection, load_mesh
+from cyxtrax.data import olaf_mesh, MeshObject
+from pathlib import Path
+import numpy as np
+
+
+# Just get the temp folder to store the atlas
+FOLDER = Path(__file__)
+TEMP_FOLDER = FOLDER.parent.parent / "dir2025"
+
+
+def main():
+    # code_embedder:A start
+    # Create the cylindrical projection interface
+    cyl_proj = CylindricalProjection()
+    obj = MeshObject(Path(r'C:\Program Files\BAM\aRTist 2.12\Data\Library\ExampleParts\Fun\Smurf.stl'),
+                     position_mm=-np.array([56., 10.5, 2.]), orientation_quat=np.array([0., 0., 0., 1.]))
+    load_mesh(cyl_proj, [obj])
+    # make an atlas
+
+    save_path = record_atlas(
+        cyl_proj,
+        [obj],
+        "smurf",
+        TEMP_FOLDER,
+        map_bounding_box=[-2, 2],
+        number_of_maps=3,
+    )
+    # code_embedder:A end
+    print(f"Atlas saved to the path: {save_path}")
+
+
+if __name__ == "__main__":
+    main()

src/cyxtrax/io/__init__.py

 from .generate_atlas import record_atlas
 from .load_maps import load_atlas
 from .save_map import save_atlas
+from .generate_single_map import record_map
 
 
-__all__ = ["record_atlas", "load_atlas", "save_atlas"]
+__all__ = ["record_atlas", "load_atlas", "save_atlas", "record_map"]

src/cyxtrax/io/generate_single_map.py

+import numpy as np
+import h5py
+from cyxtrax.simulation.artist_bridge import CylindricalProjection
+from cyxtrax.common.mesh_object import MeshObject
+from pathlib import Path
+import json
+
+
+def to_dict(mesh_object=MeshObject) -> dict:
+    return mesh_object.as_dict()
+
+
+def record_map(
+    cylindrical_projection: CylindricalProjection,
+    mesh_list: list[MeshObject],
+    atlas_name: str,
+    save_folder: Path,
+    position: np.ndarray = np.zeros((3, 1))
+) -> Path:
+    file_index = len(list(save_folder.glob("*.h5"))) + 1
+    save_path = save_folder / f"{file_index:05}_{atlas_name}.h5"
+    file = h5py.File(save_path, "a")
+    projection = file.require_dataset(
+        "maps",
+        shape=(
+            cylindrical_projection.x_px,
+            cylindrical_projection.y_px,
+            1,
+        ),  # Initial shape with third dimension as 0
+        maxshape=(
+            cylindrical_projection.x_px,
+            cylindrical_projection.y_px,
+            1,
+        ),
+        dtype=np.float32,
+    )
+
+    projection_points = file.require_dataset(
+        "positions",
+        shape=(3, 1),  # Initial shape with third dimension as 0
+        maxshape=(3, 1),
+        dtype=np.float32,
+    )
+
+    mesh_list_dict = list(map(to_dict, mesh_list))
+    print(mesh_list_dict)
+
+    file.attrs["mesh_list"] = json.dumps(mesh_list_dict)
+
+    counter = 0
+
+    image = cylindrical_projection.compute_projection(
+        position, output_full_ray_projection=True
+    )
+    # current_size = projection.shape[2]
+    # new_size = current_size + 1
+    # projection.resize((cylindrical_projection.x_px, cylindrical_projection.y_px, new_size))
+    projection[:, :, counter : counter + 1] = image[
+        :, :, np.newaxis
+    ].astype(np.float32)
+
+    # projection_points.resize((3, new_size))
+    projection_points[:, counter : counter + 1] = position.reshape((3, 1))
+    counter += 1
+    return save_path

src/cyxtrax/io/save_map.py

@@ -52,3 +52,22 @@ def save_atlas(
 
     projection[:] = maps
     projection_points[:] = map_position
+
+
+def add_rays_to_atlas(save_path: Path, ray_start: np.ndarray, ray_end: np.ndarray):
+    file = h5py.File(save_path, "a")
+    ray_endpoint = file.require_dataset(
+        "ray_endpoint",
+        shape=ray_end.shape,
+        dtype=np.float32,
+    )
+    ray_endpoint[:] = ray_end
+
+    ray_startpoint = file.require_dataset(
+        "ray_startpoint",
+        shape=ray_end.shape,
+        dtype=np.float32,
+    )
+    ray_startpoint[:] = ray_start
+
+

src/cyxtrax/simulation/artist/artist_bridge.py

+try:
+    from artist_pythonlib import API, utility, SAVEMODES  # type: ignore
+except ModuleNotFoundError:
+    from warnings import warn
+
+    warn(
+        "The module `artistlib`is not installed. The simulation module is not 100\% usable! \nInstall: https://github.com/wittlsn/aRTist-PythonLib"
+    )
+
+    def API():
+        return None
+
+    utility = None
+    SAVEMODES = None
+from pathlib import Path
+import numpy as np
+import os
+import importlib.resources
+from time import sleep
+from cyxtrax.common.mesh_object import MeshObject
+
+# !!!!!!!!!!!!!!!!!!!
+# Read carefully
+#!!!!!!!!!!!!!!!!?!!!
+
+
+# This Script assumes that the detector has a size of 6283.19 mm / 2 pi
+# and the curvature is set to the y axis with center on the source.
+
+# The source Position has to be on [0, 0, 0]
+# The detector has to be on Position [1000,0, 0] and orientation [0, 90, 0] / curvature y
+
+# The Script assumes this geometry and only moves all the objects!
+# Preview must be reseted!!!
+
+GLOBAL_COUNTER = 0
+
+
+def set_cylindrical_mode(api: API):
+    with importlib.resources.path("cyxtrax.data", "cylinder.aRTist") as file_path:
+        print(f"Load: {file_path}")
+        api.load_project(file_path)
+
+
+def set_cone_mode(api: API):
+    with importlib.resources.path("cyxtrax.data", "cone.aRTist") as file_path:
+        print(f"Load: {file_path}")
+        api.load_project(file_path)
+
+def set_cone_poly_mode(api: API):
+    with importlib.resources.path("cyxtrax.data", "cone_poly.aRTist") as file_path:
+        print(f"Load: {file_path}")
+        api.load_project(file_path)
+
+
+class CylindricalProjection:
+    def __init__(self, api: API = API()) -> None:
+        self.api = api
+        self.x_resolution_mm = np.pi
+        self.y_resolution_mm = np.pi
+        self.pixel_pitch_mm = 0.139
+        self.x_px = 2000
+        self.y_px = 2000
+        self.detector_radius_mm = 1000.0
+        self.objects = list()
+        self._cylindrical = True
+
+    @property
+    def cylindrical_mode(self) -> bool:
+        return self._cylindrical
+
+    @property
+    def cone_mode(self) -> bool:
+        return not self._cylindrical
+
+    def translate(self, position: np.ndarray):
+        # for obj in self.objects:
+        x_d, y_d, z_d = 1000, 0, 0
+        x_s, y_s, z_s = 0, 0, 0
+        x_p, y_p, z_p = position[0], position[1], position[2]
+
+        self.api.rc.send(
+            f'::PartList::Invoke {str("S")} SetPosition {str(x_s + x_p)} {str(y_s + y_p)} {str(z_s + z_p)};'
+        )
+        self.api.rc.send(
+            f'::PartList::Invoke {str("S")} SetRefPos {str(x_s + x_p)} {str(y_s + y_p)} {str(z_s + z_p)};'
+        )
+
+        self.api.rc.send(
+            f'::PartList::Invoke {str("D")} SetRefPos {str(x_s + x_p)} {str(y_s + y_p)} {str(z_s + z_p)};'
+        )
+        self.api.rc.send(
+            f'::PartList::Invoke {str("D")} SetPosition {str(x_d + x_p)} {str(y_d + y_p)} {str(z_d + z_p)};'
+        )
+        self.api.rc.send("::XDetector::SetDownCurvedView;")
+
+    def compute_projection(
+        self,
+        position: np.ndarray,
+        temp_file_path: Path = Path(r"C:\data"),
+        output_full_ray_projection: bool = True,
+    ) -> np.ndarray:
+        global GLOBAL_COUNTER
+        self.translate(position)
+        temp_file = temp_file_path / f"temp_{GLOBAL_COUNTER:01}.tiff"
+        GLOBAL_COUNTER += 1
+        GLOBAL_COUNTER = GLOBAL_COUNTER % 10
+
+        self.api.save_image(
+            temp_file, save_projection_geometry=False, save_mode=SAVEMODES.FLOAT_TIFF
+        )
+
+        while not temp_file.exists():
+            print("sleepy sleep sleep ...")
+            # The computation of the cylindrical xray tansforms takes time ... aRTist is not programmed for it.
+            # This catches problems with the computation and the remote control interfaces
+            # It is ugly ... but works ...
+            # Also if the mesh is not waterprof the triangle filter of artist seems to have problems.
+
+            sleep(2.0)
+
+        image = utility.load_projection(temp_file, load_projection_geometry=False)[0]
+        os.remove(temp_file)
+
+        if output_full_ray_projection:
+            return self.convert_rays(image)
+        else:
+            return image
+
+    def convert_rays(self, image):
+        inter_image = np.roll(image, 1000, 0)
+        image += np.flip(inter_image, 1)
+        return image
+
+    def set_cylindrical_mode(self):
+        set_cylindrical_mode(self.api)
+        self.x_px = 2000.0
+        self.y_px = 2000.0
+        self.x_resolution_mm = np.pi
+        self.y_resolution_mm = np.pi
+        self.detector_radius_mm = 1000.0
+        self._cylindrical = True
+        print(
+            "Care! The default pixel pitch etc is set. If changes are made change the state of this class ..."
+        )
+
+    def set_cone_mode(self):
+        set_cone_mode(self.api)
+        self.pixel_pitch_mm = 0.139
+        self._cylindrical = False
+        print(
+            "Care! The default pixel pitch etc is set. If changes are made change the state of this class ..."
+        )
+
+
+def load_mesh(
+    cylindrical_projection: CylindricalProjection,
+    object_list: list[MeshObject],
+    cylindrical_mode: bool = True,
+) -> bool:
+    if cylindrical_mode:
+        cylindrical_projection.set_cylindrical_mode()
+    else:
+        cylindrical_projection.set_cone_mode()
+
+    for mesh_object in object_list:
+        object_id = cylindrical_projection.api.load_part(mesh_object.object_path)
+        cylindrical_projection.api.translate(
+            object_id,
+            mesh_object.position_mm[0],
+            mesh_object.position_mm[1],
+            mesh_object.position_mm[2],
+        )
+        cylindrical_projection.api.rotate_from_quat(
+            object_id, mesh_object.orientation_quat
+        )
+
+    return True

src/cyxtrax/util/transform_to_rays.py

+from pathlib import Path
+
+from cyxtrax.io.load_maps import load_atlas
+from cyxtrax.io.save_map import add_rays_to_atlas
+import numpy as np
+
+
+def transform_to_rays(file: Path) -> np.ndarray:
+    maps, points, _ = load_atlas(file)
+    pixel_scale = (maps.shape[0]  - 1.) / maps.shape[0] 
+
+    def cylinder_coordinates(alpha: float, radius: float = 1000.) -> np.ndarray:
+        x = np.sin(alpha) * radius
+        y = np.cos(alpha) * radius
+        point = np.array([x, y, 0.])
+        return point
+    
+    def map_z(point, radius: float = 1000.):
+        point = point.reshape((1, 3))
+        z = np.linspace(-radius, radius, maps.shape[1]) * np.pi * pixel_scale
+        z_array = np.zeros((maps.shape[1], 3))
+        z_array[..., 2] = z
+        return point + z_array
+
+    
+    xyz = list(map(cylinder_coordinates, np.linspace(
+        -np.pi * pixel_scale, np.pi * pixel_scale, maps.shape[0], endpoint=False)))
+    xyz = np.array(list(map(map_z, xyz)))
+
+    start = points.swapaxes(0, 1)
+    
+    ray_end = np.expand_dims(xyz, 2) + np.expand_dims(np.expand_dims(start, 0), 0)
+    ray_end = np.roll(ray_end, 500, 0)
+
+    ray_start = np.roll(ray_end, 1000, 0)
+    ray_start = np.flip(ray_start, 1)
+    return ray_start, ray_end
+
+
+if __name__ == '__main__':
+    save_path = Path(r'.\dir2025\00004_smurf.h5')
+    maps, points, mesh_object_list = load_atlas(save_path)
+
+
+    ray_start, ray_end = transform_to_rays(save_path)
+
+    from cyxtrax.io.save_map import add_rays_to_atlas
+    from cyxtrax.mapping import map_source_2_cylinder
+
+    add_rays_to_atlas(save_path, ray_start, ray_end)
+
+    index_0, index_1, map_index = 1432, 1111, 11
+
+    source = ray_start[index_0, index_1, map_index]
+
+    values, uvw = map_source_2_cylinder(source, maps, points)
+    print(values.shape)
+    print(uvw.shape)
+
+    value = values[map_index]
+    print(uvw[:2, map_index], values[map_index], maps[index_0, index_1, map_index])
+
+    
\ No newline at end of file