Manual merge

.gitignore

@@ -5,7 +5,6 @@ CMakeCache.txt
 CMakeFiles
 CMakeScripts
 Testing
-Makefile
 cmake_install.cmake
 install_manifest.txt
 compile_commands.json
@@ -16,3 +15,12 @@ src/meow/meow_triggers/*
 src/meow/runner_data/*
 src/meow/runner_output/*
 src/meow/runner_processing/*
+
+# Compiled files
+*.so
+*.so.dSYM
+a.out
+cudafe
+
+# Ignore the $BONE_DATA symlinks, as they're only there for convinience in vscode
+data_*

.vscode/launch.json

@@ -0,0 +1,62 @@
+{
+    // Use IntelliSense to learn about possible attributes.
+    // Hover to view descriptions of existing attributes.
+    // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
+    "version": "0.2.0",
+    "configurations": [
+        {
+            "name": "Python: pre-cleanup/cylinder_surface2",
+            "type": "python",
+            "request": "launch",
+            "program": "${workspaceFolder}/pre-cleanup-src/analysis/cylinder_surface2.py",
+            "console": "integratedTerminal",
+            "args": ["770c_pag"],
+            "justMyCode": false
+        },
+        {
+            "name": "Python: Test geometry",
+            "type": "python",
+            "request": "launch",
+            "program": "${workspaceFolder}/src/test/test_geometry.py",
+            "console": "integratedTerminal",
+            "args": [],
+            "justMyCode": false
+        },
+        {
+            "name": "Python: 0400_h5tobin",
+            "type": "python",
+            "request": "launch",
+            "program": "${workspaceFolder}/src/processing_steps/0400_h5tobin.py",
+            "console": "integratedTerminal",
+            "args": ["770c_pag"],
+            "justMyCode": false
+        },
+        {
+            "name": "Python: 0500_rescale_cupy_bin",
+            "type": "python",
+            "request": "launch",
+            "program": "${workspaceFolder}/src/processing_steps/0500_rescale_cupy_bin.py",
+            "console": "integratedTerminal",
+            "args": ["770c_pag"],
+            "justMyCode": false
+        },
+        {
+            "name": "Python: 0700_implant_FoR",
+            "type": "python",
+            "request": "launch",
+            "program": "${workspaceFolder}/src/processing_steps/0700_implant_FoR.py",
+            "console": "integratedTerminal",
+            "args": ["770c_pag"],
+            "justMyCode": false
+        },
+        {
+            "name": "Python: 0800_implant_data",
+            "type": "python",
+            "request": "launch",
+            "program": "${workspaceFolder}/src/processing_steps/0800_implant_data.py",
+            "console": "integratedTerminal",
+            "args": ["770c_pag"],
+            "justMyCode": false
+        },
+    ]
+}

pre-cleanup-src/Makefile

@@ -0,0 +1,14 @@
+CXX=icpc
+CXXFLAGS += -std=c++17
+CXXFLAGS += $(shell pkg-config --cflags fmt)
+CXXFLAGS += $(shell pkg-config --cflags hdf5)
+LIBS += $(shell pkg-config --libs fmt)
+LIBS += $(shell pkg-config --libs hdf5)
+
+%.o: %.cc
+	$(CXX) $(CXXFLAGS) -c $< 
+
+h5-blockmap: h5-blockmap.o
+	echo $(LIBS)
+	$(CXX) $(CXXFLAGS) $< $(LIBS) -o $@
+

src/analysis/cylinder_surface2.py → pre-cleanup-src/analysis/cylinder_surface2.py

@@ -1,9 +1,9 @@
 #!/usr/bin/env python3
 import os, sys, h5py, numpy as np, pathlib, tqdm, vedo, matplotlib.pyplot as plt, edt, vedo.pointcloud as pc, scipy.ndimage as ndi
-sys.path.append(sys.path[0]+"/../")
+sys.path.append(sys.path[0]+"/../../src")
 from config.paths import *
-from helper_functions import *
-from pybind_kernels.geometry import cylinder_projection
+from lib.py.helpers import commandline_args
+from lib.cpp.cpu_seq.geometry import cylinder_projection
 NA = np.newaxis
 
 
@@ -25,10 +25,10 @@ def homogeneous_transform(xs, M):
 
 def np_save(path,data):
     output_dir = os.path.dirname(path)
-    pathlib.Path(output_dir).mkdir(parents=True, exist_ok=True)        
+    pathlib.Path(output_dir).mkdir(parents=True, exist_ok=True)
     np.save(path,data)
 
-    
+
 # Requires: implant-FoR
 #           soft-tissue/bone segmentation + blood analysis
 #           EDT-field
@@ -55,18 +55,18 @@ def np_save(path,data):
     print(f"Cant't read implant frame-of-reference: {e}")
     print(f"Make sure you have run segment-implant-cc.py and implant-FoR.py for {sample} at scale {mask_scale}x")
     sys.exit(-1)
-    
+
 try:
     blood_mask    = h5mask["blood/mask"][:]
     solid_implant = h5mask["implant_solid/mask"][:]
-    h5mask.close()    
+    h5mask.close()
 except Exception as e:
     print(f"Cant't read masks: {e}")
     print("Make sure you have run compute_histograms.py, generate_xx_probabilities.py, segment_from_distributions,\n"+
           "and segment-blood-cc.py")
     sys.exit(-1)
 
-    
+
 P0_binfile            = f"{binary_root}/segmented/P0/{segment_scale}x/{sample}.uint16"
 P1_binfile            = f"{binary_root}/segmented/P1/{segment_scale}x/{sample}.uint16"
 edt_binfile           = f"{binary_root}/fields/implant-edt/{mask_scale}x/{sample}.uint16"
@@ -92,4 +92,4 @@ def np_save(path,data):
                     d_min, d_max, theta_min, theta_max,
                     tuple(bbox.flatten()), tuple(Muvwp.flatten()),
                     images, counts)
-                    
+

pre-cleanup-src/config/paths.py

@@ -0,0 +1 @@
+threadripper00/paths.py

src/doitall.py → pre-cleanup-src/doitall.py

@@ -60,17 +60,21 @@
 
     11)
 	for b $(seq 0 $nblocks); do python3 histogram_processing/optimize_distributions_flat.py $sample bone_region$b edt 4 0; done
+	;;
+
+    12)
 	for b $(seq 0 $nblocks); do python3 histogram_processing/compute_probabilities_flat.py $sample bone_region$b edt 10 0; done
 	;;    
 
-    12)
+    13)
 	python3 scripts/segment-from-distributions.py $sample 0 0 bone_region optimized_distributions
 	;;
 
-    13) for m in 0 1; do python3 preprocess/rescale-cupy-bin.py $sample segmented/P$m ; done
+    14) 
+	for m in 0 1; do python3 preprocess/rescale-cupy-bin.py $sample segmented/P$m ; done
 	;;
 
-    14)
+    15)
 	python3 segmentation/segment-blod-cc.py $sample
 	;;
 

src/pybind_kernels/histograms.cc → pre-cleanup-src/pybind_kernels/cpu/histograms.cc

@@ -6,17 +6,11 @@
 #include <stdio.h>
 #include <omp.h>
 #include <chrono>
-#include <iostream>
-#include <fstream>
 #include <tqdm.h>
 using namespace std;
 namespace py = pybind11;
 
-typedef uint16_t voxel_type;
-//typedef float    field_type;
-typedef uint16_t field_type;
-typedef uint8_t mask_type;
-typedef float gauss_type;
+#include "datatypes.hh"
 
 template <typename voxel_type>
   using np_array = py::array_t<voxel_type, py::array::c_style | py::array::forcecast>;
@@ -1170,9 +1164,6 @@ void otsu(
 
 PYBIND11_MODULE(histograms, m) {
     m.doc() = "2D histogramming plugin"; // optional module docstring
-    m.def("load_slice", &load_slice);
-    m.def("append_slice", &append_slice);
-    m.def("write_slice", &write_slice);
     m.def("axis_histogram_seq_cpu",  &axis_histogram_seq_cpu);
     m.def("axis_histogram_par_cpu",  &axis_histogram_par_cpu);
     m.def("axis_histogram_par_gpu",  &axis_histogram_par_gpu);
@@ -1182,9 +1173,5 @@ PYBIND11_MODULE(histograms, m) {
     m.def("masked_minmax", &masked_minmax);
     m.def("float_minmax", &float_minmax);
     m.def("gauss_filter_par_cpu", &gauss_filter_par_cpu);
-    m.def("dilate_3d_sphere_cpu", &morphology_3d_sphere_cpu<std::bit_or<mask_type>, false>);
-    m.def("erode_3d_sphere_cpu", &morphology_3d_sphere_cpu<std::bit_and<mask_type>, true>);
-    m.def("dilate_3d_sphere_gpu", &morphology_3d_sphere_gpu<std::bit_or<mask_type>, false>);
-    m.def("erode_3d_sphere_gpu", &morphology_3d_sphere_gpu<std::bit_and<mask_type>, true>);
     m.def("otsu", &otsu);
 }

proc-steps-checklist.txt

@@ -0,0 +1,6 @@
+0100 - Haven't checked. Requires ERDA to be set up properly.
+0200 - Verified!
+0300 - Verified!
+0400 - Verified!
+0500 - Verified!
+0600 - Verified!

src/Makefile

@@ -1,14 +1,65 @@
-CXX=icpc
-CXXFLAGS += -std=c++17
-CXXFLAGS += $(shell pkg-config --cflags fmt)
-CXXFLAGS += $(shell pkg-config --cflags hdf5)
-LIBS += $(shell pkg-config --libs fmt)
-LIBS += $(shell pkg-config --libs hdf5)
+# Define constants and collections
+PYTHON = python3.10
+PYBIND_FLAGS = $(shell $(PYTHON) -m pybind11 --include)
+CXXFLAGS += $(subst -I,-isystem ,$(PYBIND_FLAGS)) # We don't care about warnings from the python headers
+PYBIND_SUFFIX = $(shell $(PYTHON)-config --extension-suffix)
 
-%.o: %.cc
-	$(CXX) $(CXXFLAGS) -c $< 
+# Detect OS for OS specific changes
+ifeq ($(shell uname -s), Darwin) # Mac OSX
+CXX = g++-12 # Use homebrew gcc, as system gcc is an alias for clang
+CXXFLAGS += -undefined dynamic_lookup # https://pybind11.readthedocs.io/en/stable/compiling.html#building-manually
+CLEANUP += $(TARGETS) $(foreach TARGET, $(TARGETS), $(TARGET).dSYM) # These are also generated on Mac
+endif
 
-h5-blockmap: h5-blockmap.o
-	echo $(LIBS)
-	$(CXX) $(CXXFLAGS) $< $(LIBS) -o $@
+CPP_FOLDER=lib/cpp
+#CXXFLAGS += -I../contrib/cpptqdm/ -Iinclude
+CXXFLAGS += -I$(CPP_FOLDER)/include -march=native -Wall -Wextra -Wfloat-equal -Wundef -Wshadow -Wuninitialized -Winit-self -shared -fPIC -g -std=c++17 -O3
+PLATFORMS=cpu_seq cpu
+cpu_seq_CXX=$(CXX)
+cpu_seq_FLAGS=-Wno-unknown-pragmas -Wno-comment -Wconversion #-Weffc++
+cpu_CXX=$(cpu_seq_CXX)
+cpu_FLAGS=$(cpu_seq_FLAGS) -fopenmp
 
+LIBS=io geometry morphology
+TARGETS = $(foreach PLATFORM, $(PLATFORMS), $(foreach LIB, $(LIBS), $(CPP_FOLDER)/$(PLATFORM)/$(LIB)$(PYBIND_SUFFIX)))
+CLEANUP = $(TARGETS) $(foreach PLATFORM, $(PLATFORMS), $(CPP_FOLDER)/$(PLATFORM)/__pycache__)
+
+# Detect if OpenACC can be used
+ifneq (, $(shell which nvc++))
+PLATFORMS += gpu
+gpu_CXX = nvc++
+gpu_FLAGS = -acc=gpu -tp=native #-Minfo=accel
+gpu_FLAGS += -Xcudafe --display_error_number # Getting the warning codes for later suppression
+gpu_FLAGS += --diag_suppress 1626 # Remove the annoying pybind warning that routine is both inline and noinline
+gpu_FLAGS += --diag_suppress 9 # Remove the warning about nested comments
+gpu_FLAGS += -Wnvlink,-w # Disable nvlink warnings
+else
+$(info OpenACC compiler nvc++ not found. Compiling without.)
+endif
+
+all: $(TARGETS) pip_install
+
+define GEN_RULE
+$(CPP_FOLDER)/$(PLATFORM)/$(LIB)$(PYBIND_SUFFIX): pybind/$(LIB)-pybind.cc $(CPP_FOLDER)/$(PLATFORM)/$(LIB).cc $(CPP_FOLDER)/include/*.hh
+	$($(PLATFORM)_CXX) $($(PLATFORM)_FLAGS) $(CXXFLAGS) -I$(CPP_FOLDER)/$(PLATFORM) $$< -o $(CPP_FOLDER)/$(PLATFORM)/$(LIB)$(PYBIND_SUFFIX)
+endef
+
+$(foreach PLATFORM, $(PLATFORMS), \
+	$(foreach LIB, $(LIBS), \
+		$(eval $(GEN_RULE)) \
+	) \
+)
+
+# TODO lightweight openacc profiling can be done with the environment variable NV_ACC_TIME=1 !!!
+
+pip_install:
+	$(PYTHON) -m pip install -r requirements.txt
+
+test: all
+	$(PYTHON) -m pytest -n auto test
+
+test_%: test/test_%.py all
+	$(PYTHON) -m pytest -n auto $<
+
+clean:
+	rm -rf $(CLEANUP) __pycache__ test/__pycache__ .pytest_cache lib/cpp/**/*.so

