I am wrapping some functions from the ITK-5.3 library in Julia, and i have the following C++ defined wrappers, specifically I have a function that takes an image and reorients the image into LPS orientation, (the filter works flawlessly on linux, macos ) but on windows, it is repeatedly making most of the pixels zero, and filling some with a garbage value.
Specifically I want help with the reorientToLPS() function
Following is the code : i am using :
#include <memory>
#include <type_traits>
// Check if std::make_unique_for_overwrite is available
#if !__cpp_lib_make_unique_for_overwrite
// Define the custom implementation in the global namespace
template<typename T>
typename std::enable_if<std::is_array<T>::value, std::unique_ptr<T>>::type
make_unique_for_overwrite(std::size_t size) {
return std::unique_ptr<T>(new typename std::remove_extent<T>::type[size]);
}
template<typename T, typename... Args>
typename std::enable_if<!std::is_array<T>::value, std::unique_ptr<T>>::type
make_unique_for_overwrite(Args&&... args) {
return std::unique_ptr<T>(new T(std::forward<Args>(args)...));
}
// Inject the custom implementation into the std namespace
namespace std {
using ::make_unique_for_overwrite;
}
#endif // !__cpp_lib_make_unique_for_overwrite
#include "jlcxx/jlcxx.hpp"
#include "jlcxx/array.hpp"
#include "jlcxx/tuple.hpp"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkImageSeriesReader.h"
#include "itkImageSeriesWriter.h"
#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkOrientImageFilter.h"
#include "itkMetaDataObject.h"
#include "itkSpatialOrientation.h"
#include "itkSpatialOrientationAdapter.h"
#include "itkImageIOBase.h"
#include <vector>
#include <stdexcept>
#include <iostream>
#include <iomanip>
#include <sstream>
#include "itkImage.h"
using ImageType = itk::Image<float, 3>;
class ITKImageWrapper {
private:
ImageType::Pointer m_Image;
public:
// Constructor to load a NIfTI file or DICOM series
ITKImageWrapper(const std::string& path, bool isDicom = false) {
if (isDicom) {
// Handle DICOM directory
auto namesGenerator = itk::GDCMSeriesFileNames::New();
namesGenerator->SetUseSeriesDetails(true);
namesGenerator->SetDirectory(path);
const auto& seriesUIDs = namesGenerator->GetSeriesUIDs();
if (seriesUIDs.empty()) {
throw std::runtime_error("No DICOM series found in the directory.");
}
// Read the first series
const std::string seriesIdentifier = seriesUIDs.begin()->c_str();
const auto& fileNames = namesGenerator->GetFileNames(seriesIdentifier);
// Read the DICOM series
auto reader = itk::ImageSeriesReader<ImageType>::New();
auto dicomIO = itk::GDCMImageIO::New();
reader->SetImageIO(dicomIO);
reader->SetFileNames(fileNames);
reader->Update();
m_Image = reader->GetOutput();
} else {
// Handle NIfTI file
auto reader = itk::ImageFileReader<ImageType>::New();
reader->SetFileName(path);
reader->Update();
m_Image = reader->GetOutput();
}
}
// Get origin
std::vector<double> getOrigin() const {
auto origin = m_Image->GetOrigin();
return {origin[0], origin[1], origin[2]};
}
// Get spacing
std::vector<double> getSpacing() const {
auto spacing = m_Image->GetSpacing();
return {spacing[0], spacing[1], spacing[2]};
}
// Get size
std::vector<std::size_t> getSize() const {
auto size = m_Image->GetLargestPossibleRegion().GetSize();
return {size[0], size[1], size[2]};
}
// Get pixel data
std::vector<float> getPixelData() const {
std::vector<float> data;
itk::ImageRegionConstIterator<ImageType> it(m_Image, m_Image->GetLargestPossibleRegion());
for (it.GoToBegin(); !it.IsAtEnd(); ++it) {
data.push_back(it.Get());
}
return data;
}
// Get direction
std::vector<double> getDirection() const {
auto direction = m_Image->GetDirection();
std::vector<double> direction_flat;
for (unsigned int i = 0; i < 3; ++i) {
for (unsigned int j = 0; j < 3; ++j) {
direction_flat.push_back(direction[i][j]);
}
}
return direction_flat;
}
void writeImage(const std::string& filename, bool isDicom = false) const {
if (isDicom) {
using OutputPixelType = signed short;
using OutputImageType = itk::Image<OutputPixelType, 2>;
using SeriesWriterType = itk::ImageSeriesWriter<ImageType, OutputImageType>;
auto writer = SeriesWriterType::New();
auto dicomIO = itk::GDCMImageIO::New();
// Configure DICOM settings
dicomIO->SetComponentType(itk::IOComponentEnum::SHORT);
writer->SetImageIO(dicomIO);
// Get image properties
ImageType::RegionType region = m_Image->GetLargestPossibleRegion();
ImageType::SizeType size = region.GetSize();
unsigned int numberOfSlices = size[2];
// Generate filenames
std::vector<std::string> outputFileNames;
for (unsigned int i = 0; i < numberOfSlices; ++i) {
std::ostringstream ss;
ss << filename << "/slice_" << std::setw(3) << std::setfill('0') << i << ".dcm";
outputFileNames.push_back(ss.str());
}
try {
writer->SetInput(m_Image);
writer->SetFileNames(outputFileNames);
writer->Update();
} catch (const itk::ExceptionObject& e) {
std::cerr << "Error writing DICOM series: " << e.what() << std::endl;
throw;
}
} else {
// NIfTI writing remains unchanged
auto writer = itk::ImageFileWriter<ImageType>::New();
writer->SetFileName(filename);
writer->SetInput(m_Image);
writer->Update();
}
}
void reorientToLPS() {
try {
// Check if image is valid
if (!m_Image) {
throw std::runtime_error("No image loaded");
}
// Debug: Print stats before
std::cout << "DEBUG: Before reorientation - checking pixel data..." << std::endl;
itk::ImageRegionConstIterator<ImageType> it_before(m_Image, m_Image->GetLargestPossibleRegion());
size_t count_before = 0;
size_t zeros_before = 0;
double sum_before = 0.0;
for (it_before.GoToBegin(); !it_before.IsAtEnd() && count_before < 1000; ++it_before, ++count_before) {
float val = it_before.Get();
if (val == 0.0f) zeros_before++;
sum_before += val;
}
std::cout << "DEBUG: First 1000 voxels - zeros: " << zeros_before << ", mean: " << sum_before/count_before << std::endl;
// Check if already in LPS orientation
auto currentDirection = m_Image->GetDirection();
std::cout << "Current direction matrix:\n" << currentDirection << std::endl;
// Check if already LPS (identity matrix means LPS)
bool isAlreadyLPS = true;
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) {
double expected = (i == j) ? 1.0 : 0.0;
if (std::abs(currentDirection[i][j] - expected) > 1e-6) {
isAlreadyLPS = false;
break;
}
}
if (!isAlreadyLPS) break;
}
if (isAlreadyLPS) {
std::cout << "Image is already in LPS orientation, skipping reorientation" << std::endl;
return;
}
using OrientFilterType = itk::OrientImageFilter<ImageType, ImageType>;
auto orientFilter = OrientFilterType::New();
orientFilter->UseImageDirectionOn();
orientFilter->SetInput(m_Image);
orientFilter->SetDesiredCoordinateOrientation(
itk::SpatialOrientationEnums::ValidCoordinateOrientations::ITK_COORDINATE_ORIENTATION_LPS);
std::cout << "Filter configured, updating..." << std::endl;
orientFilter->Update();
ImageType::Pointer tempOutput = orientFilter->GetOutput();
// Debug: Check output before disconnect
std::cout << "DEBUG: After filter, before disconnect..." << std::endl;
itk::ImageRegionConstIterator<ImageType> it_after(tempOutput, tempOutput->GetLargestPossibleRegion());
size_t count_after = 0;
size_t zeros_after = 0;
double sum_after = 0.0;
for (it_after.GoToBegin(); !it_after.IsAtEnd() && count_after < 1000; ++it_after, ++count_after) {
float val = it_after.Get();
if (val == 0.0f) zeros_after++;
sum_after += val;
}
std::cout << "DEBUG: First 1000 voxels after filter - zeros: " << zeros_after << ", mean: " << sum_after/count_after << std::endl;
// Check if filter corrupted the data
if (zeros_after > zeros_before + 100 || std::abs(sum_after) > 1e30) {
std::cerr << "WARNING: Filter corrupted data on Windows, skipping reorientation" << std::endl;
std::cerr << "This is a known issue with ITK OrientImageFilter on Windows" << std::endl;
return;
}
// If we get here, the filter worked correctly
tempOutput->DisconnectPipeline();
m_Image = tempOutput;
std::cout << "Reorientation completed successfully" << std::endl;
std::cout << "New direction matrix:\n" << m_Image->GetDirection() << std::endl;
} catch(const itk::ExceptionObject& e) {
std::cerr << "ITK Error during reorientation: " << e.what() << std::endl;
throw;
} catch(const std::exception& e) {
std::cerr << "Standard Error during reorientation: " << e.what() << std::endl;
throw;
}
}
};
JLCXX_MODULE define_julia_module(jlcxx::Module& mod)
{
mod.add_type<ITKImageWrapper>("ITKImageWrapper")
.constructor<std::string>() // For NIfTI files
.constructor<std::string, bool>() // For DICOM series
.method("getOrigin", &ITKImageWrapper::getOrigin)
.method("getSpacing", &ITKImageWrapper::getSpacing)
.method("getSize", &ITKImageWrapper::getSize)
.method("getPixelData", &ITKImageWrapper::getPixelData)
.method("getDirection", &ITKImageWrapper::getDirection)
.method("writeImage", &ITKImageWrapper::writeImage)
.method("reorientToLPS", &ITKImageWrapper::reorientToLPS);
}
I had a simpler version of the reorientToLPS function
void reorientToLPS() {
try {
using OrientFilterType = itk::OrientImageFilter<ImageType, ImageType>;
auto orientFilter = OrientFilterType::New();
orientFilter->UseImageDirectionOn();
orientFilter->SetInput(m_Image);
orientFilter->SetDesiredCoordinateOrientation(itk::SpatialOrientationEnums::ValidCoordinateOrientations::ITK_COORDINATE_ORIENTATION_LPS);
orientFilter->Update();
ImageType::Pointer tempOutput = orientFilter->GetOutput();
tempOutput->DisconnectPipeline();
m_Image = tempOutput;
}
catch(const itk::ExceptionObject& e) {
std::cerr << "Error during reorientation: " << e.what() << std::endl;
throw;
}
}
in first case , i simply skipped orientation, if the orientation was LPS or the data after orientation got corrupt
in the second version (simpler one) which i was using originally, the reorientToLPS on windows simply zeroes out all data for my image.
Following is my CMAKE file :
cmake_policy(SET CMP0025 NEW)
project(ITKIOWrapper)
cmake_minimum_required(VERSION 3.14)
set(CMAKE_CXX_STANDARD 20)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
if(APPLE)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -stdlib=libc++")
set(CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS} -stdlib=libc++")
endif()
if(WIN32)
add_definitions(-D_USE_MATH_DEFINES)
add_definitions(-DNOMINMAX)
endif()
set(CMAKE_LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DITK_LEGACY_REMOVE=OFF")
set(ITK_BUILD_TESTING OFF)
set(BUILD_TESTING OFF)
find_package(JlCxx)
include_directories(${JlCxx_INCLUDE_DIRS})
get_target_property(JlCxx_location JlCxx::cxxwrap_julia LOCATION)
get_filename_component(JlCxx_location ${JlCxx_location} DIRECTORY)
set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib;${JlCxx_location}")
message(STATUS "Found JlCxx at ${JlCxx_location}")
# Update the ITK components list
set(ITK_COMPONENTS
ITKCommon
ITKIOImageBase
ITKIOGDCM
ITKIONIFTI
ITKIONRRD
ITKIOHDF5
ITKIOMeta
ITKIOPNG
ITKIOTIFF
ITKImageIO
)
find_package(HDF5 REQUIRED COMPONENTS CXX C HL)
find_package(ITK REQUIRED COMPONENTS ${ITK_COMPONENTS})
include(${ITK_USE_FILE})
message(STATUS "Found ITK at ${ITK_USE_FILE}")
add_library(ITKIOWrapper SHARED itklib.cpp)
target_link_libraries(ITKIOWrapper
PUBLIC
JlCxx::cxxwrap_julia
${ITK_LIBRARIES}
)
target_include_directories(ITKIOWrapper
PUBLIC
${JlCxx_INCLUDE_DIRS}
${ITK_INCLUDE_DIRS}
)
# Set properties for macOS
if(APPLE)
set_target_properties(ITKIOWrapper PROPERTIES
INSTALL_RPATH "@loader_path"
BUILD_WITH_INSTALL_RPATH TRUE
INSTALL_NAME_DIR "@rpath"
)
endif()
install(TARGETS ITKIOWrapper
LIBRARY DESTINATION lib
RUNTIME DESTINATION lib
)
I am wondering if this is a bug with ITK-v5.3