From 01a9d09739f773df222eb36639c5617980d91ccd Mon Sep 17 00:00:00 2001 From: Simon Date: Sat, 2 Nov 2019 16:43:36 +0100 Subject: [PATCH 01/37] Transform IO --- src/deform_lib/CMakeLists.txt | 2 + src/deform_lib/registration/affine.cpp | 129 +++++++++++++++++++++++++ src/deform_lib/registration/affine.h | 36 +++++++ test/CMakeLists.txt | 3 +- test/test_affine.cpp | 107 ++++++++++++++++++++ 5 files changed, 276 insertions(+), 1 deletion(-) create mode 100644 src/deform_lib/registration/affine.cpp create mode 100644 src/deform_lib/registration/affine.h create mode 100644 test/test_affine.cpp diff --git a/src/deform_lib/CMakeLists.txt b/src/deform_lib/CMakeLists.txt index c8e1270..64ca77a 100644 --- a/src/deform_lib/CMakeLists.txt +++ b/src/deform_lib/CMakeLists.txt @@ -16,6 +16,8 @@ set(DEFORM_LIB_SRCS "registration/cost_functions/sub_function.h" "registration/cost_functions/unary_function.h" + "registration/affine.cpp" + "registration/affine.h" "registration/block_change_flags.cpp" "registration/block_change_flags.h" "registration/blocked_graph_cut_optimizer.h" diff --git a/src/deform_lib/registration/affine.cpp b/src/deform_lib/registration/affine.cpp new file mode 100644 index 0000000..7ed5142 --- /dev/null +++ b/src/deform_lib/registration/affine.cpp @@ -0,0 +1,129 @@ +#include "affine.h" + +#include +#include +#include + +#include + +namespace { + std::string trim_string(const std::string& str) + { + size_t i = 0, j = str.size()-1; + + while (i < j && isspace(str[i])) i++; + while (j > i && isspace(str[j])) j--; + + return std::string(str, i, j+1); + } +} + +AffineTransform::AffineTransform() : + _matrix(Matrix3x3f::Identity), + _offset{0, 0, 0} +{ +} +AffineTransform::AffineTransform(const Matrix3x3f& matrix, const float3& offset) : + _matrix(matrix), + _offset(offset) +{ +} +AffineTransform::~AffineTransform() +{ +} + +const Matrix3x3f& AffineTransform::matrix() const +{ + return _matrix; +} +const float3& AffineTransform::offset() const +{ + return _offset; +} + + +AffineTransform parse_affine_transform_file( + const std::string& filename +) +{ + // TODO: Cleanup, should probably be less strict (ordering of keys, etc) + + std::string line; + + std::ifstream f {filename, std::ifstream::in}; + if (!f.is_open()) { + throw ValidationError("Failed to open file '" + filename + "'"); + } + + std::getline(f, line); + if (line != "#Insight Transform File V1.0") { + throw ValidationError("Unsupported transform file format"); + } + + // #Transform + std::getline(f, line); + + auto split = [](const std::string& line) { + size_t p = line.find(':'); + if (p != std::string::npos) { + return std::pair( + trim_string(line.substr(0, p)), + trim_string(line.substr(p+1, std::string::npos)) + ); + } + return std::pair( + "","" + ); + }; + + // Transform: AffineTransform_double_3_3 + std::getline(f, line); + auto p = split(line); + + if (p.first != "Transform") { + throw ValidationError("Invalid file format"); + } + + if (p.second != "AffineTransform_double_3_3") { + throw ValidationError("Invalid transform name, only 'AffineTransform_double_3_3' is supported."); + } + + // Parameters: *12 + std::getline(f, line); + p = split(line); + if (p.first != "Parameters") { + throw ValidationError("Invalid file format"); + } + + Matrix3x3f matrix = Matrix3x3f::Identity; + float* pmatrix = &matrix._rows[0].x; + std::istringstream s_p(p.second); + + for (int i = 0; i < 9; ++i) { + s_p >> pmatrix[i]; + } + + float3 translation; + s_p >> translation.x; + s_p >> translation.y; + s_p >> translation.z; + + // FixedParameters: + std::getline(f, line); + p = split(line); + if (p.first != "FixedParameters") { + throw ValidationError("Invalid file format"); + } + + std::istringstream s_fp(p.second); + + float3 fixed_parameters; + s_fp >> fixed_parameters.x; + s_fp >> fixed_parameters.y; + s_fp >> fixed_parameters.z; + + float3 center = fixed_parameters; + float3 offset = translation + center - matrix * center; + + return AffineTransform(matrix, offset); +} diff --git a/src/deform_lib/registration/affine.h b/src/deform_lib/registration/affine.h new file mode 100644 index 0000000..2798199 --- /dev/null +++ b/src/deform_lib/registration/affine.h @@ -0,0 +1,36 @@ +#pragma once + +#include + +#include +#include + +class AffineTransform +{ +public: + AffineTransform(); + AffineTransform(const Matrix3x3f& matrix, const float3& offset); + ~AffineTransform(); + + inline float3 transform_point(const float3& pt) + { + return _matrix * pt + _offset; + } + + const Matrix3x3f& matrix() const; + const float3& offset() const; + +private: + Matrix3x3f _matrix; + float3 _offset; +}; + +/*! + * \brief Parses an affine transform file. + * + * Supports: + * * Insight Transform File V1.0 + */ +AffineTransform parse_affine_transform_file( + const std::string& filename +); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c6afc82..ef6b09f 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,4 +1,5 @@ -set(DEFORM_TEST_SRCS +set(DEFORM_TEST_SRCS + "test_affine.cpp" "test_arg_parser.cpp" "test_gpu.cpp" "test_gpu_cost.cpp" diff --git a/test/test_affine.cpp b/test/test_affine.cpp new file mode 100644 index 0000000..e648c74 --- /dev/null +++ b/test/test_affine.cpp @@ -0,0 +1,107 @@ +#include "catch.hpp" + +#include +#include + +#include +#include +#include +#include + +using namespace Catch; + +TEST_CASE("parse_affine") +{ + SECTION("test_parsing") + { + const char* affine_str = + "#Insight Transform File V1.0\n" + "#Transform 0\n" + "Transform: AffineTransform_double_3_3\n" + "Parameters: 1 2 3 4 5 6 7 8 9 10 11 12\n" + "FixedParameters: 13 14 15\n"; + + const std::string filename {"affine.txt"}; + std::ofstream f(filename); + assert(f.is_open()); + + f.write(affine_str, std::strlen(affine_str)); + f.close(); + + AffineTransform transform = parse_affine_transform_file(filename); + + float3 pt1 = transform.transform_point(float3{1,2,3}); + float3 pt2 = transform.transform_point(float3{-5.1,-7.2,-8.3}); + + // Expected values computed by SimpleITK + + CHECK(-49.0f == pt1.x); + CHECK(-155.0f == pt1.y); + CHECK(-261.0f == pt1.z); + + CHECK(-107.4f == pt2.x); + CHECK(-293.2f == pt2.y); + CHECK(-479.0f == pt2.z); + } + + SECTION("test_unexisting_file") + { + const std::string filename {"non_existing_file.pts"}; + + CHECK_THROWS_WITH(parse_affine_transform_file(filename), + Contains("Failed to open")); + } + + SECTION("test_invalid_file_1") + { + const char* affine_str = + "#Transform 0\n" + "Transform: AffineTransform_double_3_3\n" + "Parameters: 1 2 3 4 5 6 7 8 9 10 11 12\n" + "FixedParameters: 13 14 15\n"; + + const std::string filename {"affine.txt"}; + std::ofstream f(filename); + assert(f.is_open()); + + f.write(affine_str, std::strlen(affine_str)); + f.close(); + + CHECK_THROWS_WITH(parse_affine_transform_file(filename), + Contains("Unsupported transform file format")); + } + + SECTION("test_invalid_file_2") + { + const char* affine_str = + "#Insight Transform File V1.0\n"; + + const std::string filename {"affine.txt"}; + std::ofstream f(filename); + assert(f.is_open()); + + f.write(affine_str, std::strlen(affine_str)); + f.close(); + + CHECK_THROWS_WITH(parse_affine_transform_file(filename), + Contains("Invalid file format")); + } + + SECTION("test_invalid_file_2") + { + const char* affine_str = + "#Insight Transform File V1.0\n" + "#Transform 0\n" + "Transform: AffineTransform_double_1_1\n"; + + const std::string filename {"affine.txt"}; + std::ofstream f(filename); + assert(f.is_open()); + + f.write(affine_str, std::strlen(affine_str)); + f.close(); + + CHECK_THROWS_WITH(parse_affine_transform_file(filename), + Contains("Invalid transform name")); + } +} From 5db23758f0c53cef9ad6a4e00b220f3adbf958be Mon Sep 17 00:00:00 2001 From: Simon Date: Sat, 2 Nov 2019 17:06:48 +0100 Subject: [PATCH 02/37] Improvements to parsing --- src/deform_lib/registration/affine.cpp | 85 ++++++++++++-------------- 1 file changed, 40 insertions(+), 45 deletions(-) diff --git a/src/deform_lib/registration/affine.cpp b/src/deform_lib/registration/affine.cpp index 7ed5142..f844438 100644 --- a/src/deform_lib/registration/affine.cpp +++ b/src/deform_lib/registration/affine.cpp @@ -46,8 +46,6 @@ AffineTransform parse_affine_transform_file( const std::string& filename ) { - // TODO: Cleanup, should probably be less strict (ordering of keys, etc) - std::string line; std::ifstream f {filename, std::ifstream::in}; @@ -63,64 +61,61 @@ AffineTransform parse_affine_transform_file( // #Transform std::getline(f, line); - auto split = [](const std::string& line) { + std::map keyvalues; + while (std::getline(f, line)) + { size_t p = line.find(':'); if (p != std::string::npos) { - return std::pair( - trim_string(line.substr(0, p)), - trim_string(line.substr(p+1, std::string::npos)) - ); - } - return std::pair( - "","" - ); - }; + std::string key = trim_string(line.substr(0, p)); + std::string value = trim_string(line.substr(p+1, std::string::npos)); - // Transform: AffineTransform_double_3_3 - std::getline(f, line); - auto p = split(line); - - if (p.first != "Transform") { - throw ValidationError("Invalid file format"); + keyvalues[key] = value; + } } - if (p.second != "AffineTransform_double_3_3") { - throw ValidationError("Invalid transform name, only 'AffineTransform_double_3_3' is supported."); + if (keyvalues.find("Transform") != keyvalues.end()) { + if (keyvalues["Transform"] != "AffineTransform_double_3_3") { + throw ValidationError( + "Invalid transform name, only 'AffineTransform_double_3_3' is supported." + ); + } } - - // Parameters: *12 - std::getline(f, line); - p = split(line); - if (p.first != "Parameters") { - throw ValidationError("Invalid file format"); + else { + throw ValidationError("Invalid file format: Missing \"Transform\""); } Matrix3x3f matrix = Matrix3x3f::Identity; - float* pmatrix = &matrix._rows[0].x; - std::istringstream s_p(p.second); + float3 translation{0, 0, 0}; - for (int i = 0; i < 9; ++i) { - s_p >> pmatrix[i]; - } + if (keyvalues.find("Parameters") != keyvalues.end()) { + // Parameters: *12 - float3 translation; - s_p >> translation.x; - s_p >> translation.y; - s_p >> translation.z; + std::istringstream ss(keyvalues["Parameters"]); + + float* pmatrix = &matrix._rows[0].x; + for (int i = 0; i < 9; ++i) { + ss >> pmatrix[i]; + } - // FixedParameters: - std::getline(f, line); - p = split(line); - if (p.first != "FixedParameters") { - throw ValidationError("Invalid file format"); + ss >> translation.x; + ss >> translation.y; + ss >> translation.z; + } + else { + throw ValidationError("Invalid file format: Missing \"Parameters\""); } - std::istringstream s_fp(p.second); + float3 fixed_parameters{0, 0, 0}; - float3 fixed_parameters; - s_fp >> fixed_parameters.x; - s_fp >> fixed_parameters.y; - s_fp >> fixed_parameters.z; + // FixedParameters are optional + if (keyvalues.find("FixedParameters") != keyvalues.end()) { + // FixedParameters: + + std::istringstream ss(keyvalues["FixedParameters"]); + ss >> fixed_parameters.x; + ss >> fixed_parameters.y; + ss >> fixed_parameters.z; + } float3 center = fixed_parameters; float3 offset = translation + center - matrix * center; From 958eed5209ce646c585bb796e01f4441240188da Mon Sep 17 00:00:00 2001 From: Simon Date: Sat, 2 Nov 2019 17:41:41 +0100 Subject: [PATCH 03/37] transform_point tests --- test/test_affine.cpp | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/test/test_affine.cpp b/test/test_affine.cpp index e648c74..05d9f4d 100644 --- a/test/test_affine.cpp +++ b/test/test_affine.cpp @@ -104,4 +104,44 @@ TEST_CASE("parse_affine") CHECK_THROWS_WITH(parse_affine_transform_file(filename), Contains("Invalid transform name")); } + + SECTION("test_transform_point") + { + // Identity transform + AffineTransform transform1( + Matrix3x3f::Identity, + float3{0, 0, 0} + ); + + float3 pt1 = transform1.transform_point(float3{1, 2, 3}); + CHECK(1.0f == pt1.x); + CHECK(2.0f == pt1.y); + CHECK(3.0f == pt1.z); + + // Scaling + AffineTransform transform2( + Matrix3x3f{ + float3{2, 0, 0}, + float3{0, 3, 0}, + float3{0, 0, 4} + }, + float3{0, 0, 0} + ); + + float3 pt2 = transform2.transform_point(float3{1, 2, 3}); + CHECK(2.0f == pt2.x); + CHECK(6.0f == pt2.y); + CHECK(12.0f == pt2.z); + + // Translation + AffineTransform transform3( + Matrix3x3f::Identity, + float3{4, 6, 8} + ); + + float3 pt3 = transform3.transform_point(float3{1, 2, 3}); + CHECK(5.0f == pt3.x); + CHECK(8.0f == pt3.y); + CHECK(11.0f == pt3.z); + } } From b50d139f5b368bd07b5a188a5ed27f21619ddbb3 Mon Sep 17 00:00:00 2001 From: Simon Date: Sat, 2 Nov 2019 18:11:39 +0100 Subject: [PATCH 04/37] Extract displacement vector from affine transformation --- src/deform_lib/registration/affine.h | 6 +++++ test/test_affine.cpp | 40 ++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+) diff --git a/src/deform_lib/registration/affine.h b/src/deform_lib/registration/affine.h index 2798199..a403db6 100644 --- a/src/deform_lib/registration/affine.h +++ b/src/deform_lib/registration/affine.h @@ -17,6 +17,12 @@ class AffineTransform return _matrix * pt + _offset; } + // Returns the resulting displacement vector at the given point + inline float3 displacement_vector(const float3& pt) + { + return transform_point(pt) - pt; + } + const Matrix3x3f& matrix() const; const float3& offset() const; diff --git a/test/test_affine.cpp b/test/test_affine.cpp index 05d9f4d..4a47be9 100644 --- a/test/test_affine.cpp +++ b/test/test_affine.cpp @@ -144,4 +144,44 @@ TEST_CASE("parse_affine") CHECK(8.0f == pt3.y); CHECK(11.0f == pt3.z); } + SECTION("test_displacement_vector") + { + // Identity transform + AffineTransform transform1( + Matrix3x3f::Identity, + float3{0, 0, 0} + ); + + float3 d1 = transform1.displacement_vector(float3{2,3,4}); + CHECK(0.0f == d1.x); + CHECK(0.0f == d1.y); + CHECK(0.0f == d1.z); + + // Scaling + AffineTransform transform2( + Matrix3x3f{ + float3{2, 0, 0}, + float3{0, 3, 0}, + float3{0, 0, 4} + }, + float3{0, 0, 0} + ); + + float3 pt2 = transform2.displacement_vector(float3{1, 2, 3}); + CHECK(1.0f == pt2.x); + CHECK(4.0f == pt2.y); + CHECK(9.0f == pt2.z); + + // Translation + AffineTransform transform3( + Matrix3x3f::Identity, + float3{4, 6, 8} + ); + + float3 pt3 = transform3.displacement_vector(float3{1, 2, 3}); + CHECK(4.0f == pt3.x); + CHECK(6.0f == pt3.y); + CHECK(8.0f == pt3.z); + + } } From d1fd0fe795ff53e413fd7b3e61e883739aaa0fc7 Mon Sep 17 00:00:00 2001 From: Simon Date: Sat, 2 Nov 2019 18:20:01 +0100 Subject: [PATCH 05/37] Fix --- test/test_affine.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_affine.cpp b/test/test_affine.cpp index 4a47be9..0f4c9dd 100644 --- a/test/test_affine.cpp +++ b/test/test_affine.cpp @@ -31,7 +31,7 @@ TEST_CASE("parse_affine") AffineTransform transform = parse_affine_transform_file(filename); float3 pt1 = transform.transform_point(float3{1,2,3}); - float3 pt2 = transform.transform_point(float3{-5.1,-7.2,-8.3}); + float3 pt2 = transform.transform_point(float3{-5.1f,-7.2f,-8.3f}); // Expected values computed by SimpleITK From 802234bf0a864c0fb4ee008f458da538c7749008 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 4 Nov 2019 09:38:48 +0100 Subject: [PATCH 06/37] Adding affine transform to displacement field --- src/deform_lib/registration/affine.h | 4 ++-- src/deform_lib/registration/displacement_field.h | 8 +++++++- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/deform_lib/registration/affine.h b/src/deform_lib/registration/affine.h index a403db6..8e586d4 100644 --- a/src/deform_lib/registration/affine.h +++ b/src/deform_lib/registration/affine.h @@ -12,13 +12,13 @@ class AffineTransform AffineTransform(const Matrix3x3f& matrix, const float3& offset); ~AffineTransform(); - inline float3 transform_point(const float3& pt) + inline float3 transform_point(const float3& pt) const { return _matrix * pt + _offset; } // Returns the resulting displacement vector at the given point - inline float3 displacement_vector(const float3& pt) + inline float3 displacement(const float3& pt) const { return transform_point(pt) - pt; } diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index 0b97adb..897a3f2 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -3,6 +3,7 @@ #include #include +#include "affine.h" #include "settings.h" /** Wrapper for a displacement field image. @@ -56,7 +57,10 @@ class DisplacementField // Returns coordinates in world space inline float3 transform_index(const int3& p) const { - return _df.index2point(p) + get(p); + + void set_affine_transform(const AffineTransform& transform) + { + _affine = transform; } dim3 size() const @@ -101,4 +105,6 @@ class DisplacementField Settings::UpdateRule _update_rule; stk::VolumeFloat3 _df; + AffineTransform _affine; + }; \ No newline at end of file From ca1a215a803c50d39356531770d1cc08bc69a0f6 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 4 Nov 2019 10:11:51 +0100 Subject: [PATCH 07/37] Refactoring, utilize a update field for updates --- .../blocked_graph_cut_optimizer.h | 7 ++- .../blocked_graph_cut_optimizer.inl | 27 +++------ .../registration/displacement_field.h | 57 +++++++++---------- .../registration/registration_engine.cpp | 26 +++++++-- 4 files changed, 62 insertions(+), 55 deletions(-) diff --git a/src/deform_lib/registration/blocked_graph_cut_optimizer.h b/src/deform_lib/registration/blocked_graph_cut_optimizer.h index 116b576..d737cc5 100644 --- a/src/deform_lib/registration/blocked_graph_cut_optimizer.h +++ b/src/deform_lib/registration/blocked_graph_cut_optimizer.h @@ -28,7 +28,8 @@ class BlockedGraphCutOptimizer UnaryFunction& unary_fn, Regularizer& binary_fn, float3 step_size, - DisplacementField& df + DisplacementField& df, + Settings::UpdateRule update_rule ); private: @@ -39,8 +40,8 @@ class BlockedGraphCutOptimizer const int3& block_dims, const int3& block_offset, const float3& delta, // delta in [mm] - DisplacementField& df, - DisplacementField& df_tmp + const DisplacementField& df, + DisplacementField& update_field ); double calculate_energy( diff --git a/src/deform_lib/registration/blocked_graph_cut_optimizer.inl b/src/deform_lib/registration/blocked_graph_cut_optimizer.inl index 37c3e15..492fd8e 100644 --- a/src/deform_lib/registration/blocked_graph_cut_optimizer.inl +++ b/src/deform_lib/registration/blocked_graph_cut_optimizer.inl @@ -32,21 +32,13 @@ void BlockedGraphCutOptimizer UnaryFunction& unary_fn, Regularizer& binary_fn, float3 step_size, - DisplacementField& df + DisplacementField& df, + Settings::UpdateRule update_rule ) { dim3 dims = df.size(); - // Create a buffer displacement field for when doing composition - DisplacementField df_buffer; - if (df.update_rule() == Settings::UpdateRule_Compositive) { - df_buffer = df.clone(); - } - else { - // Only need cloning for compositive updates, since Additive updates do independent - // updates. - df_buffer = df; - } + DisplacementField update_field(df.size()); // Setting the block size to (0, 0, 0) will disable blocking and run the whole volume int3 block_dims = _block_size; @@ -156,7 +148,7 @@ void BlockedGraphCutOptimizer block_offset, delta, df, - df_buffer + update_field ); if (block_changed) @@ -165,9 +157,7 @@ void BlockedGraphCutOptimizer change_flags.set_block(block_p, block_changed, use_shift == 1); } - if (df.update_rule() == Settings::UpdateRule_Compositive) { - df.copy_from(df_buffer); - } + df.update(update_field, (update_rule == Settings::UpdateRule_Compositive)); } } } @@ -198,8 +188,8 @@ bool BlockedGraphCutOptimizer const int3& block_dims, const int3& block_offset, const float3& delta, // delta in mm - DisplacementField& df, - DisplacementField& df_tmp + const DisplacementField& df, + DisplacementField& update_field ) { dim3 dims = df.size(); @@ -338,8 +328,7 @@ bool BlockedGraphCutOptimizer } if (solver.get_var(sub_x, sub_y, sub_z) == 1) { - int3 p{gx, gy, gz}; - df_tmp.set(p, df.get(p, delta)); + update_field.set(int3{gx, gy, gz}, delta); changed_flag = true; } } diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index 897a3f2..6293892 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -14,15 +14,15 @@ class DisplacementField { public: - DisplacementField() : _update_rule(Settings::UpdateRule_Additive) {} - DisplacementField( - const stk::VolumeFloat3& df, - Settings::UpdateRule update_rule = Settings::UpdateRule_Additive - ) : - _update_rule(update_rule), + DisplacementField() {} + DisplacementField(const stk::VolumeFloat3& df) : _df(df) { } + DisplacementField(const dim3& dims) : + _df(dims, float3{0, 0, 0}) + { + } ~DisplacementField() {} inline float3 get(const int3& p) const @@ -33,19 +33,19 @@ class DisplacementField // delta : Delta in world space (mm) inline float3 get(const int3& p, const float3& delta) const { - if (_update_rule == Settings::UpdateRule_Compositive) { - // Convert delta to fixed image space - // TODO: What about orientation? - float3 fp { - p.x + delta.x / _df.spacing().x, - p.y + delta.y / _df.spacing().y, - p.z + delta.z / _df.spacing().z - }; - return _df.linear_at(fp, stk::Border_Replicate) + delta; - } - else /*(_update_rule == Settings::UpdateRule_Additive)*/ { + // if (_update_rule == Settings::UpdateRule_Compositive) { + // // Convert delta to fixed image space + // // TODO: What about orientation? + // float3 fp { + // p.x + delta.x / _df.spacing().x, + // p.y + delta.y / _df.spacing().y, + // p.z + delta.z / _df.spacing().z + // }; + // return _df.linear_at(fp, stk::Border_Replicate) + delta; + // } + // else /*(_update_rule == Settings::UpdateRule_Additive)*/ { return _df(p) + delta; - } + //} } inline void set(const int3& p, const float3& d) @@ -57,6 +57,15 @@ class DisplacementField // Returns coordinates in world space inline float3 transform_index(const int3& p) const { + return float3{float(p.x), float(p.y), float(p.z)} + _df(p); + } + + void update(const DisplacementField& update_field, bool composite) + { + for (int3 p : _df.size()) { + _df(p) += update_field.get(p); + } + } void set_affine_transform(const AffineTransform& transform) { @@ -68,11 +77,6 @@ class DisplacementField return _df.size(); } - DisplacementField clone() const - { - return DisplacementField(_df.clone(), _update_rule); - } - void copy_from(const DisplacementField& other) { _df.copy_from(other._df); @@ -96,15 +100,10 @@ class DisplacementField return _df.valid(); } - Settings::UpdateRule update_rule() const - { - return _update_rule; - } private: - Settings::UpdateRule _update_rule; stk::VolumeFloat3 _df; AffineTransform _affine; -}; \ No newline at end of file +}; diff --git a/src/deform_lib/registration/registration_engine.cpp b/src/deform_lib/registration/registration_engine.cpp index f500a19..ac8f829 100644 --- a/src/deform_lib/registration/registration_engine.cpp +++ b/src/deform_lib/registration/registration_engine.cpp @@ -618,7 +618,7 @@ stk::Volume RegistrationEngine::execute() #endif for (int l = _settings.num_pyramid_levels-1; l >= 0; --l) { - DisplacementField df(_deformation_pyramid.volume(l), _settings.update_rule); + DisplacementField df(_deformation_pyramid.volume(l)); std::vector neighborhood = determine_neighborhood(l); @@ -648,7 +648,13 @@ stk::Volume RegistrationEngine::execute() _settings.levels[l].block_energy_epsilon, _settings.levels[l].max_iteration_count ); - optimizer.execute(unary_fn, binary_fn, _settings.levels[l].step_size, df); + optimizer.execute( + unary_fn, + binary_fn, + _settings.levels[l].step_size, + df, + _settings.update_rule + ); } #if defined(DF_ENABLE_GCO) else if (_settings.solver == Settings::Solver_GCO) { @@ -658,7 +664,13 @@ stk::Volume RegistrationEngine::execute() _settings.levels[l].block_energy_epsilon, _settings.levels[l].max_iteration_count ); - optimizer.execute(unary_fn, binary_fn, _settings.levels[l].step_size, df); + optimizer.execute( + unary_fn, + binary_fn, + _settings.levels[l].step_size, + df, + _settings.update_rule + ); } #endif #if defined(DF_ENABLE_GRIDCUT) @@ -669,7 +681,13 @@ stk::Volume RegistrationEngine::execute() _settings.levels[l].block_energy_epsilon, _settings.levels[l].max_iteration_count ); - optimizer.execute(unary_fn, binary_fn, _settings.levels[l].step_size, df); + optimizer.execute( + unary_fn, + binary_fn, + _settings.levels[l].step_size, + df, + _settings.update_rule + ); } #endif From 4d053fc256ab28a00f2dfa168094362733a1dbec Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 4 Nov 2019 10:18:37 +0100 Subject: [PATCH 08/37] Reset update field --- src/deform_lib/registration/blocked_graph_cut_optimizer.inl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/deform_lib/registration/blocked_graph_cut_optimizer.inl b/src/deform_lib/registration/blocked_graph_cut_optimizer.inl index 492fd8e..aa4980b 100644 --- a/src/deform_lib/registration/blocked_graph_cut_optimizer.inl +++ b/src/deform_lib/registration/blocked_graph_cut_optimizer.inl @@ -158,6 +158,8 @@ void BlockedGraphCutOptimizer } df.update(update_field, (update_rule == Settings::UpdateRule_Compositive)); + // TODO: Reset + update_field = DisplacementField(update_field.size()); } } } From 1afa03c0be7309adc9269232e2ccd5a4e701ed1a Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 4 Nov 2019 10:25:02 +0100 Subject: [PATCH 09/37] Woops --- src/deform_lib/registration/displacement_field.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index 6293892..f7eb7d0 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -57,7 +57,7 @@ class DisplacementField // Returns coordinates in world space inline float3 transform_index(const int3& p) const { - return float3{float(p.x), float(p.y), float(p.z)} + _df(p); + return _df.index2point(p) + _df(p); } void update(const DisplacementField& update_field, bool composite) From 811b2d419e94fb01fe64c0abefad31b1b658d3e4 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 4 Nov 2019 10:36:24 +0100 Subject: [PATCH 10/37] Composition --- src/deform_lib/registration/displacement_field.h | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index f7eb7d0..3e17d74 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -63,7 +63,19 @@ class DisplacementField void update(const DisplacementField& update_field, bool composite) { for (int3 p : _df.size()) { - _df(p) += update_field.get(p); + float3 d = update_field.get(p); + + if (composite) { + float3 fp { + p.x + d.x / _df.spacing().x, + p.y + d.y / _df.spacing().y, + p.z + d.z / _df.spacing().z + }; + _df(p) = _df.linear_at(fp, stk::Border_Replicate) + d; + } + else { + _df(p) += d; + } } } From b962009e58b252a30f1ecc8c03754a3fd130ad31 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 4 Nov 2019 10:52:08 +0100 Subject: [PATCH 11/37] Composite --- src/deform_lib/registration/displacement_field.h | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index 3e17d74..0830c48 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -63,18 +63,15 @@ class DisplacementField void update(const DisplacementField& update_field, bool composite) { for (int3 p : _df.size()) { - float3 d = update_field.get(p); - if (composite) { - float3 fp { - p.x + d.x / _df.spacing().x, - p.y + d.y / _df.spacing().y, - p.z + d.z / _df.spacing().z - }; - _df(p) = _df.linear_at(fp, stk::Border_Replicate) + d; + float3 p1 = _df.index2point(p); + float3 p2 = p1 + update_field.get(p); + float3 p3 = p2 + _df.linear_at_point(p2, stk::Border_Replicate); + + _df(p) = p3 - p1; } else { - _df(p) += d; + _df(p) += update_field.get(p); } } } From 11d4833cb283c2eaba7019eb35ce2101147d4ff5 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 4 Nov 2019 11:22:12 +0100 Subject: [PATCH 12/37] Composition in cost function --- .../blocked_graph_cut_optimizer.h | 1 + .../blocked_graph_cut_optimizer.inl | 6 +++-- .../registration/displacement_field.h | 23 ++++++++----------- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/deform_lib/registration/blocked_graph_cut_optimizer.h b/src/deform_lib/registration/blocked_graph_cut_optimizer.h index d737cc5..8549aeb 100644 --- a/src/deform_lib/registration/blocked_graph_cut_optimizer.h +++ b/src/deform_lib/registration/blocked_graph_cut_optimizer.h @@ -41,6 +41,7 @@ class BlockedGraphCutOptimizer const int3& block_offset, const float3& delta, // delta in [mm] const DisplacementField& df, + Settings::UpdateRule update_rule, DisplacementField& update_field ); diff --git a/src/deform_lib/registration/blocked_graph_cut_optimizer.inl b/src/deform_lib/registration/blocked_graph_cut_optimizer.inl index aa4980b..8b9b213 100644 --- a/src/deform_lib/registration/blocked_graph_cut_optimizer.inl +++ b/src/deform_lib/registration/blocked_graph_cut_optimizer.inl @@ -148,6 +148,7 @@ void BlockedGraphCutOptimizer block_offset, delta, df, + update_rule, update_field ); @@ -191,6 +192,7 @@ bool BlockedGraphCutOptimizer const int3& block_offset, const float3& delta, // delta in mm const DisplacementField& df, + Settings::UpdateRule update_rule, DisplacementField& update_field ) { @@ -222,7 +224,7 @@ bool BlockedGraphCutOptimizer } float3 d1 = df.get(p); - float3 d1d = df.get(p, delta); + float3 d1d = df.get(p, delta, update_rule == Settings::UpdateRule_Compositive); double f0 = unary_fn(p, d1); double f1 = unary_fn(p, d1d); @@ -275,7 +277,7 @@ bool BlockedGraphCutOptimizer #define ADD_STEP(x_, y_, z_) \ int3 step = int3{x_, y_, z_}; \ float3 d2 = df.get(p+step); \ - float3 d2d = df.get(p+step, delta); \ + float3 d2d = df.get(p+step, delta, update_rule == Settings::UpdateRule_Compositive); \ double f00 = binary_fn(p, d1, d2, step); \ double f01 = binary_fn(p, d1, d2d, step); \ double f10 = binary_fn(p, d1d, d2, step); \ diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index 0830c48..7fbd34e 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -31,21 +31,18 @@ class DisplacementField } // delta : Delta in world space (mm) - inline float3 get(const int3& p, const float3& delta) const + inline float3 get(const int3& p, const float3& delta, bool composite) const { - // if (_update_rule == Settings::UpdateRule_Compositive) { - // // Convert delta to fixed image space - // // TODO: What about orientation? - // float3 fp { - // p.x + delta.x / _df.spacing().x, - // p.y + delta.y / _df.spacing().y, - // p.z + delta.z / _df.spacing().z - // }; - // return _df.linear_at(fp, stk::Border_Replicate) + delta; - // } - // else /*(_update_rule == Settings::UpdateRule_Additive)*/ { + if (composite) { + float3 p1 = _df.index2point(p); + float3 p2 = p1 + delta; + float3 p3 = p2 + _df.linear_at_point(p2, stk::Border_Replicate); + + return p3 - p1; + } + else /*(_update_rule == Settings::UpdateRule_Additive)*/ { return _df(p) + delta; - //} + } } inline void set(const int3& p, const float3& d) From 63c05ecec324ac388f86d889197b0eccde0f3e50 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 4 Nov 2019 13:41:41 +0100 Subject: [PATCH 13/37] Reset update field --- src/deform_lib/registration/blocked_graph_cut_optimizer.inl | 3 +-- src/deform_lib/registration/displacement_field.h | 4 ++++ 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/src/deform_lib/registration/blocked_graph_cut_optimizer.inl b/src/deform_lib/registration/blocked_graph_cut_optimizer.inl index 8b9b213..604a6c6 100644 --- a/src/deform_lib/registration/blocked_graph_cut_optimizer.inl +++ b/src/deform_lib/registration/blocked_graph_cut_optimizer.inl @@ -159,8 +159,7 @@ void BlockedGraphCutOptimizer } df.update(update_field, (update_rule == Settings::UpdateRule_Compositive)); - // TODO: Reset - update_field = DisplacementField(update_field.size()); + update_field.fill(float3{0, 0, 0}); } } } diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index 7fbd34e..9196040 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -71,6 +71,10 @@ class DisplacementField _df(p) += update_field.get(p); } } + + void fill(const float3& v) + { + _df.fill(v); } void set_affine_transform(const AffineTransform& transform) From ee27440a143401b501e96a9fecbe9f73be775588 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 4 Nov 2019 13:43:32 +0100 Subject: [PATCH 14/37] Threaded update --- src/deform_lib/registration/displacement_field.h | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index 9196040..805b090 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -59,7 +59,13 @@ class DisplacementField void update(const DisplacementField& update_field, bool composite) { - for (int3 p : _df.size()) { + dim3 dims = update_field.size(); + + #pragma omp parallel for + for (int z = 0; z < (int)dims.z; ++z) { + for (int y = 0; y < (int)dims.y; ++y) { + for (int x = 0; x < (int)dims.x; ++x) { + int3 p {x, y, z}; if (composite) { float3 p1 = _df.index2point(p); float3 p2 = p1 + update_field.get(p); @@ -70,7 +76,8 @@ class DisplacementField else { _df(p) += update_field.get(p); } - } + }}} + } void fill(const float3& v) { From 41f4de9e54d9b5f7f7b5882df3ff782ef2207ec6 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 4 Nov 2019 13:43:54 +0100 Subject: [PATCH 15/37] Test --- .../registration/displacement_field.h | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index 805b090..1cd76c9 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -33,16 +33,16 @@ class DisplacementField // delta : Delta in world space (mm) inline float3 get(const int3& p, const float3& delta, bool composite) const { - if (composite) { - float3 p1 = _df.index2point(p); - float3 p2 = p1 + delta; - float3 p3 = p2 + _df.linear_at_point(p2, stk::Border_Replicate); - - return p3 - p1; - } - else /*(_update_rule == Settings::UpdateRule_Additive)*/ { + // if (composite) { + // float3 p1 = _df.index2point(p); + // float3 p2 = p1 + delta; + // float3 p3 = p2 + _df.linear_at_point(p2, stk::Border_Replicate); + + // return p3 - p1; + // } + //else /*(_update_rule == Settings::UpdateRule_Additive)*/ { return _df(p) + delta; - } + //} } inline void set(const int3& p, const float3& d) From 28293b420234db9f6bf65021c0ab19747cd8e9b2 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 4 Nov 2019 14:16:34 +0100 Subject: [PATCH 16/37] Affine transform --- .../registration/displacement_field.h | 36 ++++++++++++------- 1 file changed, 23 insertions(+), 13 deletions(-) diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index 1cd76c9..89cb6b4 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -25,26 +25,35 @@ class DisplacementField } ~DisplacementField() {} + // p : Index in displacement field inline float3 get(const int3& p) const { - return _df(p); + return get(_df.index2point(p));//_df(p); + } + + // p : Point in world space + inline float3 get(const float3& p) const + { + return _df.linear_at_point(_affine.transform_point(p), stk::Border_Replicate); } // delta : Delta in world space (mm) inline float3 get(const int3& p, const float3& delta, bool composite) const { - // if (composite) { - // float3 p1 = _df.index2point(p); - // float3 p2 = p1 + delta; - // float3 p3 = p2 + _df.linear_at_point(p2, stk::Border_Replicate); - - // return p3 - p1; - // } - //else /*(_update_rule == Settings::UpdateRule_Additive)*/ { - return _df(p) + delta; - //} + if (composite) { + float3 p1 = _df.index2point(p); + float3 p2 = p1 + delta; + float3 p3 = p2 + get(p2);//_df.linear_at_point(p2, stk::Border_Replicate); + + return p3 - p1; + } + else /*(_update_rule == Settings::UpdateRule_Additive)*/ { + return get(p) + delta; + } } + // Sets the displacement at given index p to value d. + // This modifies the displacement field directly with no regards to the affine transform inline void set(const int3& p, const float3& d) { _df(p) = d; @@ -54,7 +63,8 @@ class DisplacementField // Returns coordinates in world space inline float3 transform_index(const int3& p) const { - return _df.index2point(p) + _df(p); + float3 fp = _df.index2point(p); + return fp + get(fp); } void update(const DisplacementField& update_field, bool composite) @@ -69,7 +79,7 @@ class DisplacementField if (composite) { float3 p1 = _df.index2point(p); float3 p2 = p1 + update_field.get(p); - float3 p3 = p2 + _df.linear_at_point(p2, stk::Border_Replicate); + float3 p3 = p2 + get(p2); _df(p) = p3 - p1; } From 95e8bf7adea71870a2594064ba982d5a83d691b9 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 4 Nov 2019 14:56:32 +0100 Subject: [PATCH 17/37] Fix --- src/deform_lib/registration/affine.cpp | 4 +--- .../registration/displacement_field.h | 21 ++++++++++++++----- 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/src/deform_lib/registration/affine.cpp b/src/deform_lib/registration/affine.cpp index f844438..804ef98 100644 --- a/src/deform_lib/registration/affine.cpp +++ b/src/deform_lib/registration/affine.cpp @@ -42,9 +42,7 @@ const float3& AffineTransform::offset() const } -AffineTransform parse_affine_transform_file( - const std::string& filename -) +AffineTransform parse_affine_transform_file(const std::string& filename) { std::string line; diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index 89cb6b4..e28fb29 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -25,16 +25,20 @@ class DisplacementField } ~DisplacementField() {} + // Returns displacement at index p // p : Index in displacement field inline float3 get(const int3& p) const { - return get(_df.index2point(p));//_df(p); + float3 p2 = _affine.transform_point(transform_index(p)); + return p2 - float3{(float)p.x, (float)p.y, (float)p.z}; } + // Returns displacement at point p // p : Point in world space inline float3 get(const float3& p) const { - return _df.linear_at_point(_affine.transform_point(p), stk::Border_Replicate); + float3 p2 = _affine.transform_point(transform_point(p)); + return p2 - p; } // delta : Delta in world space (mm) @@ -53,18 +57,25 @@ class DisplacementField } // Sets the displacement at given index p to value d. - // This modifies the displacement field directly with no regards to the affine transform + // This modifies the displacement field directly with no regards to the affine + // transformation inline void set(const int3& p, const float3& d) { _df(p) = d; } + // p : Point in world space + // Returns coordinates in world space + inline float3 transform_point(const float3& p) const + { + return p + _df.linear_at_point(p, stk::Border_Replicate); + } + // p : Index in displacement field // Returns coordinates in world space inline float3 transform_index(const int3& p) const { - float3 fp = _df.index2point(p); - return fp + get(fp); + return _df.index2point(p) + _df(p); } void update(const DisplacementField& update_field, bool composite) From 6ca0f5ca1810887e9b1b0f7e587f751e72ffb16f Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 4 Nov 2019 16:13:23 +0100 Subject: [PATCH 18/37] Bugfix --- src/deform_lib/registration/displacement_field.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index e28fb29..0319cb6 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -30,7 +30,7 @@ class DisplacementField inline float3 get(const int3& p) const { float3 p2 = _affine.transform_point(transform_index(p)); - return p2 - float3{(float)p.x, (float)p.y, (float)p.z}; + return p2 - _df.index2point(p); } // Returns displacement at point p @@ -47,7 +47,7 @@ class DisplacementField if (composite) { float3 p1 = _df.index2point(p); float3 p2 = p1 + delta; - float3 p3 = p2 + get(p2);//_df.linear_at_point(p2, stk::Border_Replicate); + float3 p3 = p2 + _df.linear_at_point(p2, stk::Border_Replicate); return p3 - p1; } From b2cf29a06c0267e894f1112e5a8b517e68b09d1b Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 4 Nov 2019 16:35:56 +0100 Subject: [PATCH 19/37] Apply affine in transform_x --- src/deform_lib/registration/displacement_field.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index 0319cb6..a6781e9 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -29,16 +29,14 @@ class DisplacementField // p : Index in displacement field inline float3 get(const int3& p) const { - float3 p2 = _affine.transform_point(transform_index(p)); - return p2 - _df.index2point(p); + return transform_index(p) - _df.index2point(p); } // Returns displacement at point p // p : Point in world space inline float3 get(const float3& p) const { - float3 p2 = _affine.transform_point(transform_point(p)); - return p2 - p; + return transform_point(p) - p; } // delta : Delta in world space (mm) @@ -68,14 +66,16 @@ class DisplacementField // Returns coordinates in world space inline float3 transform_point(const float3& p) const { - return p + _df.linear_at_point(p, stk::Border_Replicate); + return _affine.transform_point( + p + _df.linear_at_point(p, stk::Border_Replicate) + ); } // p : Index in displacement field // Returns coordinates in world space inline float3 transform_index(const int3& p) const { - return _df.index2point(p) + _df(p); + return _affine.transform_point(_df.index2point(p) + _df(p)); } void update(const DisplacementField& update_field, bool composite) From 5c5ca5be65e13cf700168f129568cd557c7089b6 Mon Sep 17 00:00:00 2001 From: Simon Date: Tue, 5 Nov 2019 12:32:10 +0100 Subject: [PATCH 20/37] Temp --- src/deform/transform.cpp | 2 +- .../blocked_graph_cut_optimizer.inl | 7 ++- .../registration/displacement_field.h | 43 ++++++++++++++----- 3 files changed, 40 insertions(+), 12 deletions(-) diff --git a/src/deform/transform.cpp b/src/deform/transform.cpp index cc1810f..b80f130 100644 --- a/src/deform/transform.cpp +++ b/src/deform/transform.cpp @@ -11,7 +11,7 @@ bool TransformCommand::_parse_arguments(void) { _args.add_positional("command", "registration, transform, regularize, jacobian"); _args.add_positional("source", "Path to the image you want to transform"); - _args.add_positional("displacement", "Path to the displacement field used to transform"); + _args.add_positional("transform", "Path to the transform file, could be either affine or a displacement field"); _args.add_positional("output", "Path to the resulting file"); _args.add_option("interp", "-i, --interp", "Interpolation to use, either 'nn' or 'linear' (default)"); diff --git a/src/deform_lib/registration/blocked_graph_cut_optimizer.inl b/src/deform_lib/registration/blocked_graph_cut_optimizer.inl index 604a6c6..18be1ab 100644 --- a/src/deform_lib/registration/blocked_graph_cut_optimizer.inl +++ b/src/deform_lib/registration/blocked_graph_cut_optimizer.inl @@ -331,7 +331,12 @@ bool BlockedGraphCutOptimizer } if (solver.get_var(sub_x, sub_y, sub_z) == 1) { - update_field.set(int3{gx, gy, gz}, delta); + if (update_rule == Settings::UpdateRule_Compositive) { + update_field.set(int3{gx, gy, gz}, delta); + } + else { + update_field.set(int3{gx, gy, gz}, delta); + } changed_flag = true; } } diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index a6781e9..9150b50 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -15,8 +15,12 @@ class DisplacementField { public: DisplacementField() {} - DisplacementField(const stk::VolumeFloat3& df) : - _df(df) + DisplacementField( + const stk::VolumeFloat3& df, + const AffineTransform& affine = AffineTransform() + ) : + _df(df), + _affine(affine) { } DisplacementField(const dim3& dims) : @@ -45,7 +49,7 @@ class DisplacementField if (composite) { float3 p1 = _df.index2point(p); float3 p2 = p1 + delta; - float3 p3 = p2 + _df.linear_at_point(p2, stk::Border_Replicate); + float3 p3 = transform_point(p2); return p3 - p1; } @@ -66,22 +70,41 @@ class DisplacementField // Returns coordinates in world space inline float3 transform_point(const float3& p) const { - return _affine.transform_point( + + + float3 d1 = p + _df.linear_at_point(p, stk::Border_Replicate); + + float3 d2 = _affine.transform_point( p + _df.linear_at_point(p, stk::Border_Replicate) ); + + // if (abs(d1.x - d2.x) > 0.000000001f + // || abs(d1.y - d2.y) > 0.000000001f + // || abs(d1.z - d2.z) > 0.000000001f) { + // LOG(Info) << "pt: " << d1 << " != " << d2; + // LOG(Info) << p; + // } + + return d2; } // p : Index in displacement field // Returns coordinates in world space inline float3 transform_index(const int3& p) const { - return _affine.transform_point(_df.index2point(p) + _df(p)); + + float3 d1 = _df.index2point(p) + _df(p); + float3 d2 = _affine.transform_point(_df.index2point(p) + _df(p)); + + return d2; } void update(const DisplacementField& update_field, bool composite) { dim3 dims = update_field.size(); + DisplacementField buffer = this->clone(); + #pragma omp parallel for for (int z = 0; z < (int)dims.z; ++z) { for (int y = 0; y < (int)dims.y; ++y) { @@ -90,7 +113,7 @@ class DisplacementField if (composite) { float3 p1 = _df.index2point(p); float3 p2 = p1 + update_field.get(p); - float3 p3 = p2 + get(p2); + float3 p3 = buffer.transform_point(p2); _df(p) = p3 - p1; } @@ -110,14 +133,14 @@ class DisplacementField _affine = transform; } - dim3 size() const + DisplacementField clone() const { - return _df.size(); + return DisplacementField(_df.clone(), _affine); } - void copy_from(const DisplacementField& other) + dim3 size() const { - _df.copy_from(other._df); + return _df.size(); } // Volume containing the displacements only From e99cbc33c7ba49f24c7bbf8d95868f4b327be155 Mon Sep 17 00:00:00 2001 From: Simon Date: Wed, 6 Nov 2019 13:49:11 +0100 Subject: [PATCH 21/37] Remove `regularize_initial_displacement` --- README.md | 5 --- .../blocked_graph_cut_optimizer.h | 4 +-- .../blocked_graph_cut_optimizer.inl | 4 +-- .../registration/cost_functions/regularizer.h | 12 +------ .../gpu/cost_functions/binary_function.cu | 31 +++++++++---------- .../gpu/cost_functions/binary_function.h | 10 ------ .../registration/gpu_registration_engine.cpp | 10 ------ .../registration/registration_engine.cpp | 26 ++++------------ src/deform_lib/registration/settings.cpp | 5 +-- src/deform_lib/registration/settings.h | 4 --- test/test_gpu_cost.cpp | 1 - test/test_settings.cpp | 3 -- 12 files changed, 26 insertions(+), 89 deletions(-) diff --git a/README.md b/README.md index 5e81559..abb864e 100644 --- a/README.md +++ b/README.md @@ -104,7 +104,6 @@ To perform a registration using the standalone executable ```yaml pyramid_levels: 6 pyramid_stop_level: 0 -regularize_initial_displacement: false solver: gco constraints_weight: 1000.0 landmarks_weight: 1.0 @@ -159,10 +158,6 @@ halves the resolution of the input volumes. Setting `pyramid_stop_level` to > 0 specifies that the registration should not be run on the original resolution (level 0). -`regularize_initial_displacement` controls whether to include the initial -displacement in the regularization. This is disabled by default, meaning the -initial displacement have no effect on the regularization. - `solver` selects which solver to use for the energy minimization. Available solvers are `gco`, `gridcut`, and `icm`. Note: deform needs to be compiled with `DF_ENABLE_GCO` and `DF_ENABLE_GRIDCUT` for `gco` and `gridcut` to be enabled. diff --git a/src/deform_lib/registration/blocked_graph_cut_optimizer.h b/src/deform_lib/registration/blocked_graph_cut_optimizer.h index 8549aeb..2bae47d 100644 --- a/src/deform_lib/registration/blocked_graph_cut_optimizer.h +++ b/src/deform_lib/registration/blocked_graph_cut_optimizer.h @@ -28,8 +28,8 @@ class BlockedGraphCutOptimizer UnaryFunction& unary_fn, Regularizer& binary_fn, float3 step_size, - DisplacementField& df, - Settings::UpdateRule update_rule + Settings::UpdateRule update_rule, + DisplacementField& df ); private: diff --git a/src/deform_lib/registration/blocked_graph_cut_optimizer.inl b/src/deform_lib/registration/blocked_graph_cut_optimizer.inl index 18be1ab..f42a308 100644 --- a/src/deform_lib/registration/blocked_graph_cut_optimizer.inl +++ b/src/deform_lib/registration/blocked_graph_cut_optimizer.inl @@ -32,8 +32,8 @@ void BlockedGraphCutOptimizer UnaryFunction& unary_fn, Regularizer& binary_fn, float3 step_size, - DisplacementField& df, - Settings::UpdateRule update_rule + Settings::UpdateRule update_rule, + DisplacementField& df ) { dim3 dims = df.size(); diff --git a/src/deform_lib/registration/cost_functions/regularizer.h b/src/deform_lib/registration/cost_functions/regularizer.h index 3cfa2f8..2c38675 100644 --- a/src/deform_lib/registration/cost_functions/regularizer.h +++ b/src/deform_lib/registration/cost_functions/regularizer.h @@ -41,14 +41,6 @@ struct Regularizer _spacing = spacing; } - // Sets the initial displacement for this registration level. This will be - // the reference when computing the regularization energy. Any displacement - // identical to the initial displacement will result in zero energy. - void set_initial_displacement(const DisplacementField& initial) - { - _initial = initial; - } - #ifdef DF_ENABLE_REGULARIZATION_WEIGHT_MAP void set_weight_map(const stk::VolumeFloat& map) { _weight_map = map; } #endif // DF_ENABLE_REGULARIZATION_WEIGHT_MAP @@ -67,7 +59,7 @@ struct Regularizer }; // The diff should be relative to the initial displacement diff - float3 diff = (def0-_initial.get(p)) - (def1-_initial.get(p+step)); + float3 diff = def0 - def1; float dist_squared = stk::norm2(diff); float step_squared = stk::norm2(step_in_mm); @@ -94,8 +86,6 @@ struct Regularizer float _half_exponent; float3 _spacing; - DisplacementField _initial; - #ifdef DF_ENABLE_REGULARIZATION_WEIGHT_MAP stk::VolumeFloat _weight_map; #endif // DF_ENABLE_REGULARIZATION_WEIGHT_MAP diff --git a/src/deform_lib/registration/gpu/cost_functions/binary_function.cu b/src/deform_lib/registration/gpu/cost_functions/binary_function.cu index c97ec56..9e6931e 100644 --- a/src/deform_lib/registration/gpu/cost_functions/binary_function.cu +++ b/src/deform_lib/registration/gpu/cost_functions/binary_function.cu @@ -30,7 +30,6 @@ __device__ float4 energy( template __global__ void regularizer_kernel( TDisplacementField df, - TDisplacementField initial_df, float4 delta, float weight, float scale, @@ -62,8 +61,8 @@ __global__ void regularizer_kernel( // Cost ordered as E00, E01, E10, E11 - float4 d0 = df.get(p) - initial_df.get(p); - float4 d1 = df.get(p, delta) - initial_df.get(p); + float4 d0 = df.get(p); + float4 d1 = df.get(p, delta); float4 o_x = {0, 0, 0, 0}; float4 o_y = {0, 0, 0, 0}; @@ -71,8 +70,8 @@ __global__ void regularizer_kernel( if (gx + 1 < (int) df.size().x) { int3 step {1, 0, 0}; - float4 dn0 = df.get(p+step) - initial_df.get(p+step); - float4 dn1 = df.get(p+step, delta) - initial_df.get(p+step); + float4 dn0 = df.get(p+step); + float4 dn1 = df.get(p+step, delta); o_x = energy( d0, @@ -85,8 +84,8 @@ __global__ void regularizer_kernel( } if (gy + 1 < (int) df.size().y) { int3 step {0, 1, 0}; - float4 dn0 = df.get(p+step) - initial_df.get(p+step); - float4 dn1 = df.get(p+step, delta) - initial_df.get(p+step); + float4 dn0 = df.get(p+step); + float4 dn1 = df.get(p+step, delta); o_y = energy( d0, @@ -99,8 +98,8 @@ __global__ void regularizer_kernel( } if (gz + 1 < (int) df.size().z) { int3 step {0, 0, 1}; - float4 dn0 = df.get(p+step) - initial_df.get(p+step); - float4 dn1 = df.get(p+step, delta) - initial_df.get(p+step); + float4 dn0 = df.get(p+step); + float4 dn1 = df.get(p+step, delta); o_z = energy( d0, @@ -120,8 +119,8 @@ __global__ void regularizer_kernel( if (x == 0 && gx != 0) { int3 step {-1, 0, 0}; - float4 dn0 = df.get(p+step) - initial_df.get(p+step); - float4 dn1 = df.get(p+step, delta) - initial_df.get(p+step); + float4 dn0 = df.get(p+step); + float4 dn1 = df.get(p+step, delta); float4 e = energy( d0, @@ -143,8 +142,8 @@ __global__ void regularizer_kernel( if (y == 0 && gy != 0) { int3 step {0, -1, 0}; - float4 dn0 = df.get(p+step) - initial_df.get(p+step); - float4 dn1 = df.get(p+step, delta) - initial_df.get(p+step); + float4 dn0 = df.get(p+step); + float4 dn1 = df.get(p+step, delta); float4 e = energy( d0, @@ -163,8 +162,8 @@ __global__ void regularizer_kernel( if (z == 0 && gz != 0) { int3 step {0, 0, -1}; - float4 dn0 = df.get(p+step) - initial_df.get(p+step); - float4 dn1 = df.get(p+step, delta) - initial_df.get(p+step); + float4 dn0 = df.get(p+step); + float4 dn1 = df.get(p+step, delta); float4 e = energy( d0, @@ -227,7 +226,6 @@ void GpuBinaryFunction::operator()( regularizer_kernel> <<>>( df, - _initial, delta4, _weight, _scale, @@ -244,7 +242,6 @@ void GpuBinaryFunction::operator()( regularizer_kernel> <<>>( df, - _initial, delta4, _weight, _scale, diff --git a/src/deform_lib/registration/gpu/cost_functions/binary_function.h b/src/deform_lib/registration/gpu/cost_functions/binary_function.h index 4108f55..c0ffcea 100644 --- a/src/deform_lib/registration/gpu/cost_functions/binary_function.h +++ b/src/deform_lib/registration/gpu/cost_functions/binary_function.h @@ -35,14 +35,6 @@ class GpuBinaryFunction _spacing = spacing; } - // Sets the initial displacement for this registration level. This will be - // the reference when computing the regularization energy. Any displacement - // identical to the initial displacement will result in zero energy. - void set_initial_displacement(const GpuDisplacementField& initial) - { - _initial = initial; - } - // Computes the regularization cost in three directions (x+, y+, z+), with and without // applied delta. Results are stored into the three provided cost volumes (of type float2) // df : Displacement field @@ -69,7 +61,5 @@ class GpuBinaryFunction float _scale; float _half_exponent; float3 _spacing; - - GpuDisplacementField _initial; }; diff --git a/src/deform_lib/registration/gpu_registration_engine.cpp b/src/deform_lib/registration/gpu_registration_engine.cpp index f8240f6..c388a00 100644 --- a/src/deform_lib/registration/gpu_registration_engine.cpp +++ b/src/deform_lib/registration/gpu_registration_engine.cpp @@ -154,16 +154,6 @@ void GpuRegistrationEngine::build_binary_function(int level, GpuBinaryFunction& binary_fn.set_regularization_weight(_settings.levels[level].regularization_weight); binary_fn.set_regularization_scale(_settings.levels[level].regularization_scale); binary_fn.set_regularization_exponent(_settings.levels[level].regularization_exponent); - - stk::GpuVolume df = _deformation_pyramid.volume(level); - if (!_settings.regularize_initial_displacement) { - // Clone the def, because the current copy will be changed when executing the optimizer - binary_fn.set_initial_displacement(df.clone()); - } - else { - stk::VolumeFloat4 zeros(df.size(), float4{0,0,0,0}); - binary_fn.set_initial_displacement(stk::GpuVolume(zeros)); - } } GpuRegistrationEngine::GpuRegistrationEngine(const Settings& settings) : diff --git a/src/deform_lib/registration/registration_engine.cpp b/src/deform_lib/registration/registration_engine.cpp index ac8f829..4df0f52 100644 --- a/src/deform_lib/registration/registration_engine.cpp +++ b/src/deform_lib/registration/registration_engine.cpp @@ -217,20 +217,6 @@ void RegistrationEngine::build_regularizer(int level, Regularizer& binary_fn) if (_regularization_weight_map.volume(level).valid()) binary_fn.set_weight_map(_regularization_weight_map.volume(level)); #endif - - stk::Volume df = _deformation_pyramid.volume(level); - if (!_settings.regularize_initial_displacement) { - // Clone the def, because the current copy will be changed when executing the optimizer - binary_fn.set_initial_displacement(DisplacementField(df.clone())); - } - else { - binary_fn.set_initial_displacement(stk::VolumeFloat3(df.size(), float3{0,0,0})); - } - - #ifdef DF_ENABLE_REGULARIZATION_WEIGHT_MAP - if (_regularization_weight_map.volume(level).valid()) - binary_fn.set_weight_map(_regularization_weight_map.volume(level)); - #endif } void RegistrationEngine::build_unary_function(int level, UnaryFunction& unary_fn) @@ -652,8 +638,8 @@ stk::Volume RegistrationEngine::execute() unary_fn, binary_fn, _settings.levels[l].step_size, - df, - _settings.update_rule + _settings.update_rule, + df ); } #if defined(DF_ENABLE_GCO) @@ -668,8 +654,8 @@ stk::Volume RegistrationEngine::execute() unary_fn, binary_fn, _settings.levels[l].step_size, - df, - _settings.update_rule + _settings.update_rule, + df ); } #endif @@ -685,8 +671,8 @@ stk::Volume RegistrationEngine::execute() unary_fn, binary_fn, _settings.levels[l].step_size, - df, - _settings.update_rule + _settings.update_rule, + df ); } #endif diff --git a/src/deform_lib/registration/settings.cpp b/src/deform_lib/registration/settings.cpp index cb8ee58..b73459f 100644 --- a/src/deform_lib/registration/settings.cpp +++ b/src/deform_lib/registration/settings.cpp @@ -21,8 +21,6 @@ pyramid_levels: 6 pyramid_stop_level: 0 -regularize_initial_displacement: false - solver: gco update_rule: additive @@ -372,7 +370,6 @@ void print_registration_settings(const Settings& settings, std::ostream& s) s << "pyramid_stop_level = " << settings.pyramid_stop_level << std::endl; s << "num_pyramid_levels = " << settings.num_pyramid_levels << std::endl; s << "landmarks_stop_level = " << settings.landmarks_stop_level << std::endl; - s << "regularize_initial_displacement = " << settings.regularize_initial_displacement << std::endl; s << "solver = " << solver_to_str(settings.solver) << std::endl; s << "update_rule = " << update_rule_to_str(settings.update_rule) << std::endl; @@ -442,7 +439,7 @@ void parse_registration_settings(const std::string& str, Settings& settings) settings.landmarks_stop_level = value.as(); } else if (key == "regularize_initial_displacement") { - settings.regularize_initial_displacement = value.as(); + LOG(Warning) << "'regularize_initial_displacement' is no longer supported"; } else if (key == "solver") { std::string solver = value.as(); diff --git a/src/deform_lib/registration/settings.h b/src/deform_lib/registration/settings.h index f38bc41..4841223 100644 --- a/src/deform_lib/registration/settings.h +++ b/src/deform_lib/registration/settings.h @@ -134,9 +134,6 @@ struct Settings std::vector image_slots; - // Should the regularization also apply on the initial displacement? - bool regularize_initial_displacement; - // What solver to use for the energy minimization Solver solver; @@ -147,7 +144,6 @@ struct Settings pyramid_stop_level(0), num_pyramid_levels(6), landmarks_stop_level(0), - regularize_initial_displacement(false), #if defined(DF_ENABLE_GCO) solver(Solver_GCO), #elif defined(DF_ENABLE_GRIDCUT) diff --git a/test/test_gpu_cost.cpp b/test/test_gpu_cost.cpp index 06ae1cb..1106e2f 100644 --- a/test/test_gpu_cost.cpp +++ b/test/test_gpu_cost.cpp @@ -159,7 +159,6 @@ TEST_CASE("gpu_cost_regularizer", "") fill_data(df); Regularizer cpu_fn(1.0f); - cpu_fn.set_initial_displacement(zero_df); for (int z = 0; z < (int)dims.z; ++z) { for (int y = 0; y < (int)dims.y; ++y) { diff --git a/test/test_settings.cpp b/test/test_settings.cpp index 4ce54b2..de6660a 100644 --- a/test/test_settings.cpp +++ b/test/test_settings.cpp @@ -22,8 +22,6 @@ regularization_exponent: 1.5 constraints_weight: 1234.1234 -regularize_initial_displacement: true - solver: icm update_rule: compositive @@ -126,7 +124,6 @@ TEST_CASE("parse_registration_settings") REQUIRE(settings.pyramid_stop_level == 2); REQUIRE(settings.num_pyramid_levels == 4); - REQUIRE(settings.regularize_initial_displacement == true); REQUIRE(settings.solver == Settings::Solver_ICM); REQUIRE(settings.update_rule == Settings::UpdateRule_Compositive); From f9255d3a7d1aaa3df22ac246451d192c34758d1c Mon Sep 17 00:00:00 2001 From: Simon Date: Wed, 6 Nov 2019 13:59:27 +0100 Subject: [PATCH 22/37] Fixed test --- test/test_affine.cpp | 40 ---------------------------------------- 1 file changed, 40 deletions(-) diff --git a/test/test_affine.cpp b/test/test_affine.cpp index 0f4c9dd..8e350e0 100644 --- a/test/test_affine.cpp +++ b/test/test_affine.cpp @@ -144,44 +144,4 @@ TEST_CASE("parse_affine") CHECK(8.0f == pt3.y); CHECK(11.0f == pt3.z); } - SECTION("test_displacement_vector") - { - // Identity transform - AffineTransform transform1( - Matrix3x3f::Identity, - float3{0, 0, 0} - ); - - float3 d1 = transform1.displacement_vector(float3{2,3,4}); - CHECK(0.0f == d1.x); - CHECK(0.0f == d1.y); - CHECK(0.0f == d1.z); - - // Scaling - AffineTransform transform2( - Matrix3x3f{ - float3{2, 0, 0}, - float3{0, 3, 0}, - float3{0, 0, 4} - }, - float3{0, 0, 0} - ); - - float3 pt2 = transform2.displacement_vector(float3{1, 2, 3}); - CHECK(1.0f == pt2.x); - CHECK(4.0f == pt2.y); - CHECK(9.0f == pt2.z); - - // Translation - AffineTransform transform3( - Matrix3x3f::Identity, - float3{4, 6, 8} - ); - - float3 pt3 = transform3.displacement_vector(float3{1, 2, 3}); - CHECK(4.0f == pt3.x); - CHECK(6.0f == pt3.y); - CHECK(8.0f == pt3.z); - - } } From 49183f5959683d845049bd846cca962c6e501514 Mon Sep 17 00:00:00 2001 From: Simon Date: Wed, 6 Nov 2019 17:13:40 +0100 Subject: [PATCH 23/37] Optional affine transform added to transform command --- src/deform/transform.cpp | 22 ++++++++++++++----- src/deform_lib/CMakeLists.txt | 4 ++-- .../{affine.cpp => affine_transform.cpp} | 4 ++-- .../{affine.h => affine_transform.h} | 0 .../registration/displacement_field.h | 2 +- src/deform_lib/registration/transform.cpp | 4 ++-- test/test_affine.cpp | 2 +- 7 files changed, 24 insertions(+), 14 deletions(-) rename src/deform_lib/registration/{affine.cpp => affine_transform.cpp} (98%) rename src/deform_lib/registration/{affine.h => affine_transform.h} (100%) diff --git a/src/deform/transform.cpp b/src/deform/transform.cpp index b80f130..7153a46 100644 --- a/src/deform/transform.cpp +++ b/src/deform/transform.cpp @@ -3,6 +3,8 @@ #include #include "deform_lib/arg_parser.h" +#include "deform_lib/registration/affine_transform.h" +#include "deform_lib/registration/displacement_field.h" #include "deform_lib/registration/transform.h" #include "deform/command.h" @@ -11,10 +13,11 @@ bool TransformCommand::_parse_arguments(void) { _args.add_positional("command", "registration, transform, regularize, jacobian"); _args.add_positional("source", "Path to the image you want to transform"); - _args.add_positional("transform", "Path to the transform file, could be either affine or a displacement field"); + _args.add_positional("displacement", "Path to the displacement field"); _args.add_positional("output", "Path to the resulting file"); _args.add_option("interp", "-i, --interp", "Interpolation to use, either 'nn' or 'linear' (default)"); + _args.add_option("affine", "-a, --affine", "Path to an (optional) affine transformation file"); return _args.parse(); } @@ -41,18 +44,25 @@ int TransformCommand::_execute(void) LOG(Info) << "Interpolation method: " << ((interp == transform::Interp_Linear) ? "linear" : "nn"); LOG(Info) << "Input: '" << _args.positional("source") << "'"; - LOG(Info) << "Displacement: '" << _args.positional("displacement") << "'"; + LOG(Info) << "Displacement field: '" << _args.positional("displacement") << "'"; + + if (_args.is_set("affine")) + LOG(Info) << "Affine transform: '" << _args.option("affine") << "'"; stk::Volume src = stk::read_volume(_args.positional("source").c_str()); if (!src.valid()) return EXIT_FAILURE; - stk::Volume def = stk::read_volume(_args.positional("displacement").c_str()); - if (!def.valid()) + stk::Volume df = stk::read_volume(_args.positional("displacement").c_str()); + if (!df.valid()) return EXIT_FAILURE; - ASSERT(def.voxel_type() == stk::Type_Float3); + ASSERT(df.voxel_type() == stk::Type_Float3); - stk::Volume result = transform_volume(src, def, interp); + AffineTransform affine; + if (_args.is_set("affine")) + affine = parse_affine_transform_file(_args.option("affine")); + + stk::Volume result = transform_volume(src, DisplacementField(df, affine), interp); if (!result.valid()) return EXIT_FAILURE; diff --git a/src/deform_lib/CMakeLists.txt b/src/deform_lib/CMakeLists.txt index 64ca77a..33ef5e8 100644 --- a/src/deform_lib/CMakeLists.txt +++ b/src/deform_lib/CMakeLists.txt @@ -16,8 +16,8 @@ set(DEFORM_LIB_SRCS "registration/cost_functions/sub_function.h" "registration/cost_functions/unary_function.h" - "registration/affine.cpp" - "registration/affine.h" + "registration/affine_transform.cpp" + "registration/affine_transform.h" "registration/block_change_flags.cpp" "registration/block_change_flags.h" "registration/blocked_graph_cut_optimizer.h" diff --git a/src/deform_lib/registration/affine.cpp b/src/deform_lib/registration/affine_transform.cpp similarity index 98% rename from src/deform_lib/registration/affine.cpp rename to src/deform_lib/registration/affine_transform.cpp index 804ef98..f0816d8 100644 --- a/src/deform_lib/registration/affine.cpp +++ b/src/deform_lib/registration/affine_transform.cpp @@ -1,4 +1,4 @@ -#include "affine.h" +#include "affine_transform.h" #include #include @@ -89,7 +89,7 @@ AffineTransform parse_affine_transform_file(const std::string& filename) // Parameters: *12 std::istringstream ss(keyvalues["Parameters"]); - + float* pmatrix = &matrix._rows[0].x; for (int i = 0; i < 9; ++i) { ss >> pmatrix[i]; diff --git a/src/deform_lib/registration/affine.h b/src/deform_lib/registration/affine_transform.h similarity index 100% rename from src/deform_lib/registration/affine.h rename to src/deform_lib/registration/affine_transform.h diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index 9150b50..430bd41 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -3,7 +3,7 @@ #include #include -#include "affine.h" +#include "affine_transform.h" #include "settings.h" /** Wrapper for a displacement field image. diff --git a/src/deform_lib/registration/transform.cpp b/src/deform_lib/registration/transform.cpp index 1dc7e9f..2c8901c 100644 --- a/src/deform_lib/registration/transform.cpp +++ b/src/deform_lib/registration/transform.cpp @@ -68,8 +68,8 @@ namespace for (int z = 0; z < int(dims.z); ++z) { for (int y = 0; y < int(dims.y); ++y) { for (int x = 0; x < int(dims.x); ++x) { - const float3 moving_p = df.transform_index(int3{x,y,z}); - out(x, y, z) = src.linear_at_point(moving_p, stk::Border_Constant); + const float3 pt = df.transform_index(int3{x,y,z}); + out(x, y, z) = src.linear_at_point(pt, stk::Border_Constant); } } } diff --git a/test/test_affine.cpp b/test/test_affine.cpp index 8e350e0..697f80f 100644 --- a/test/test_affine.cpp +++ b/test/test_affine.cpp @@ -1,6 +1,6 @@ #include "catch.hpp" -#include +#include #include #include From cd0cc3cab4c6a8d7579276c5e7ec51e980c61f60 Mon Sep 17 00:00:00 2001 From: Simon Date: Wed, 6 Nov 2019 17:24:00 +0100 Subject: [PATCH 24/37] Cleanup --- src/deform_lib/CMakeLists.txt | 2 + .../registration/displacement_field.cpp | 67 ++++++++++++ .../registration/displacement_field.h | 103 +++--------------- 3 files changed, 87 insertions(+), 85 deletions(-) create mode 100644 src/deform_lib/registration/displacement_field.cpp diff --git a/src/deform_lib/CMakeLists.txt b/src/deform_lib/CMakeLists.txt index 33ef5e8..35ec827 100644 --- a/src/deform_lib/CMakeLists.txt +++ b/src/deform_lib/CMakeLists.txt @@ -22,6 +22,8 @@ set(DEFORM_LIB_SRCS "registration/block_change_flags.h" "registration/blocked_graph_cut_optimizer.h" "registration/blocked_graph_cut_optimizer.inl" + "registration/displacement_field.cpp" + "registration/displacement_field.h" "registration/landmarks.cpp" "registration/landmarks.h" "registration/registration_engine.cpp" diff --git a/src/deform_lib/registration/displacement_field.cpp b/src/deform_lib/registration/displacement_field.cpp new file mode 100644 index 0000000..9ea5291 --- /dev/null +++ b/src/deform_lib/registration/displacement_field.cpp @@ -0,0 +1,67 @@ +#include "displacement_field.h" + + + +DisplacementField::DisplacementField() {} +DisplacementField::DisplacementField( + const stk::VolumeFloat3& df, + const AffineTransform& affine +) : + _df(df), + _affine(affine) +{ +} +DisplacementField::DisplacementField(const dim3& dims) : + _df(dims, float3{0, 0, 0}) +{ +} +DisplacementField::~DisplacementField() {} + +void DisplacementField::update(const DisplacementField& update_field, bool composite) +{ + dim3 dims = update_field.size(); + + DisplacementField buffer = this->clone(); + + #pragma omp parallel for + for (int z = 0; z < (int)dims.z; ++z) { + for (int y = 0; y < (int)dims.y; ++y) { + for (int x = 0; x < (int)dims.x; ++x) { + int3 p {x, y, z}; + if (composite) { + float3 p1 = _df.index2point(p); + float3 p2 = p1 + update_field.get(p); + float3 p3 = buffer.transform_point(p2); + + _df(p) = p3 - p1; + } + else { + _df(p) += update_field.get(p); + } + }}} +} +void DisplacementField::fill(const float3& v) +{ + _df.fill(v); +} +DisplacementField DisplacementField::clone() const +{ + return DisplacementField(_df.clone(), _affine); +} +dim3 DisplacementField::size() const +{ + return _df.size(); +} +const stk::VolumeFloat3& DisplacementField::volume() const +{ + return _df; +} +stk::VolumeFloat3& DisplacementField::volume() +{ + return _df; +} +bool DisplacementField::valid() const +{ + return _df.valid(); +} + diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index 430bd41..01c72c7 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -14,20 +14,14 @@ class DisplacementField { public: - DisplacementField() {} + DisplacementField(); DisplacementField( const stk::VolumeFloat3& df, const AffineTransform& affine = AffineTransform() - ) : - _df(df), - _affine(affine) - { - } - DisplacementField(const dim3& dims) : - _df(dims, float3{0, 0, 0}) - { - } - ~DisplacementField() {} + ); + // Creates a new identity displacement field of size dims + DisplacementField(const dim3& dims); + ~DisplacementField(); // Returns displacement at index p // p : Index in displacement field @@ -53,7 +47,7 @@ class DisplacementField return p3 - p1; } - else /*(_update_rule == Settings::UpdateRule_Additive)*/ { + else { return get(p) + delta; } } @@ -70,101 +64,40 @@ class DisplacementField // Returns coordinates in world space inline float3 transform_point(const float3& p) const { - - - float3 d1 = p + _df.linear_at_point(p, stk::Border_Replicate); - - float3 d2 = _affine.transform_point( + return _affine.transform_point( p + _df.linear_at_point(p, stk::Border_Replicate) ); - - // if (abs(d1.x - d2.x) > 0.000000001f - // || abs(d1.y - d2.y) > 0.000000001f - // || abs(d1.z - d2.z) > 0.000000001f) { - // LOG(Info) << "pt: " << d1 << " != " << d2; - // LOG(Info) << p; - // } - - return d2; } // p : Index in displacement field // Returns coordinates in world space inline float3 transform_index(const int3& p) const { - - float3 d1 = _df.index2point(p) + _df(p); - float3 d2 = _affine.transform_point(_df.index2point(p) + _df(p)); - - return d2; + return _affine.transform_point(_df.index2point(p) + _df(p)); } - void update(const DisplacementField& update_field, bool composite) - { - dim3 dims = update_field.size(); - - DisplacementField buffer = this->clone(); - - #pragma omp parallel for - for (int z = 0; z < (int)dims.z; ++z) { - for (int y = 0; y < (int)dims.y; ++y) { - for (int x = 0; x < (int)dims.x; ++x) { - int3 p {x, y, z}; - if (composite) { - float3 p1 = _df.index2point(p); - float3 p2 = p1 + update_field.get(p); - float3 p3 = buffer.transform_point(p2); - - _df(p) = p3 - p1; - } - else { - _df(p) += update_field.get(p); - } - }}} - } - - void fill(const float3& v) - { - _df.fill(v); - } + void update(const DisplacementField& update_field, bool composite); - void set_affine_transform(const AffineTransform& transform) - { - _affine = transform; - } + // Fills the displacement field component with vector v + void fill(const float3& v); - DisplacementField clone() const - { - return DisplacementField(_df.clone(), _affine); - } + // Clones the displacement field and any affine transformation + DisplacementField clone() const; - dim3 size() const - { - return _df.size(); - } + // Size of the displacement field + dim3 size() const; // Volume containing the displacements only - const stk::VolumeFloat3& volume() const - { - return _df; - } + const stk::VolumeFloat3& volume() const; // Volume containing the displacements only - stk::VolumeFloat3& volume() - { - return _df; - } + stk::VolumeFloat3& volume(); // Returns true if the volume is allocated and ready for use - bool valid() const - { - return _df.valid(); - } - + bool valid() const; private: stk::VolumeFloat3 _df; - AffineTransform _affine; }; From 4e168b93a9d75eb798a284054da300620c94eb0c Mon Sep 17 00:00:00 2001 From: Simon Date: Thu, 7 Nov 2019 10:34:47 +0100 Subject: [PATCH 25/37] Pretty --- src/deform_lib/registration/registration_engine.cpp | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/deform_lib/registration/registration_engine.cpp b/src/deform_lib/registration/registration_engine.cpp index 4df0f52..7f9c8fd 100644 --- a/src/deform_lib/registration/registration_engine.cpp +++ b/src/deform_lib/registration/registration_engine.cpp @@ -47,14 +47,12 @@ namespace ASSERT(df.size() == mask.size() && df.size() == values.size()); dim3 dims = df.size(); for (int z = 0; z < int(dims.z); ++z) { - for (int y = 0; y < int(dims.y); ++y) { - for (int x = 0; x < int(dims.x); ++x) { - if (mask(x, y, z) > 0) { - df.set(int3{x,y,z}, values(x, y, z)); - } - } + for (int y = 0; y < int(dims.y); ++y) { + for (int x = 0; x < int(dims.x); ++x) { + if (mask(x, y, z) > 0) { + df.set(int3{x,y,z}, values(x, y, z)); } - } + }}} } template From 90c1fc56e61b3bb8b74607ab589776f9f93b4818 Mon Sep 17 00:00:00 2001 From: Simon Date: Thu, 7 Nov 2019 14:59:40 +0100 Subject: [PATCH 26/37] Affine transform added to command-line interface --- src/deform/registration.cpp | 13 ++++++++++++- src/deform_lib/registration/registration.cpp | 3 +++ src/deform_lib/registration/registration.h | 3 +++ src/deform_lib/registration/registration_engine.cpp | 4 ++++ src/deform_lib/registration/registration_engine.h | 3 +++ 5 files changed, 25 insertions(+), 1 deletion(-) diff --git a/src/deform/registration.cpp b/src/deform/registration.cpp index 1d9b3b1..32da9b0 100644 --- a/src/deform/registration.cpp +++ b/src/deform/registration.cpp @@ -2,6 +2,7 @@ #include +#include #include #include #include @@ -24,7 +25,8 @@ bool RegistrationCommand::_parse_arguments(void) _args.add_option("fixed_mask", "-fm, --fixed-mask", "Path to the fixed image mask"); _args.add_option("moving_mask", "-mm, --moving-mask", "Path to the moving image mask"); _args.add_group(); - _args.add_option("init_deform", "-d0", "Path to the initial displacement field"); + _args.add_option("init_deform", "-d0", "Path to an initial displacement field"); + _args.add_option("affine", "-a,--affine", "Path to an initial affine transformation"); _args.add_group(); _args.add_option("fixed_points", "-fp, --fixed-points", "Path to the fixed landmark points"); _args.add_option("moving_points", "-mp, --moving-points", "Path to the moving landmark points"); @@ -157,6 +159,14 @@ int RegistrationCommand::_execute(void) return EXIT_FAILURE; } + std::string affine_file = _args.get("affine", ""); + LOG(Info) << "Affine transform: '" << affine_file << "'"; + + AffineTransform affine_transform; + if (!affine_file.empty()) { + affine_transform = parse_affine_transform_file(affine_file); + } + // Constraints std::string constraint_mask_file = _args.get("constraint_mask", ""); std::string constraint_values_file = _args.get("constraint_values", ""); @@ -229,6 +239,7 @@ int RegistrationCommand::_execute(void) fixed_landmarks, moving_landmarks, initial_displacement, + affine_transform, constraint_mask, constraint_values, #ifdef DF_ENABLE_REGULARIZATION_WEIGHT_MAP diff --git a/src/deform_lib/registration/registration.cpp b/src/deform_lib/registration/registration.cpp index a48c13f..2f93848 100644 --- a/src/deform_lib/registration/registration.cpp +++ b/src/deform_lib/registration/registration.cpp @@ -91,6 +91,7 @@ stk::Volume registration( const std::vector& fixed_landmarks, const std::vector& moving_landmarks, const stk::Volume& initial_deformation, + const AffineTransform& affine_transform, const stk::Volume& constraint_mask, const stk::Volume& constraint_values, #ifdef DF_ENABLE_REGULARIZATION_WEIGHT_MAP @@ -195,6 +196,8 @@ stk::Volume registration( engine.set_initial_deformation(initial_deformation); } + engine.set_affine_transform(affine_transform); + // Constraints if (constraint_mask.valid() && constraint_values.valid()) { validate_volume_properties(constraint_mask, fixed_ref, "constraint mask"); diff --git a/src/deform_lib/registration/registration.h b/src/deform_lib/registration/registration.h index bd3c018..bad8c92 100644 --- a/src/deform_lib/registration/registration.h +++ b/src/deform_lib/registration/registration.h @@ -6,6 +6,8 @@ #include +class AffineTransform; + /*! * \brief Ensure that the input volume properties (size, origin, * spacing, direction) match the reference volume ones. @@ -34,6 +36,7 @@ stk::Volume registration( const std::vector& fixed_landmarks, const std::vector& moving_landmarks, const stk::Volume& initial_deformation, + const AffineTransform& affine_transform, const stk::Volume& constraint_mask, const stk::Volume& constraint_values, #ifdef DF_ENABLE_REGULARIZATION_WEIGHT_MAP diff --git a/src/deform_lib/registration/registration_engine.cpp b/src/deform_lib/registration/registration_engine.cpp index 7f9c8fd..2919801 100644 --- a/src/deform_lib/registration/registration_engine.cpp +++ b/src/deform_lib/registration/registration_engine.cpp @@ -538,6 +538,10 @@ void RegistrationEngine::set_initial_deformation(const stk::Volume& def) _deformation_pyramid.build_from_base(def, filters::downsample_vectorfield_by_2); #endif } +void RegistrationEngine::set_affine_transform(const AffineTransform& affine_transform) +{ + _affine_transform = affine_transform; +} void RegistrationEngine::set_image_pair( int i, const stk::Volume& fixed, diff --git a/src/deform_lib/registration/registration_engine.h b/src/deform_lib/registration/registration_engine.h index 8ee8303..b2201c8 100644 --- a/src/deform_lib/registration/registration_engine.h +++ b/src/deform_lib/registration/registration_engine.h @@ -18,6 +18,8 @@ class RegistrationEngine ~RegistrationEngine(); void set_initial_deformation(const stk::Volume& def); + void set_affine_transform(const AffineTransform& affine_transform); + void set_image_pair( int i, const stk::Volume& fixed, @@ -81,6 +83,7 @@ class RegistrationEngine VolumePyramid _regularization_weight_map; #endif // DF_ENABLE_REGULARIZATION_WEIGHT_MAP + AffineTransform _affine_transform; VolumePyramid _constraints_pyramid; VolumePyramid _constraints_mask_pyramid; From 751c3a62225300b6163178bbc3fc00d50f01b4ae Mon Sep 17 00:00:00 2001 From: Simon Date: Fri, 8 Nov 2019 15:56:53 +0100 Subject: [PATCH 27/37] fix --- src/deform_lib/registration/registration.cpp | 3 +++ src/deform_lib/registration/registration_engine.h | 4 +++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/deform_lib/registration/registration.cpp b/src/deform_lib/registration/registration.cpp index 2f93848..62b2082 100644 --- a/src/deform_lib/registration/registration.cpp +++ b/src/deform_lib/registration/registration.cpp @@ -247,6 +247,7 @@ stk::Volume registration( const std::vector& fixed_landmarks, const std::vector& moving_landmarks, const stk::Volume& initial_deformation, + const AffineTransform& affine_transform, const stk::Volume& constraint_mask, const stk::Volume& constraint_values, #ifdef DF_ENABLE_REGULARIZATION_WEIGHT_MAP @@ -270,6 +271,7 @@ stk::Volume registration( fixed_landmarks, moving_landmarks, initial_deformation, + affine_transform, constraint_mask, constraint_values, #ifdef DF_ENABLE_REGULARIZATION_WEIGHT_MAP @@ -288,6 +290,7 @@ stk::Volume registration( fixed_landmarks, moving_landmarks, initial_deformation, + affine_transform, constraint_mask, constraint_values, #ifdef DF_ENABLE_REGULARIZATION_WEIGHT_MAP diff --git a/src/deform_lib/registration/registration_engine.h b/src/deform_lib/registration/registration_engine.h index b2201c8..6a94d6d 100644 --- a/src/deform_lib/registration/registration_engine.h +++ b/src/deform_lib/registration/registration_engine.h @@ -1,6 +1,7 @@ #pragma once #include "../config.h" +#include "affine_transform.h" #include "settings.h" #include "volume_pyramid.h" @@ -23,7 +24,8 @@ class RegistrationEngine void set_image_pair( int i, const stk::Volume& fixed, - const stk::Volume& moving); + const stk::Volume& moving + ); #ifdef DF_ENABLE_REGULARIZATION_WEIGHT_MAP void set_regularization_weight_map(const stk::Volume& map); From 5b107ff364e5fc734cf47d805b89672016e9474c Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 11 Nov 2019 13:30:42 +0100 Subject: [PATCH 28/37] AffineTransform added to python interface --- src/python_wrapper/_pydeform.cpp | 60 ++++++++++++++++++++++++++++++-- 1 file changed, 57 insertions(+), 3 deletions(-) diff --git a/src/python_wrapper/_pydeform.cpp b/src/python_wrapper/_pydeform.cpp index b459f04..c689121 100644 --- a/src/python_wrapper/_pydeform.cpp +++ b/src/python_wrapper/_pydeform.cpp @@ -13,6 +13,8 @@ #include #include +#include +#include #include #include #include @@ -256,6 +258,9 @@ moving_landmarks: np.ndarray initial_displacement: stk.Volume Initial guess of the displacement field. +affine_transform: AffineTransform + Initial affine transformation + constraint_mask: stk.Volume Boolean mask for the constraints on the displacement. Requires to provide `constraint_values`. @@ -305,6 +310,7 @@ stk::Volume registration_wrapper( const py::object& fixed_landmarks, const py::object& moving_landmarks, const stk::Volume& initial_displacement, + const AffineTransform& affine_transform, const stk::Volume& constraint_mask, const stk::Volume& constraint_values, #ifdef DF_ENABLE_REGULARIZATION_WEIGHT_MAP @@ -412,6 +418,7 @@ stk::Volume registration_wrapper( fixed_landmarks_, moving_landmarks_, initial_displacement, + affine_transform, constraint_mask, constraint_values, #ifdef DF_ENABLE_REGULARIZATION_WEIGHT_MAP @@ -449,6 +456,9 @@ interpolator: pydeform.Interpolator `pydeform.Interpolator.Linear` or `pydeform.Interpolator.NearestNeighbour`. +affine_transform: AffineTransform + Optional affine transformation + Returns ------- stk.Volume @@ -459,9 +469,10 @@ stk.Volume stk::Volume transform_wrapper( const stk::Volume& src, const stk::Volume& df, - transform::Interp interp) + transform::Interp interp, + const AffineTransform& affine_transform) { - return transform_volume(src, df, interp); + return transform_volume(src, DisplacementField(df, affine_transform), interp); } std::string jacobian_docstring = @@ -528,6 +539,40 @@ stk::Volume regularization_wrapper( ); } + +AffineTransform make_affine_transform( + const py::array_t& matrix, + const py::array_t& offset +) +{ + AffineTransform t; + + if (matrix.ndim() != 2 || + matrix.shape(0) != 3 || + matrix.shape(1) != 3) { + throw py::value_error("Invalid shape of affine matrix, expected (3, 3)."); + } + + Matrix3x3f m; + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + m(i, j) = (float)matrix.at(i, j); + }} + + if (offset.ndim() != 1 || + offset.shape(0) != 3) { + throw py::value_error("Invalid shape of affine offset, expected (3)."); + } + + float3 o { + (float)offset.at(0), + (float)offset.at(1), + (float)offset.at(2) + }; + + return AffineTransform(m, o); +} + PYBIND11_MODULE(_pydeform, m) { m.attr("__version__") = GIT_VERSION_TAG; @@ -558,6 +603,13 @@ PYBIND11_MODULE(_pydeform, m) .value("Fatal", stk::LogLevel::Fatal, "Report only fatal errors") .export_values(); + py::class_(m, "AffineTransform") + .def(py::init<>()) + .def(py::init(&make_affine_transform), + py::arg("matrix"), + py::arg("offset") + ); + m.def("register", ®istration_wrapper, registration_docstring.c_str(), @@ -568,6 +620,7 @@ PYBIND11_MODULE(_pydeform, m) py::arg("fixed_landmarks") = py::none(), py::arg("moving_landmarks") = py::none(), py::arg("initial_displacement") = stk::Volume(), + py::arg("affine_transform") = AffineTransform(), py::arg("constraint_mask") = stk::Volume(), py::arg("constraint_values") = stk::Volume(), #ifdef DF_ENABLE_REGULARIZATION_WEIGHT_MAP @@ -586,7 +639,8 @@ PYBIND11_MODULE(_pydeform, m) transform_docstring.c_str(), py::arg("image"), py::arg("displacement"), - py::arg("interpolator") = transform::Interp_Linear + py::arg("interpolator") = transform::Interp_Linear, + py::arg("affine_transform") = AffineTransform() ); m.def("jacobian", From bb0e1cd522c49e360bd70d916edc41c353d65c76 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 11 Nov 2019 13:53:21 +0100 Subject: [PATCH 29/37] Apply affine transformation into resulting displacement field --- .../registration/displacement_field.cpp | 19 +++++++++++++++++++ .../registration/displacement_field.h | 7 +++++++ src/deform_lib/registration/registration.cpp | 4 ++-- .../registration/registration_engine.cpp | 6 +++++- 4 files changed, 33 insertions(+), 3 deletions(-) diff --git a/src/deform_lib/registration/displacement_field.cpp b/src/deform_lib/registration/displacement_field.cpp index 9ea5291..ed0b2d5 100644 --- a/src/deform_lib/registration/displacement_field.cpp +++ b/src/deform_lib/registration/displacement_field.cpp @@ -65,3 +65,22 @@ bool DisplacementField::valid() const return _df.valid(); } +stk::VolumeFloat3 compute_displacement_field( + const stk::VolumeFloat3& vector_field, + const AffineTransform& affine +) +{ + DisplacementField df(vector_field, affine); + + dim3 dims = df.size(); + stk::VolumeFloat3 out(dims); + + #pragma omp parallel for + for (int z = 0; z < (int)dims.z; ++z) { + for (int y = 0; y < (int)dims.y; ++y) { + for (int x = 0; x < (int)dims.x; ++x) { + out(x,y,z) = df.get(int3{x, y, z}); + }}} + + return out; +} diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index 01c72c7..c2ed3eb 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -101,3 +101,10 @@ class DisplacementField AffineTransform _affine; }; + +// Computes the composite of a vector field and an affine transform: +// u'(x) <- A(x + u(x)) + b - x +stk::VolumeFloat3 compute_displacement_field( + const stk::VolumeFloat3& vector_field, + const AffineTransform& affine +); diff --git a/src/deform_lib/registration/registration.cpp b/src/deform_lib/registration/registration.cpp index 62b2082..d8b08f5 100644 --- a/src/deform_lib/registration/registration.cpp +++ b/src/deform_lib/registration/registration.cpp @@ -227,7 +227,7 @@ stk::Volume registration( using namespace std::chrono; auto t_start = high_resolution_clock::now(); - stk::Volume def = engine.execute(); + stk::Volume df = engine.execute(); auto t_end = high_resolution_clock::now(); int elapsed = int(round(duration_cast>(t_end - t_start).count())); LOG(Info) << "Registration completed in " @@ -235,7 +235,7 @@ stk::Volume registration( << ":" << std::right << std::setw(2) << std::setfill('0') << elapsed % 60; - return def; + return df; } stk::Volume registration( diff --git a/src/deform_lib/registration/registration_engine.cpp b/src/deform_lib/registration/registration_engine.cpp index 2919801..75af8cb 100644 --- a/src/deform_lib/registration/registration_engine.cpp +++ b/src/deform_lib/registration/registration_engine.cpp @@ -705,7 +705,11 @@ stk::Volume RegistrationEngine::execute() } } - return _deformation_pyramid.volume(0); + // Composite affine transform and displacement field + return compute_displacement_field( + _deformation_pyramid.volume(0), + _affine_transform + ); } stk::Volume RegistrationEngine::deformation_field(int level) { From 52043b5c77a52a6a7e189d0dcd854d27b8cc17d8 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 11 Nov 2019 14:39:59 +0100 Subject: [PATCH 30/37] Copy meta data --- src/deform_lib/registration/displacement_field.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/deform_lib/registration/displacement_field.cpp b/src/deform_lib/registration/displacement_field.cpp index ed0b2d5..4d5d3d4 100644 --- a/src/deform_lib/registration/displacement_field.cpp +++ b/src/deform_lib/registration/displacement_field.cpp @@ -70,10 +70,11 @@ stk::VolumeFloat3 compute_displacement_field( const AffineTransform& affine ) { - DisplacementField df(vector_field, affine); - - dim3 dims = df.size(); + dim3 dims = vector_field.size(); stk::VolumeFloat3 out(dims); + out.copy_meta_from(vector_field); + + DisplacementField df(vector_field, affine); #pragma omp parallel for for (int z = 0; z < (int)dims.z; ++z) { From 64b6e7df2b6003d357316bbf418986a8824b5d81 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 11 Nov 2019 19:12:22 +0100 Subject: [PATCH 31/37] Error handling --- src/deform/registration.cpp | 8 +++++++- src/deform/transform.cpp | 11 +++++++++-- 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/src/deform/registration.cpp b/src/deform/registration.cpp index 32da9b0..9609131 100644 --- a/src/deform/registration.cpp +++ b/src/deform/registration.cpp @@ -164,7 +164,13 @@ int RegistrationCommand::_execute(void) AffineTransform affine_transform; if (!affine_file.empty()) { - affine_transform = parse_affine_transform_file(affine_file); + try{ + affine_transform = parse_affine_transform_file(affine_file); + } + catch (ValidationError& e) { + LOG(Error) << e.what(); + return EXIT_FAILURE; + } } // Constraints diff --git a/src/deform/transform.cpp b/src/deform/transform.cpp index 7153a46..acb251e 100644 --- a/src/deform/transform.cpp +++ b/src/deform/transform.cpp @@ -59,8 +59,15 @@ int TransformCommand::_execute(void) ASSERT(df.voxel_type() == stk::Type_Float3); AffineTransform affine; - if (_args.is_set("affine")) - affine = parse_affine_transform_file(_args.option("affine")); + if (_args.is_set("affine")) { + try{ + affine = parse_affine_transform_file(_args.option("affine")); + } + catch (ValidationError& e) { + LOG(Error) << e.what(); + return EXIT_FAILURE; + } + } stk::Volume result = transform_volume(src, DisplacementField(df, affine), interp); if (!result.valid()) From c0db00d171803228fb4f53019e14a06ee3bf9216 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 11 Nov 2019 20:30:41 +0100 Subject: [PATCH 32/37] Fix --- .../registration/blocked_graph_cut_optimizer.inl | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/deform_lib/registration/blocked_graph_cut_optimizer.inl b/src/deform_lib/registration/blocked_graph_cut_optimizer.inl index f42a308..b46945d 100644 --- a/src/deform_lib/registration/blocked_graph_cut_optimizer.inl +++ b/src/deform_lib/registration/blocked_graph_cut_optimizer.inl @@ -331,12 +331,7 @@ bool BlockedGraphCutOptimizer } if (solver.get_var(sub_x, sub_y, sub_z) == 1) { - if (update_rule == Settings::UpdateRule_Compositive) { - update_field.set(int3{gx, gy, gz}, delta); - } - else { - update_field.set(int3{gx, gy, gz}, delta); - } + update_field.set(int3{gx, gy, gz}, delta); changed_flag = true; } } From 275d944274d35aac498362a20a982d3c37c93ea4 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 11 Nov 2019 20:30:52 +0100 Subject: [PATCH 33/37] Actually append the affine transformation --- src/deform_lib/registration/registration_engine.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/deform_lib/registration/registration_engine.cpp b/src/deform_lib/registration/registration_engine.cpp index 75af8cb..71a2c26 100644 --- a/src/deform_lib/registration/registration_engine.cpp +++ b/src/deform_lib/registration/registration_engine.cpp @@ -606,7 +606,7 @@ stk::Volume RegistrationEngine::execute() #endif for (int l = _settings.num_pyramid_levels-1; l >= 0; --l) { - DisplacementField df(_deformation_pyramid.volume(l)); + DisplacementField df(_deformation_pyramid.volume(l), _affine_transform); std::vector neighborhood = determine_neighborhood(l); From a14cfdc333eb848d8a4fb52e291b0a00d97d483b Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 11 Nov 2019 20:58:05 +0100 Subject: [PATCH 34/37] Fixed composition with affine transform --- .../registration/displacement_field.cpp | 21 +++++++++++++++---- .../registration/displacement_field.h | 16 +++++++++----- 2 files changed, 28 insertions(+), 9 deletions(-) diff --git a/src/deform_lib/registration/displacement_field.cpp b/src/deform_lib/registration/displacement_field.cpp index 4d5d3d4..7b4bda7 100644 --- a/src/deform_lib/registration/displacement_field.cpp +++ b/src/deform_lib/registration/displacement_field.cpp @@ -29,13 +29,26 @@ void DisplacementField::update(const DisplacementField& update_field, bool compo for (int x = 0; x < (int)dims.x; ++x) { int3 p {x, y, z}; if (composite) { - float3 p1 = _df.index2point(p); - float3 p2 = p1 + update_field.get(p); - float3 p3 = buffer.transform_point(p2); + // Transform consists of two components, affine and displacement vector + // V(x) = Ax + b + // W(x) = x + u(x) + // + // The composition is computed as + // V(W(x)) = A(x + u(x)) + b + // + // Applying delta, produces following eq + // V(W'(x)) = A(x + u'(x)) + b = A(x + u(x + delta) + delta) + b + // + // Applying the delta only requires us to modify u(x) + // u'(x) = u(x + delta) + delta - _df(p) = p3 - p1; + float3 delta = update_field.get(p); + float3 p1 = _df.index2point(p); + _df(p) = _df.linear_at_point(p1 + delta, stk::Border_Replicate) + delta; } else { + // For additive updates we simply add the delta to the displacement + _df(p) += update_field.get(p); } }}} diff --git a/src/deform_lib/registration/displacement_field.h b/src/deform_lib/registration/displacement_field.h index c2ed3eb..dde02a3 100644 --- a/src/deform_lib/registration/displacement_field.h +++ b/src/deform_lib/registration/displacement_field.h @@ -41,14 +41,18 @@ class DisplacementField inline float3 get(const int3& p, const float3& delta, bool composite) const { if (composite) { + // T(x) = A(x + u(x + delta) + delta) + b float3 p1 = _df.index2point(p); - float3 p2 = p1 + delta; - float3 p3 = transform_point(p2); - - return p3 - p1; + float3 p2 = _affine.transform_point( + p1 + _df.linear_at_point(p1 + delta, stk::Border_Replicate) + delta + ); + return p2 - p1; } else { - return get(p) + delta; + // T(x) = A(x + u(x) + delta) + b + float3 p1 = _df.index2point(p); + float3 p2 = _affine.transform_point(p1 + _df(p) + delta); + return p2 - p1; } } @@ -64,6 +68,7 @@ class DisplacementField // Returns coordinates in world space inline float3 transform_point(const float3& p) const { + // T(x) = A(x + u(x)) + b return _affine.transform_point( p + _df.linear_at_point(p, stk::Border_Replicate) ); @@ -73,6 +78,7 @@ class DisplacementField // Returns coordinates in world space inline float3 transform_index(const int3& p) const { + // T(x) = A(x + u(x)) + b return _affine.transform_point(_df.index2point(p) + _df(p)); } From 71c07cf29376c8aed1fa5758830d0e110f9e614d Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 11 Nov 2019 22:11:25 +0100 Subject: [PATCH 35/37] WIP: Affine on GPU --- src/deform_lib/CMakeLists.txt | 1 + .../gpu/gpu_displacement_field.cpp | 49 +++++++++++++ .../registration/gpu/gpu_displacement_field.h | 68 ++++++------------- 3 files changed, 71 insertions(+), 47 deletions(-) create mode 100644 src/deform_lib/registration/gpu/gpu_displacement_field.cpp diff --git a/src/deform_lib/CMakeLists.txt b/src/deform_lib/CMakeLists.txt index 35ec827..e3af445 100644 --- a/src/deform_lib/CMakeLists.txt +++ b/src/deform_lib/CMakeLists.txt @@ -84,6 +84,7 @@ if (DF_USE_CUDA) "registration/gpu/cost_functions/cost_function.h" "registration/gpu/cost_functions/unary_function.h" + "registration/gpu/gpu_displacement_field.cpp" "registration/gpu/gpu_displacement_field.h" "registration/gpu/transform.cu" "registration/gpu/transform.h" diff --git a/src/deform_lib/registration/gpu/gpu_displacement_field.cpp b/src/deform_lib/registration/gpu/gpu_displacement_field.cpp new file mode 100644 index 0000000..375f115 --- /dev/null +++ b/src/deform_lib/registration/gpu/gpu_displacement_field.cpp @@ -0,0 +1,49 @@ +#include "gpu_displacement_field.h" + +GpuDisplacementField::GpuDisplacementField() +{ +} +GpuDisplacementField::GpuDisplacementField( + const stk::GpuVolume& df, + const AffineTransform& affine_transform +) : + _df(df), _affine_transform(affine_transform) +{ + ASSERT(_df.voxel_type() == stk::Type_Float4); +} +GpuDisplacementField::~GpuDisplacementField() +{ +} +GpuDisplacementField GpuDisplacementField::clone() +{ + return GpuDisplacementField(_df.clone(), _affine_transform); +} +void GpuDisplacementField::copy_from(const GpuDisplacementField& other) +{ + _df.copy_from(other._df); + _affine_transform = other._affine_transform; +} +dim3 GpuDisplacementField::size() const +{ + return _df.size(); +} +const stk::GpuVolume& GpuDisplacementField::volume() const +{ + return _df; +} +const AffineTransform& GpuDisplacementField::affine_transform() const +{ + return _affine_transform; +} +const float3& GpuDisplacementField::origin() const +{ + return _df.origin(); +} +const float3& GpuDisplacementField::spacing() const +{ + return _df.spacing(); +} +const Matrix3x3f& GpuDisplacementField::direction() const +{ + return _df.direction(); +} diff --git a/src/deform_lib/registration/gpu/gpu_displacement_field.h b/src/deform_lib/registration/gpu/gpu_displacement_field.h index 8feffc0..ecb480f 100644 --- a/src/deform_lib/registration/gpu/gpu_displacement_field.h +++ b/src/deform_lib/registration/gpu/gpu_displacement_field.h @@ -6,6 +6,7 @@ #include #include +#include "../affine_transform.h" #include "../settings.h" /** Wrapper for a displacement field image. @@ -16,63 +17,33 @@ class GpuDisplacementField { public: - GpuDisplacementField() : _update_rule(Settings::UpdateRule_Additive) {} + GpuDisplacementField(); GpuDisplacementField( const stk::GpuVolume& df, - Settings::UpdateRule update_rule = Settings::UpdateRule_Additive - ) : - _update_rule(update_rule), - _df(df) - { - ASSERT(_df.voxel_type() == stk::Type_Float4); - } - ~GpuDisplacementField() {} + const AffineTransform& affine_transform = AffineTransform() + ); + ~GpuDisplacementField(); - GpuDisplacementField clone() - { - return GpuDisplacementField(_df.clone(), _update_rule); - } + // Clones the displacement field and underlying data + GpuDisplacementField clone(); - void copy_from(const GpuDisplacementField& other) - { - _df.copy_from(other.volume()); - } + // Copies another displacement field into this + void copy_from(const GpuDisplacementField& other); - dim3 size() const - { - return _df.size(); - } + // Returns the dimensions of the displacement field + dim3 size() const; // Volume containing the displacements only - const stk::GpuVolume& volume() const - { - return _df; - } + const stk::GpuVolume& volume() const; + const AffineTransform& affine_transform() const; - const float3& origin() const - { - return _df.origin(); - } - - const float3& spacing() const - { - return _df.spacing(); - } - - const Matrix3x3f& direction() const - { - return _df.direction(); - } - - Settings::UpdateRule update_rule() const - { - return _update_rule; - } + const float3& origin() const; + const float3& spacing() const; + const Matrix3x3f& direction() const; private: - Settings::UpdateRule _update_rule; stk::GpuVolume _df; - + AffineTransform _affine_transform; }; #ifdef __CUDACC__ @@ -116,8 +87,10 @@ template class DisplacementField { public: - DisplacementField(const GpuDisplacementField& df) : + DisplacementField(const GpuDisplacementField& df, + const AffineTransform& affine_transform) : _df(df.volume()), + _affine_transform(affine_transform), _dims(df.size()), _origin(df.origin()), _spacing(df.spacing()), @@ -173,6 +146,7 @@ class DisplacementField private: stk::cuda::VolumePtr _df; + AffineTransform _affine_transform; dim3 _dims; float3 _origin; From b7abaffc3a51d91a11b0ea5c917493710171f94f Mon Sep 17 00:00:00 2001 From: Simon Date: Wed, 13 Nov 2019 22:43:58 +0100 Subject: [PATCH 36/37] GPU Implementation of affine initialization --- src/deform_lib/CMakeLists.txt | 1 + .../registration/affine_transform.h | 5 +- .../gpu/cost_functions/binary_function.cu | 5 +- .../gpu/cost_functions/binary_function.h | 1 + .../gpu/cost_functions/cost_function.h | 1 + .../gpu/cost_functions/cost_function_kernel.h | 12 ++- .../gpu/cost_functions/cross_correlation.cu | 11 ++- .../gpu/cost_functions/cross_correlation.h | 1 + .../gpu/cost_functions/landmarks.cu | 11 ++- .../gpu/cost_functions/landmarks.h | 1 + .../gpu/cost_functions/squared_distance.cu | 11 ++- .../gpu/cost_functions/squared_distance.h | 1 + .../gpu/cost_functions/unary_function.h | 12 ++- .../gpu/gpu_displacement_field.cu | 51 +++++++++++++ .../registration/gpu/gpu_displacement_field.h | 74 +++++++++++-------- src/deform_lib/registration/gpu/transform.cu | 2 +- .../registration/gpu_registration_engine.cpp | 21 +++++- .../registration/gpu_registration_engine.h | 7 +- .../hybrid_graph_cut_optimizer.cu | 44 ++++++----- .../registration/hybrid_graph_cut_optimizer.h | 3 + .../hybrid_graph_cut_optimizer.inl | 9 ++- src/deform_lib/registration/registration.cpp | 2 +- .../registration/registration_engine.cpp | 4 +- .../registration/registration_engine.h | 2 +- 24 files changed, 221 insertions(+), 71 deletions(-) create mode 100644 src/deform_lib/registration/gpu/gpu_displacement_field.cu diff --git a/src/deform_lib/CMakeLists.txt b/src/deform_lib/CMakeLists.txt index e3af445..86da1f3 100644 --- a/src/deform_lib/CMakeLists.txt +++ b/src/deform_lib/CMakeLists.txt @@ -85,6 +85,7 @@ if (DF_USE_CUDA) "registration/gpu/cost_functions/unary_function.h" "registration/gpu/gpu_displacement_field.cpp" + "registration/gpu/gpu_displacement_field.cu" "registration/gpu/gpu_displacement_field.h" "registration/gpu/transform.cu" "registration/gpu/transform.h" diff --git a/src/deform_lib/registration/affine_transform.h b/src/deform_lib/registration/affine_transform.h index 8e586d4..396ec45 100644 --- a/src/deform_lib/registration/affine_transform.h +++ b/src/deform_lib/registration/affine_transform.h @@ -5,6 +5,8 @@ #include #include +#include + class AffineTransform { public: @@ -12,8 +14,9 @@ class AffineTransform AffineTransform(const Matrix3x3f& matrix, const float3& offset); ~AffineTransform(); - inline float3 transform_point(const float3& pt) const + CUDA_HOST_DEVICE inline float3 transform_point(const float3& pt) const { + // Used by cuda::DisplacementField return _matrix * pt + _offset; } diff --git a/src/deform_lib/registration/gpu/cost_functions/binary_function.cu b/src/deform_lib/registration/gpu/cost_functions/binary_function.cu index 9e6931e..75fc78a 100644 --- a/src/deform_lib/registration/gpu/cost_functions/binary_function.cu +++ b/src/deform_lib/registration/gpu/cost_functions/binary_function.cu @@ -186,6 +186,7 @@ void GpuBinaryFunction::operator()( const float3& delta, const int3& offset, const int3& dims, + Settings::UpdateRule update_rule, stk::GpuVolume& cost_x, stk::GpuVolume& cost_y, stk::GpuVolume& cost_z, @@ -222,7 +223,7 @@ void GpuBinaryFunction::operator()( 0 }; - if (df.update_rule() == Settings::UpdateRule_Compositive) { + if (update_rule == Settings::UpdateRule_Compositive) { regularizer_kernel> <<>>( df, @@ -238,7 +239,7 @@ void GpuBinaryFunction::operator()( cost_z ); } - else if (df.update_rule() == Settings::UpdateRule_Additive) { + else if (update_rule == Settings::UpdateRule_Additive) { regularizer_kernel> <<>>( df, diff --git a/src/deform_lib/registration/gpu/cost_functions/binary_function.h b/src/deform_lib/registration/gpu/cost_functions/binary_function.h index c0ffcea..c56f311 100644 --- a/src/deform_lib/registration/gpu/cost_functions/binary_function.h +++ b/src/deform_lib/registration/gpu/cost_functions/binary_function.h @@ -50,6 +50,7 @@ class GpuBinaryFunction const float3& delta, const int3& offset, const int3& dims, + Settings::UpdateRule update_rule, stk::GpuVolume& cost_x, stk::GpuVolume& cost_y, stk::GpuVolume& cost_z, diff --git a/src/deform_lib/registration/gpu/cost_functions/cost_function.h b/src/deform_lib/registration/gpu/cost_functions/cost_function.h index cc0c20b..52c8f13 100644 --- a/src/deform_lib/registration/gpu/cost_functions/cost_function.h +++ b/src/deform_lib/registration/gpu/cost_functions/cost_function.h @@ -32,6 +32,7 @@ struct GpuCostFunction float weight, const int3& offset, const int3& dims, + Settings::UpdateRule update_rule, stk::GpuVolume& cost_acc, stk::cuda::Stream& stream ) = 0; diff --git a/src/deform_lib/registration/gpu/cost_functions/cost_function_kernel.h b/src/deform_lib/registration/gpu/cost_functions/cost_function_kernel.h index c7221cc..ab4c582 100644 --- a/src/deform_lib/registration/gpu/cost_functions/cost_function_kernel.h +++ b/src/deform_lib/registration/gpu/cost_functions/cost_function_kernel.h @@ -72,8 +72,12 @@ struct CostFunctionKernel } } + int3 p {x,y,z}; + float4 d = df.get(p, delta); + float3 pt = df.index2point(p) + float3{d.x, d.y, d.z}; + const float3 moving_p = _inv_moving_direction - * (df.transform_index({x,y,z}, delta) - _moving_origin) + * (pt - _moving_origin) * _inv_moving_spacing; // Check if the moving voxels are masked out @@ -87,7 +91,6 @@ struct CostFunctionKernel } } - float4 d = df.get(int3{x, y, z}, delta); float c = _impl( _fixed, _moving, @@ -155,6 +158,7 @@ void invoke_cost_function_kernel( const int3& offset, const int3& dims, GpuDisplacementField& df, + Settings::UpdateRule update_rule, stk::cuda::Stream& stream ) { @@ -183,7 +187,7 @@ void invoke_cost_function_kernel( float4 d4 { delta.x, delta.y, delta.z, 0 }; - if (df.update_rule() == Settings::UpdateRule_Compositive) { + if (update_rule == Settings::UpdateRule_Compositive) { cost_function_kernel<<>>( kernel, offset, @@ -193,7 +197,7 @@ void invoke_cost_function_kernel( 1 ); } - else if (df.update_rule() == Settings::UpdateRule_Additive) { + else if (update_rule == Settings::UpdateRule_Additive) { cost_function_kernel<<>>( kernel, offset, diff --git a/src/deform_lib/registration/gpu/cost_functions/cross_correlation.cu b/src/deform_lib/registration/gpu/cost_functions/cross_correlation.cu index 4e151e3..355ec19 100644 --- a/src/deform_lib/registration/gpu/cost_functions/cross_correlation.cu +++ b/src/deform_lib/registration/gpu/cost_functions/cross_correlation.cu @@ -95,6 +95,7 @@ void GpuCostFunction_NCC::cost( float weight, const int3& offset, const int3& dims, + Settings::UpdateRule update_rule, stk::GpuVolume& cost_acc, stk::cuda::Stream& stream ) @@ -117,6 +118,14 @@ void GpuCostFunction_NCC::cost( cost_acc ); - invoke_cost_function_kernel(kernel, delta, offset, dims, df, stream); + invoke_cost_function_kernel( + kernel, + delta, + offset, + dims, + df, + update_rule, + stream + ); } diff --git a/src/deform_lib/registration/gpu/cost_functions/cross_correlation.h b/src/deform_lib/registration/gpu/cost_functions/cross_correlation.h index 2b48068..0d6afe1 100644 --- a/src/deform_lib/registration/gpu/cost_functions/cross_correlation.h +++ b/src/deform_lib/registration/gpu/cost_functions/cross_correlation.h @@ -22,6 +22,7 @@ struct GpuCostFunction_NCC : public GpuCostFunction float weight, const int3& offset, const int3& dims, + Settings::UpdateRule update_rule, stk::GpuVolume& cost_acc, stk::cuda::Stream& stream ); diff --git a/src/deform_lib/registration/gpu/cost_functions/landmarks.cu b/src/deform_lib/registration/gpu/cost_functions/landmarks.cu index 31622c3..d9b34dc 100644 --- a/src/deform_lib/registration/gpu/cost_functions/landmarks.cu +++ b/src/deform_lib/registration/gpu/cost_functions/landmarks.cu @@ -104,6 +104,7 @@ void GpuCostFunction_Landmarks::cost( float weight, const int3& offset, const int3& dims, + Settings::UpdateRule update_rule, stk::GpuVolume& cost_acc, stk::cuda::Stream& stream ) @@ -128,7 +129,15 @@ void GpuCostFunction_Landmarks::cost( cost_acc ); - invoke_cost_function_kernel(kernel, delta, offset, dims, df, stream); + invoke_cost_function_kernel( + kernel, + delta, + offset, + dims, + df, + update_rule, + stream + ); } diff --git a/src/deform_lib/registration/gpu/cost_functions/landmarks.h b/src/deform_lib/registration/gpu/cost_functions/landmarks.h index eb4ab0b..204be00 100644 --- a/src/deform_lib/registration/gpu/cost_functions/landmarks.h +++ b/src/deform_lib/registration/gpu/cost_functions/landmarks.h @@ -30,6 +30,7 @@ struct GpuCostFunction_Landmarks : public GpuCostFunction float weight, const int3& offset, const int3& dims, + Settings::UpdateRule update_rule, stk::GpuVolume& cost_acc, stk::cuda::Stream& stream ); diff --git a/src/deform_lib/registration/gpu/cost_functions/squared_distance.cu b/src/deform_lib/registration/gpu/cost_functions/squared_distance.cu index faf7cd2..86e4cc4 100644 --- a/src/deform_lib/registration/gpu/cost_functions/squared_distance.cu +++ b/src/deform_lib/registration/gpu/cost_functions/squared_distance.cu @@ -36,6 +36,7 @@ void GpuCostFunction_SSD::cost( float weight, const int3& offset, const int3& dims, + Settings::UpdateRule update_rule, stk::GpuVolume& cost_acc, stk::cuda::Stream& stream ) @@ -58,6 +59,14 @@ void GpuCostFunction_SSD::cost( cost_acc ); - invoke_cost_function_kernel(kernel, delta, offset, dims, df, stream); + invoke_cost_function_kernel( + kernel, + delta, + offset, + dims, + df, + update_rule, + stream + ); } diff --git a/src/deform_lib/registration/gpu/cost_functions/squared_distance.h b/src/deform_lib/registration/gpu/cost_functions/squared_distance.h index 1ffefe6..ace45dc 100644 --- a/src/deform_lib/registration/gpu/cost_functions/squared_distance.h +++ b/src/deform_lib/registration/gpu/cost_functions/squared_distance.h @@ -19,6 +19,7 @@ struct GpuCostFunction_SSD : public GpuCostFunction float weight, const int3& offset, const int3& dims, + Settings::UpdateRule update_rule, stk::GpuVolume& cost_acc, stk::cuda::Stream& stream ); diff --git a/src/deform_lib/registration/gpu/cost_functions/unary_function.h b/src/deform_lib/registration/gpu/cost_functions/unary_function.h index 550f29c..1aba000 100644 --- a/src/deform_lib/registration/gpu/cost_functions/unary_function.h +++ b/src/deform_lib/registration/gpu/cost_functions/unary_function.h @@ -26,12 +26,22 @@ class GpuUnaryFunction const float3& delta, const int3& offset, const int3& dims, + Settings::UpdateRule update_rule, stk::GpuVolume& cost_acc, stk::cuda::Stream& stream ) { for (auto& fn : _functions) { - fn.function->cost(df, delta, fn.weight, offset, dims, cost_acc, stream); + fn.function->cost( + df, + delta, + fn.weight, + offset, + dims, + update_rule, + cost_acc, + stream + ); } } diff --git a/src/deform_lib/registration/gpu/gpu_displacement_field.cu b/src/deform_lib/registration/gpu/gpu_displacement_field.cu new file mode 100644 index 0000000..3224154 --- /dev/null +++ b/src/deform_lib/registration/gpu/gpu_displacement_field.cu @@ -0,0 +1,51 @@ +#include "gpu_displacement_field.h" + +#include +#include + +namespace cuda { + using namespace stk::cuda; +} + +__global__ void compute_displacement_field_kernel( + cuda::DisplacementField<> df_in, + cuda::VolumePtr df_out +) +{ + int x = blockIdx.x*blockDim.x + threadIdx.x; + int y = blockIdx.y*blockDim.y + threadIdx.y; + int z = blockIdx.z*blockDim.z + threadIdx.z; + + if (x >= df_in.size().x || + y >= df_in.size().y || + z >= df_in.size().z) + { + return; + } + df_out(x, y, z) = df_in.get(int3{x, y, z}); +} + +stk::GpuVolume cuda::compute_displacement_field( + const stk::GpuVolume& vector_field, + const AffineTransform& affine +) +{ + ASSERT(vector_field.voxel_type() == stk::Type_Float4); + + stk::GpuVolume out(vector_field.size(), stk::Type_Float4); + + dim3 dims = vector_field.size(); + dim3 block_size {32,32,1}; + dim3 grid_size { + (dims.x + block_size.x - 1) / block_size.x, + (dims.y + block_size.y - 1) / block_size.y, + (dims.z + block_size.z - 1) / block_size.z + }; + GpuDisplacementField df(vector_field, affine); + compute_displacement_field_kernel<<>>( + cuda::DisplacementField<>(df), + out + ); + return out; +} + diff --git a/src/deform_lib/registration/gpu/gpu_displacement_field.h b/src/deform_lib/registration/gpu/gpu_displacement_field.h index ecb480f..b7e43ce 100644 --- a/src/deform_lib/registration/gpu/gpu_displacement_field.h +++ b/src/deform_lib/registration/gpu/gpu_displacement_field.h @@ -46,51 +46,53 @@ class GpuDisplacementField AffineTransform _affine_transform; }; -#ifdef __CUDACC__ - namespace cuda { +#ifdef __CUDACC__ + struct CompositiveUpdate { - __device__ float4 operator()( + __device__ float3 operator()( const stk::cuda::VolumePtr& df, const dim3& dims, const float3& inv_spacing, - int x, int y, int z, + const int3& index, const float4& delta ) { // Convert delta from mm to image space - return stk::cuda::linear_at_clamp( + float4 o = stk::cuda::linear_at_clamp( df, dims, - x + delta.x * inv_spacing.x, - y + delta.y * inv_spacing.y, - z + delta.z * inv_spacing.z + index.x + delta.x * inv_spacing.x, + index.y + delta.y * inv_spacing.y, + index.z + delta.z * inv_spacing.z ) + delta; + return float3{o.x, o.y, o.z}; } }; struct AdditiveUpdate { - __device__ float4 operator()( + __device__ float3 operator()( const stk::cuda::VolumePtr& df, const dim3& /*dims*/, const float3& /*inv_spacing*/, - int x, int y, int z, + const int3& index, const float4& delta ) { - return df(x, y, z) + delta; + float4 o = df(index.x, index.y, index.z) + delta; + return float3{o.x, o.y, o.z}; } }; + template class DisplacementField { public: - DisplacementField(const GpuDisplacementField& df, - const AffineTransform& affine_transform) : + DisplacementField(const GpuDisplacementField& df) : _df(df.volume()), - _affine_transform(affine_transform), + _affine(df.affine_transform()), _dims(df.size()), _origin(df.origin()), _spacing(df.spacing()), @@ -106,31 +108,30 @@ class DisplacementField __device__ float4 get(const int3& p) const { - return _df(p.x, p.y, p.z); + float3 out = transform_index(p) - index2point(p); + return float4{out.x, out.y, out.z, 0}; } __device__ float4 get(const int3& p, const float4& delta) const { + float3 p1 = index2point(p); + TUpdateFn fn; - return fn(_df, _dims, _inv_spacing, p.x, p.y, p.z, delta); - } + float3 p2 = p1 + fn(_df, _dims, _inv_spacing, p, delta); + float3 p3 = _affine.transform_point(p2); - // p : Index in displacement field - // Returns coordinates in world space - __device__ float3 transform_index(const int3& p) const - { - float4 d = get(p); - return index2point(p) + float3{d.x, d.y, d.z}; + float3 out = p3 - p1; + return float4{out.x, out.y, out.z, 0}; } - // Transforms the index with a given delta // p : Index in displacement field - // delta : Delta to apply to the transformation // Returns coordinates in world space - __device__ float3 transform_index(const int3& p, const float4& delta) const + __device__ float3 transform_index(const int3& p) const { - float4 d = get(p, delta); - return index2point(p) + float3{d.x, d.y, d.z}; + float4 d = _df(p.x, p.y, p.z); + return _affine.transform_point( + index2point(p) + float3{d.x, d.y, d.z} + ); } __device__ float3 index2point(const int3& index) const @@ -138,7 +139,7 @@ class DisplacementField float3 xyz {float(index.x), float(index.y), float(index.z)}; return _origin + _direction * (_spacing * xyz); } - + __device__ const dim3& size() const { return _dims; @@ -146,13 +147,24 @@ class DisplacementField private: stk::cuda::VolumePtr _df; - AffineTransform _affine_transform; + ::AffineTransform _affine; dim3 _dims; float3 _origin; float3 _spacing; float3 _inv_spacing; Matrix3x3f _direction; + Matrix3x3f _inv_direction; }; -} // namespace cuda + #endif // __CUDACC__ + + +// Computes the composite of a vector field and an affine transform: +// u'(x) <- A(x + u(x)) + b - x +stk::GpuVolume compute_displacement_field( + const stk::GpuVolume& vector_field, + const AffineTransform& affine +); + +} // namespace cuda diff --git a/src/deform_lib/registration/gpu/transform.cu b/src/deform_lib/registration/gpu/transform.cu index ffadc1d..4b24a77 100644 --- a/src/deform_lib/registration/gpu/transform.cu +++ b/src/deform_lib/registration/gpu/transform.cu @@ -95,7 +95,7 @@ static void run_nn_kernel( transform_kernel_nn<<>>( \ src, \ src.size(), \ - df, \ + cuda::DisplacementField<>(df), \ src.origin(), \ inv_moving_spacing, \ src.inverse_direction(), \ diff --git a/src/deform_lib/registration/gpu_registration_engine.cpp b/src/deform_lib/registration/gpu_registration_engine.cpp index c388a00..40ea33b 100644 --- a/src/deform_lib/registration/gpu_registration_engine.cpp +++ b/src/deform_lib/registration/gpu_registration_engine.cpp @@ -15,6 +15,7 @@ #include "gpu/cost_functions/landmarks.h" #include "gpu/cost_functions/squared_distance.h" #include "gpu/cost_functions/unary_function.h" +#include "gpu/gpu_displacement_field.h" #include "hybrid_graph_cut_optimizer.h" #include "../filters/gpu/resample.h" @@ -173,7 +174,7 @@ GpuRegistrationEngine::~GpuRegistrationEngine() { } -void GpuRegistrationEngine::set_initial_deformation(const stk::Volume& def) +void GpuRegistrationEngine::set_initial_displacement_field(const stk::Volume& def) { // GPU prefers float4 over float3 ASSERT(def.voxel_type() == stk::Type_Float4 || def.voxel_type() == stk::Type_Float3); @@ -188,6 +189,10 @@ void GpuRegistrationEngine::set_initial_deformation(const stk::Volume& def) _deformation_pyramid.build_from_base(gpu_def, filters::gpu::downsample_vectorfield_by_2); } +void GpuRegistrationEngine::set_affine_transform(const AffineTransform& affine_transform) +{ + _affine_transform = affine_transform; +} void GpuRegistrationEngine::set_image_pair( int i, const stk::Volume& fixed, @@ -272,11 +277,11 @@ stk::Volume GpuRegistrationEngine::execute() initial.set_spacing(base.spacing()); initial.set_direction(base.direction()); - set_initial_deformation(initial); + set_initial_displacement_field(initial); } for (int l = _settings.num_pyramid_levels-1; l >= 0; --l) { - GpuDisplacementField df(_deformation_pyramid.volume(l), _settings.update_rule); + GpuDisplacementField df(_deformation_pyramid.volume(l), _affine_transform); std::vector neighborhood = determine_neighborhood(l); @@ -293,6 +298,7 @@ stk::Volume GpuRegistrationEngine::execute() HybridGraphCutOptimizer> optimizer( neighborhood, _settings.levels[l], + _settings.update_rule, unary_fn, binary_fn, df, @@ -306,6 +312,7 @@ stk::Volume GpuRegistrationEngine::execute() HybridGraphCutOptimizer> optimizer( neighborhood, _settings.levels[l], + _settings.update_rule, unary_fn, binary_fn, df, @@ -320,6 +327,7 @@ stk::Volume GpuRegistrationEngine::execute() HybridGraphCutOptimizer> optimizer( neighborhood, _settings.levels[l], + _settings.update_rule unary_fn, binary_fn, df, @@ -347,7 +355,12 @@ stk::Volume GpuRegistrationEngine::execute() } } - return volume_float4_to_float3(_deformation_pyramid.volume(0).download()); + return volume_float4_to_float3( + cuda::compute_displacement_field( + _deformation_pyramid.volume(0), + _affine_transform + ).download() + ); } std::vector GpuRegistrationEngine::determine_neighborhood(int level) const { diff --git a/src/deform_lib/registration/gpu_registration_engine.h b/src/deform_lib/registration/gpu_registration_engine.h index cc13876..b18f621 100644 --- a/src/deform_lib/registration/gpu_registration_engine.h +++ b/src/deform_lib/registration/gpu_registration_engine.h @@ -1,6 +1,7 @@ #pragma once #include "../config.h" +#include "affine_transform.h" #include "settings.h" #include "gpu_volume_pyramid.h" #include "worker_pool.h" @@ -19,7 +20,9 @@ class GpuRegistrationEngine GpuRegistrationEngine(const Settings& settings); ~GpuRegistrationEngine(); - void set_initial_deformation(const stk::Volume& def); + void set_initial_displacement_field(const stk::Volume& def); + void set_affine_transform(const AffineTransform& affine_transform); + void set_image_pair( int i, const stk::Volume& fixed, @@ -69,6 +72,8 @@ class GpuRegistrationEngine std::vector _fixed_landmarks; std::vector _moving_landmarks; + AffineTransform _affine_transform; + WorkerPool _worker_pool; std::vector _stream_pool; }; diff --git a/src/deform_lib/registration/hybrid_graph_cut_optimizer.cu b/src/deform_lib/registration/hybrid_graph_cut_optimizer.cu index 40df7a4..c56b896 100644 --- a/src/deform_lib/registration/hybrid_graph_cut_optimizer.cu +++ b/src/deform_lib/registration/hybrid_graph_cut_optimizer.cu @@ -14,25 +14,29 @@ namespace cuda { template __global__ void apply_displacement_delta_kernel( - cuda::DisplacementField df_in, - cuda::VolumePtr df_out, + cuda::VolumePtr df_in, + dim3 dims, cuda::VolumePtr labels, - float4 delta + float4 delta, + float3 inv_spacing, + cuda::VolumePtr df_out ) { int x = blockIdx.x*blockDim.x + threadIdx.x; int y = blockIdx.y*blockDim.y + threadIdx.y; int z = blockIdx.z*blockDim.z + threadIdx.z; - if (x >= df_in.size().x || - y >= df_in.size().y || - z >= df_in.size().z) + if (x >= dims.x || + y >= dims.y || + z >= dims.z) { return; } if (labels(x,y,z)) { - df_out(x,y,z) = df_in.get(int3{x, y, z}, delta); + UpdateFn fn; + float3 d = fn(df_in, dims, inv_spacing, int3{x,y,z}, delta); + df_out(x,y,z) = {d.x, d.y, d.z}; } } @@ -41,6 +45,7 @@ void apply_displacement_delta( GpuDisplacementField& df_out, stk::GpuVolume& labels, const float3& delta, + Settings::UpdateRule update_rule, cuda::Stream stream ) { @@ -58,24 +63,28 @@ void apply_displacement_delta( 1.0f / df_in.spacing().z }; - if (df_in.update_rule() == Settings::UpdateRule_Additive) { + if (update_rule == Settings::UpdateRule_Additive) { // In and out buffer for displacement field in the additive case can // be the same, since all updates are guaranteed to be independent. - apply_displacement_delta_kernel + apply_displacement_delta_kernel <<>>( - cuda::DisplacementField(df_out), df_out.volume(), + dims, labels, - float4{delta.x, delta.y, delta.z, 0.0f} + float4{delta.x, delta.y, delta.z, 0.0f}, + inv_spacing, + df_out.volume() ); } - else if (df_in.update_rule() == Settings::UpdateRule_Compositive) { - apply_displacement_delta_kernel + else if (update_rule == Settings::UpdateRule_Compositive) { + apply_displacement_delta_kernel <<>>( - cuda::DisplacementField(df_in), - df_out.volume(), + df_in.volume(), + dims, labels, - float4{delta.x, delta.y, delta.z, 0.0f} + float4{delta.x, delta.y, delta.z, 0.0f}, + inv_spacing, + df_out.volume() ); } CUDA_CHECK_ERRORS(cudaPeekAtLastError()); @@ -166,7 +175,7 @@ double calculate_energy( cuda::Stream& stream = stk::cuda::Stream::null(); // Update rule doesn't matter in this case since we don't want the energy for a move. - unary_fn(df, {0,0,0}, begin, end, unary_cost, stream); + unary_fn(df, {0,0,0}, begin, end, Settings::UpdateRule_Additive, unary_cost, stream); // Compute binary terms binary_fn( @@ -174,6 +183,7 @@ double calculate_energy( {0, 0, 0}, begin, end, + Settings::UpdateRule_Additive, binary_cost_x, binary_cost_y, binary_cost_z, diff --git a/src/deform_lib/registration/hybrid_graph_cut_optimizer.h b/src/deform_lib/registration/hybrid_graph_cut_optimizer.h index 1ef7048..ce46415 100644 --- a/src/deform_lib/registration/hybrid_graph_cut_optimizer.h +++ b/src/deform_lib/registration/hybrid_graph_cut_optimizer.h @@ -28,6 +28,7 @@ class HybridGraphCutOptimizer HybridGraphCutOptimizer( const std::vector& neighborhood, const Settings::Level& settings, + Settings::UpdateRule update_rule, GpuUnaryFunction& unary_fn, GpuBinaryFunction& binary_fn, GpuDisplacementField& df, @@ -92,6 +93,7 @@ class HybridGraphCutOptimizer std::vector _neighborhood; const Settings::Level& _settings; + Settings::UpdateRule _update_rule; WorkerPool& _worker_pool; std::vector& _stream_pool; @@ -141,6 +143,7 @@ void apply_displacement_delta( GpuDisplacementField& df_out, stk::GpuVolume& labels, const float3& delta, + Settings::UpdateRule update_rule, stk::cuda::Stream stream ); diff --git a/src/deform_lib/registration/hybrid_graph_cut_optimizer.inl b/src/deform_lib/registration/hybrid_graph_cut_optimizer.inl index 6bcae2e..d449783 100644 --- a/src/deform_lib/registration/hybrid_graph_cut_optimizer.inl +++ b/src/deform_lib/registration/hybrid_graph_cut_optimizer.inl @@ -15,6 +15,7 @@ template HybridGraphCutOptimizer::HybridGraphCutOptimizer( const std::vector& neighborhood, const Settings::Level& settings, + Settings::UpdateRule update_rule, GpuUnaryFunction& unary_fn, GpuBinaryFunction& binary_fn, GpuDisplacementField& df, @@ -22,6 +23,7 @@ HybridGraphCutOptimizer::HybridGraphCutOptimizer( std::vector& stream_pool) : _neighborhood(neighborhood), _settings(settings), + _update_rule(update_rule), _worker_pool(worker_pool), _stream_pool(stream_pool), _unary_fn(unary_fn), @@ -34,7 +36,7 @@ HybridGraphCutOptimizer::HybridGraphCutOptimizer( // for the previous iteration since each update may depend on neighbouring // updates. As compared to additive were we only have d(x) = d(x) + delta. - if (_df.update_rule() == Settings::UpdateRule_Compositive) { + if (_update_rule == Settings::UpdateRule_Compositive) { _df_tmp = _df.clone(); } @@ -287,7 +289,7 @@ size_t HybridGraphCutOptimizer::dispatch_blocks() // updates. As compared to additive were we only have d(x) = d(x) + delta. GpuDisplacementField df_in = _df; - if (_df.update_rule() == Settings::UpdateRule_Compositive) { + if (_update_rule == Settings::UpdateRule_Compositive) { _df_tmp.copy_from(_df); df_in = _df_tmp; } @@ -297,6 +299,7 @@ size_t HybridGraphCutOptimizer::dispatch_blocks() _df, _gpu_labels, _current_delta, + _update_rule, stk::cuda::Stream::null() ); } @@ -400,6 +403,7 @@ void HybridGraphCutOptimizer::block_cost_task( _current_delta, block.begin, block_dims, + _update_rule, _gpu_unary_cost, stream ); @@ -419,6 +423,7 @@ void HybridGraphCutOptimizer::block_cost_task( _current_delta, block.begin, block_dims, + _update_rule, _gpu_binary_cost_x, _gpu_binary_cost_y, _gpu_binary_cost_z, diff --git a/src/deform_lib/registration/registration.cpp b/src/deform_lib/registration/registration.cpp index d8b08f5..566cbcd 100644 --- a/src/deform_lib/registration/registration.cpp +++ b/src/deform_lib/registration/registration.cpp @@ -193,7 +193,7 @@ stk::Volume registration( // Initial deformation if (initial_deformation.valid()) { validate_volume_properties(initial_deformation, fixed_ref, "initial deformation field"); - engine.set_initial_deformation(initial_deformation); + engine.set_initial_displacement_field(initial_deformation); } engine.set_affine_transform(affine_transform); diff --git a/src/deform_lib/registration/registration_engine.cpp b/src/deform_lib/registration/registration_engine.cpp index 71a2c26..2f454fa 100644 --- a/src/deform_lib/registration/registration_engine.cpp +++ b/src/deform_lib/registration/registration_engine.cpp @@ -527,7 +527,7 @@ void RegistrationEngine::set_moving_mask(const stk::VolumeFloat& moving_mask) _moving_mask_pyramid.set_level_count(_settings.num_pyramid_levels); _moving_mask_pyramid.build_from_base(moving_mask, filters::downsample_volume_by_2); } -void RegistrationEngine::set_initial_deformation(const stk::Volume& def) +void RegistrationEngine::set_initial_displacement_field(const stk::Volume& def) { ASSERT(def.voxel_type() == stk::Type_Float3); // Only single-precision supported for now ASSERT(_settings.num_pyramid_levels); @@ -598,7 +598,7 @@ stk::Volume RegistrationEngine::execute() stk::VolumeFloat3 initial(base.size(), float3{0, 0, 0}); initial.copy_meta_from(base); - set_initial_deformation(initial); + set_initial_displacement_field(initial); } #ifdef DF_OUTPUT_DEBUG_VOLUMES diff --git a/src/deform_lib/registration/registration_engine.h b/src/deform_lib/registration/registration_engine.h index 6a94d6d..a201843 100644 --- a/src/deform_lib/registration/registration_engine.h +++ b/src/deform_lib/registration/registration_engine.h @@ -18,7 +18,7 @@ class RegistrationEngine RegistrationEngine(const Settings& settings); ~RegistrationEngine(); - void set_initial_deformation(const stk::Volume& def); + void set_initial_displacement_field(const stk::Volume& def); void set_affine_transform(const AffineTransform& affine_transform); void set_image_pair( From 347c5a4e9a6817d29f2bc34e25fc286768a50f73 Mon Sep 17 00:00:00 2001 From: Simon Date: Thu, 14 Nov 2019 09:24:33 +0100 Subject: [PATCH 37/37] Copy meta data --- src/deform_lib/registration/gpu/gpu_displacement_field.cu | 1 + 1 file changed, 1 insertion(+) diff --git a/src/deform_lib/registration/gpu/gpu_displacement_field.cu b/src/deform_lib/registration/gpu/gpu_displacement_field.cu index 3224154..1b708eb 100644 --- a/src/deform_lib/registration/gpu/gpu_displacement_field.cu +++ b/src/deform_lib/registration/gpu/gpu_displacement_field.cu @@ -33,6 +33,7 @@ stk::GpuVolume cuda::compute_displacement_field( ASSERT(vector_field.voxel_type() == stk::Type_Float4); stk::GpuVolume out(vector_field.size(), stk::Type_Float4); + out.copy_meta_from(vector_field); dim3 dims = vector_field.size(); dim3 block_size {32,32,1};