src/config/paths.py

@@ -0,0 +1,12 @@
+data_root = "/data"
+fast_root = "/data_fast"
+
+hdf5_root        = f"{data_root}/MAXIBONE/Goats/tomograms"
+hdf5_root_fast   = f"{fast_root}/MAXIBONE/Goats/tomograms"
+binary_root      = f"{hdf5_root}/binary"
+binary_root_fast = f"{hdf5_root_fast}/binary"
+
+esrf_data_local     = f"{hdf5_root}/ESRF/"
+esrf_data_sftp      =  "/XNS/XrayImaging/MiG/manjula.esci.nbi.dk.2_localhost/"
+esrf_implants_root  = f"{esrf_data_local}/esrf_dental_implants_april_2013/"
+esrf_granules_root  = f"{esrf_data_local}/esrf_dental_granules_july_2012/"

src/exploration/print_cpp_type_ids.cpp

@@ -0,0 +1,24 @@
+#include<iostream>
+#include<stdint.h>
+
+int main() {
+    /*
+    This class is used to print out the code of the type.
+    This is handy when debugging the templated type at runtime.
+    */
+
+    std::cout << "int8 " << typeid(int8_t).name() << std::endl;
+    std::cout << "int16 " << typeid(int16_t).name() << std::endl;
+    std::cout << "int32 " << typeid(int32_t).name() << std::endl;
+    std::cout << "int64 " << typeid(int64_t).name() << std::endl;
+    std::cout << "int128 " << typeid(__int128_t).name() << std::endl;
+
+    std::cout << "uint8 " << typeid(uint8_t).name() << std::endl;
+    std::cout << "uint16 " << typeid(uint16_t).name() << std::endl;
+    std::cout << "uint32 " << typeid(uint32_t).name() << std::endl;
+    std::cout << "uint64 " << typeid(uint64_t).name() << std::endl;
+    std::cout << "uint128 " << typeid(__uint128_t).name() << std::endl;
+
+    std::cout << "float " << typeid(float).name() << std::endl;
+    std::cout << "double " << typeid(double).name() << std::endl;
+}