From 8e0ee6a864ec106dbaaa383ea4f72317a104df3f Mon Sep 17 00:00:00 2001 From: mattcieslak Date: Sun, 18 Jan 2026 11:58:25 -0500 Subject: [PATCH 01/11] add chen mode --- auto_track.cpp | 152 ++++++++++++++++++++++++++++++++++++------------- 1 file changed, 113 insertions(+), 39 deletions(-) diff --git a/auto_track.cpp b/auto_track.cpp index 605f852b..8d1c0e83 100644 --- a/auto_track.cpp +++ b/auto_track.cpp @@ -1,6 +1,8 @@ #include #include #include +#include +#include #include "auto_track.h" #include "ui_auto_track.h" #include "libs/dsi/image_model.hpp" @@ -8,6 +10,22 @@ #include "libs/tracking/tracking_thread.hpp" #include extern std::vector fa_template_list; +static bool find_string_case_insensitive(const std::string & str1, const std::string & str2) +{ + auto it = std::search( + str1.begin(), str1.end(), + str2.begin(), str2.end(), + [](char ch1, char ch2) { return std::toupper(ch1) == std::toupper(ch2); } + ); + return (it != str1.end() ); +} +static bool is_selected(std::vector& selected_tracts,const std::string& tract_name) +{ + for(const auto& each: selected_tracts) + if(find_string_case_insensitive(tract_name,each)) + return true; + return false; +} auto_track::auto_track(QWidget *parent) : QMainWindow(parent), ui(new Ui::auto_track) @@ -105,12 +123,23 @@ int trk_post(tipl::program_option& po, std::string tract_file_name,bool output_track); std::string run_auto_track(tipl::program_option& po,const std::vector& file_list,int& prog) { + const bool chen_mode = po.get("chen_mode",0); std::string trk_format = po.get("trk_format","tt.gz"); std::string tolerance_string = po.get("tolerance","22,26,30"); float yield_rate = po.get("yield_rate",0.00001f); size_t yield_check_count = 10.0f/yield_rate; bool overwrite = po.get("overwrite",0); uint32_t thread_count = tipl::max_thread_count; + if(chen_mode) + { + po.set("template",po.get("template",0)); + po.set("track_voxel_ratio",po.get("track_voxel_ratio",2.0f)); + po.set("tip_iteration",po.get("tip_iteration",4)); + po.set("use_roi",po.get("use_roi",0)); + po.set("check_ending",po.get("check_ending",1)); + if(!po.has("tractography_atlas")) + po.set("tractography_atlas",0); + } std::vector tolerance; { if(!po.has("export")) @@ -135,7 +164,11 @@ std::string run_auto_track(tipl::program_option& po,const std::vector { std::shared_ptr fib(new fib_data); set_template(fib,po); - auto list = fib->get_tractography_all_levels(); + if(chen_mode && !fib->load_track_atlas(false)) + return fib->error_msg; + auto list = chen_mode ? fib->tractography_name_list : fib->get_tractography_all_levels(); + if(list.empty()) + return "no tractography atlas available in the selected template"; { std::string labels; for(auto each : list) @@ -147,47 +180,60 @@ std::string run_auto_track(tipl::program_option& po,const std::vector tipl::out() << "available track_ids in current template: " << labels; } auto selections = tipl::split(po.get("track_id","Arcuate,Cingulum,Aslant,InferiorFronto,InferiorLongitudinal,SuperiorLongitudinal,Uncinate,Fornix,Corticos,ThalamicR,Optic,Lemniscus,Reticular,Corpus"),','); - std::vector selected(list.size()); - std::vector backup_subcomponents; - for(const auto& each : selections) + if(chen_mode) { - auto sep_count = std::count(each.begin(),each.end(),'_'); - for(size_t i = 0;i < list.size();++i) + std::vector selected; + for(const auto& each : selections) + selected.push_back(each); + if(selected.empty()) + selected = list; + for(const auto& name : list) + if(is_selected(selected,name)) + tract_name_list.push_back(name); + } + else + { + std::vector selected(list.size()); + std::vector backup_subcomponents; + for(const auto& each : selections) { - if(selected[i]) - continue; - if(tipl::equal_case_insensitive(list[i],each)) - selected[i] = true; - if(tipl::contains_case_insensitive(list[i],each)) + auto sep_count = std::count(each.begin(),each.end(),'_'); + for(size_t i = 0;i < list.size();++i) { - if(std::count(list[i].begin(),list[i].end(),'_') < 2) // not subbundle then contain also work + if(selected[i]) + continue; + if(tipl::equal_case_insensitive(list[i],each)) selected[i] = true; - else - backup_subcomponents.push_back(i); + if(tipl::contains_case_insensitive(list[i],each)) + { + if(std::count(list[i].begin(),list[i].end(),'_') < 2) // not subbundle then contain also work + selected[i] = true; + else + backup_subcomponents.push_back(i); + } } } - } - if(std::all_of(selected.begin(), selected.end(), [](bool s){return !s; })) - { - tipl::out() << "no primary bundle matches. select subcomponents..."; - for(auto each : backup_subcomponents) - selected[each] = true; - } - - std::string selected_list; - for(size_t i = 0;i < list.size();++i) - { - if(selected[i]) + if(std::all_of(selected.begin(), selected.end(), [](bool s){return !s; })) { - if(!selected_list.empty()) - selected_list += ","; - selected_list += list[i]; - tract_name_list.push_back(list[i]); + tipl::out() << "no primary bundle matches. select subcomponents..."; + for(auto each : backup_subcomponents) + selected[each] = true; } + + for(size_t i = 0;i < list.size();++i) + if(selected[i]) + tract_name_list.push_back(list[i]); } if(tract_name_list.empty()) return "cannot find any tract matching --track_id"; + std::string selected_list; + for(size_t i = 0;i < tract_name_list.size();++i) + { + if(i) + selected_list += ","; + selected_list += tract_name_list[i]; + } tipl::out() << "selected tracts: " << selected_list; } @@ -257,7 +303,7 @@ std::string run_auto_track(tipl::program_option& po,const std::vector { ThreadData thread(handle); { - if(!handle->load_track_atlas(true/*symmetric*/)) + if(!handle->load_track_atlas(!chen_mode/*symmetric*/)) return handle->error_msg + " at " + fib_file_name; if (po.has("threshold_index") && !handle->dir.set_tracking_index(po.get("threshold_index"))) @@ -269,19 +315,47 @@ std::string run_auto_track(tipl::program_option& po,const std::vector thread.param.step_size = po.get("step_size",thread.param.step_size); thread.param.smooth_fraction = po.get("smoothing",thread.param.smooth_fraction); - auto minmax = handle->get_track_minmax_length(tract_name); - thread.param.min_length = handle->vs[0]*std::max(tolerance[tracking_iteration], - minmax.first-2.0f*tolerance[tracking_iteration])/handle->tract_atlas_jacobian; - thread.param.max_length = handle->vs[0]*(minmax.second+2.0f*tolerance[tracking_iteration])/handle->tract_atlas_jacobian; + if(chen_mode) + { + auto track_ids = handle->get_track_ids(tract_name); + float min_l = 0.0f,max_l = 0.0f; + for(size_t idx = 0;idx < track_ids.size();++idx) + { + auto id = track_ids[idx]; + if(idx == 0) + { + min_l = handle->tract_atlas_min_length[id]; + max_l = handle->tract_atlas_max_length[id]; + } + else + { + min_l = std::min(min_l,handle->tract_atlas_min_length[id]); + max_l = std::max(max_l,handle->tract_atlas_max_length[id]); + } + } + thread.param.min_length = handle->vs[0]*std::max(tolerance[tracking_iteration], + min_l-2.0f*tolerance[tracking_iteration])/handle->tract_atlas_jacobian; + thread.param.max_length = handle->vs[0]*(max_l+2.0f*tolerance[tracking_iteration])/handle->tract_atlas_jacobian; + } + else + { + auto minmax = handle->get_track_minmax_length(tract_name); + thread.param.min_length = handle->vs[0]*std::max(tolerance[tracking_iteration], + minmax.first-2.0f*tolerance[tracking_iteration])/handle->tract_atlas_jacobian; + thread.param.max_length = handle->vs[0]*(minmax.second+2.0f*tolerance[tracking_iteration])/handle->tract_atlas_jacobian; + } tipl::out() << "min_length(mm): " << thread.param.min_length << std::endl; tipl::out() << "max_length(mm): " << thread.param.max_length << std::endl; - thread.param.tip_iteration = po.get("tip_iteration",32); - thread.param.check_ending = po.get("check_ending",1); - thread.param.track_voxel_ratio = po.get("track_voxel_ratio",thread.param.track_voxel_ratio); + thread.param.tip_iteration = po.get("tip_iteration",chen_mode ? 4 : 32); + if(chen_mode) + thread.param.check_ending = po.get("check_ending",1) && !tipl::contains_case_insensitive(tract_name,"Cingulum"); + else + thread.param.check_ending = po.get("check_ending",1); + thread.param.track_voxel_ratio = po.get("track_voxel_ratio",chen_mode ? 2.0f : thread.param.track_voxel_ratio); } { thread.roi_mgr->use_auto_track = true; - thread.roi_mgr->use_roi = po.get("use_roi",1); + thread.roi_mgr->use_roi = po.get("use_roi",chen_mode ? 0 : 1); thread.roi_mgr->tract_name = tract_name; thread.roi_mgr->tolerance_dis_in_icbm152_mm = tolerance[tracking_iteration]; } From abbbc57ca9fc9a72397800c57a60eeb9ed675374 Mon Sep 17 00:00:00 2001 From: mattcieslak Date: Sun, 18 Jan 2026 12:43:52 -0500 Subject: [PATCH 02/11] don't change tract selection to chen compat --- auto_track.cpp | 85 +++++++++++++++----------------------------------- 1 file changed, 25 insertions(+), 60 deletions(-) diff --git a/auto_track.cpp b/auto_track.cpp index 8d1c0e83..01f0e06d 100644 --- a/auto_track.cpp +++ b/auto_track.cpp @@ -1,8 +1,6 @@ #include #include #include -#include -#include #include "auto_track.h" #include "ui_auto_track.h" #include "libs/dsi/image_model.hpp" @@ -10,22 +8,6 @@ #include "libs/tracking/tracking_thread.hpp" #include extern std::vector fa_template_list; -static bool find_string_case_insensitive(const std::string & str1, const std::string & str2) -{ - auto it = std::search( - str1.begin(), str1.end(), - str2.begin(), str2.end(), - [](char ch1, char ch2) { return std::toupper(ch1) == std::toupper(ch2); } - ); - return (it != str1.end() ); -} -static bool is_selected(std::vector& selected_tracts,const std::string& tract_name) -{ - for(const auto& each: selected_tracts) - if(find_string_case_insensitive(tract_name,each)) - return true; - return false; -} auto_track::auto_track(QWidget *parent) : QMainWindow(parent), ui(new Ui::auto_track) @@ -164,11 +146,7 @@ std::string run_auto_track(tipl::program_option& po,const std::vector { std::shared_ptr fib(new fib_data); set_template(fib,po); - if(chen_mode && !fib->load_track_atlas(false)) - return fib->error_msg; - auto list = chen_mode ? fib->tractography_name_list : fib->get_tractography_all_levels(); - if(list.empty()) - return "no tractography atlas available in the selected template"; + auto list = fib->get_tractography_all_levels(); { std::string labels; for(auto each : list) @@ -180,51 +158,38 @@ std::string run_auto_track(tipl::program_option& po,const std::vector tipl::out() << "available track_ids in current template: " << labels; } auto selections = tipl::split(po.get("track_id","Arcuate,Cingulum,Aslant,InferiorFronto,InferiorLongitudinal,SuperiorLongitudinal,Uncinate,Fornix,Corticos,ThalamicR,Optic,Lemniscus,Reticular,Corpus"),','); - if(chen_mode) - { - std::vector selected; - for(const auto& each : selections) - selected.push_back(each); - if(selected.empty()) - selected = list; - for(const auto& name : list) - if(is_selected(selected,name)) - tract_name_list.push_back(name); - } - else + std::vector selected(list.size()); + std::vector backup_subcomponents; + for(const auto& each : selections) { - std::vector selected(list.size()); - std::vector backup_subcomponents; - for(const auto& each : selections) + auto sep_count = std::count(each.begin(),each.end(),'_'); + for(size_t i = 0;i < list.size();++i) { - auto sep_count = std::count(each.begin(),each.end(),'_'); - for(size_t i = 0;i < list.size();++i) + if(selected[i]) + continue; + if(tipl::equal_case_insensitive(list[i],each)) + selected[i] = true; + if(tipl::contains_case_insensitive(list[i],each)) { - if(selected[i]) - continue; - if(tipl::equal_case_insensitive(list[i],each)) + if(std::count(list[i].begin(),list[i].end(),'_') < 2) // not subbundle then contain also work selected[i] = true; - if(tipl::contains_case_insensitive(list[i],each)) - { - if(std::count(list[i].begin(),list[i].end(),'_') < 2) // not subbundle then contain also work - selected[i] = true; - else - backup_subcomponents.push_back(i); - } + else + backup_subcomponents.push_back(i); } } + } - if(std::all_of(selected.begin(), selected.end(), [](bool s){return !s; })) - { - tipl::out() << "no primary bundle matches. select subcomponents..."; - for(auto each : backup_subcomponents) - selected[each] = true; - } - - for(size_t i = 0;i < list.size();++i) - if(selected[i]) - tract_name_list.push_back(list[i]); + if(std::all_of(selected.begin(), selected.end(), [](bool s){return !s; })) + { + tipl::out() << "no primary bundle matches. select subcomponents..."; + for(auto each : backup_subcomponents) + selected[each] = true; } + + for(size_t i = 0;i < list.size();++i) + if(selected[i]) + tract_name_list.push_back(list[i]); + if(tract_name_list.empty()) return "cannot find any tract matching --track_id"; std::string selected_list; From abdf699244e4bcbcca8ef5c2995a5c653ca6dac3 Mon Sep 17 00:00:00 2001 From: mattcieslak Date: Sun, 18 Jan 2026 16:48:46 -0500 Subject: [PATCH 03/11] attempt at a pro file --- dsi_studio.pro | 188 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 188 insertions(+) create mode 100644 dsi_studio.pro diff --git a/dsi_studio.pro b/dsi_studio.pro new file mode 100644 index 00000000..7486c267 --- /dev/null +++ b/dsi_studio.pro @@ -0,0 +1,188 @@ +QT += core \ + gui \ + opengl \ + charts \ + network + +greaterThan(QT_MAJOR_VERSION, 5): QT += widgets +equals(QT_MAJOR_VERSION, 6): DEFINES += QT6_PATCH + +CONFIG += c++17 +# CONFIG += console + +TARGET = dsi_studio +TEMPLATE = app + +INCLUDEPATH += ./plot +DEFINES += TIPL_USE_QT + +linux* { +QMAKE_CXXFLAGS += -fpermissive -Wno-sign-compare +LIBS += -lGLU -lz +} + +win32* { +INCLUDEPATH += ../include $$[QT_INSTALL_HEADERS]/QtZlib +QMAKE_CXXFLAGS += -wd4244 -wd4267 -wd4018 +LIBS += -lOpenGL32 -lGlu32 +RC_ICONS = dsi_studio.ico +} + +mac { +QMAKE_MACOSX_DEPLOYMENT_TARGET = 10.15 +INCLUDEPATH += /Users/admin/include +ICON = dsi_studio.icns +LIBS += -lz +} + +INCLUDEPATH += libs \ + libs/dsi \ + libs/tracking \ + libs/mapping \ + dicom + +HEADERS += mainwindow.h \ + console.h \ + dicom/dicom_parser.h \ + dicom/dwi_header.hpp \ + libs/dsi/tessellated_icosahedron.hpp \ + libs/dsi/odf_process.hpp \ + libs/dsi/image_model.hpp \ + libs/dsi/gqi_process.hpp \ + libs/dsi/gqi_mni_reconstruction.hpp \ + libs/dsi/dti_process.hpp \ + libs/dsi/basic_voxel.hpp \ + SliceModel.h \ + tracking/tracking_window.h \ + reconstruction/reconstruction_window.h \ + tracking/slice_view_scene.h \ + opengl/glwidget.h \ + opengl/renderingtablewidget.h \ + opengl/tract_render.hpp \ + opengl/region_render.hpp \ + libs/tracking/tracking_method.hpp \ + libs/tracking/roi.hpp \ + libs/tracking/fib_data.hpp \ + libs/tracking/basic_process.hpp \ + libs/tracking/tract_cluster.hpp \ + tracking/region/regiontablewidget.h \ + tracking/region/Regions.h \ + libs/tracking/tract_model.hpp \ + tracking/tract/tracttablewidget.h \ + qcolorcombobox.h \ + libs/tracking/tracking_thread.hpp \ + libs/mapping/atlas.hpp \ + view_image.h \ + manual_alignment.h \ + tracking/tract_report.hpp \ + tracking/color_bar_dialog.hpp \ + tracking/connectivity_matrix_dialog.h \ + tracking/atlasdialog.h \ + filebrowser.h \ + qcompletelineedit.h \ + libs/mapping/connectometry_db.hpp \ + connectometry/createdbdialog.h \ + connectometry/match_db.h \ + connectometry/db_window.h \ + connectometry/group_connectometry.hpp \ + connectometry/group_connectometry_analysis.h \ + regtoolbox.h \ + auto_track.h \ + tracking/device.h \ + tracking/devicetablewidget.h \ + mac_filesystem.hpp \ + libs/dsi/hist_process.hpp \ + xnat_dialog.h + +FORMS += mainwindow.ui \ + console.ui \ + tracking/tracking_window.ui \ + reconstruction/reconstruction_window.ui \ + dicom/dicom_parser.ui \ + view_image.ui \ + manual_alignment.ui \ + tracking/tract_report.ui \ + tracking/color_bar_dialog.ui \ + tracking/connectivity_matrix_dialog.ui \ + tracking/atlasdialog.ui \ + filebrowser.ui \ + connectometry/createdbdialog.ui \ + connectometry/match_db.ui \ + connectometry/db_window.ui \ + connectometry/group_connectometry.ui \ + regtoolbox.ui \ + auto_track.ui \ + xnat_dialog.ui + +RESOURCES += \ + icons.qrc + +SOURCES += main.cpp \ + console.cpp \ + mainwindow.cpp \ + dicom/dicom_parser.cpp \ + dicom/dwi_header.cpp \ + libs/dsi/dsi_interface_imp.cpp \ + libs/dsi/gqi_process.cpp \ + libs/tracking/tract_cluster.cpp \ + SliceModel.cpp \ + tracking/tracking_window.cpp \ + tracking/tracking_window_action.cpp \ + reconstruction/reconstruction_window.cpp \ + tracking/slice_view_scene.cpp \ + opengl/glwidget.cpp \ + opengl/tract_render.cpp \ + opengl/region_render.cpp \ + opengl/renderingtablewidget.cpp \ + tracking/region/regiontablewidget.cpp \ + tracking/region/Regions.cpp \ + libs/tracking/tract_model.cpp \ + tracking/tract/tracttablewidget.cpp \ + qcolorcombobox.cpp \ + cmd/trk.cpp \ + cmd/rec.cpp \ + cmd/src.cpp \ + libs/mapping/atlas.cpp \ + cmd/ana.cpp \ + view_image.cpp \ + manual_alignment.cpp \ + tracking/tract_report.cpp \ + tracking/color_bar_dialog.cpp \ + cmd/exp.cpp \ + tracking/connectivity_matrix_dialog.cpp \ + libs/dsi/tessellated_icosahedron.cpp \ + cmd/atl.cpp \ + tracking/atlasdialog.cpp \ + cmd/cnt.cpp \ + cmd/vis.cpp \ + filebrowser.cpp \ + qcompletelineedit.cpp \ + libs/tracking/fib_data.cpp \ + libs/tracking/tracking_thread.cpp \ + cmd/ren.cpp \ + libs/mapping/connectometry_db.cpp \ + connectometry/createdbdialog.cpp \ + connectometry/match_db.cpp \ + connectometry/db_window.cpp \ + connectometry/group_connectometry.cpp \ + connectometry/group_connectometry_analysis.cpp \ + regtoolbox.cpp \ + cmd/cnn.cpp \ + cmd/qc.cpp \ + libs/dsi/basic_voxel.cpp \ + libs/dsi/image_model.cpp \ + cmd/reg.cpp \ + auto_track.cpp \ + cmd/atk.cpp \ + tracking/device.cpp \ + tracking/devicetablewidget.cpp \ + cmd/xnat.cpp \ + xnat_dialog.cpp + +OTHER_FILES += \ + options.txt \ + dicom_tag.txt \ + FreeSurferColorLUT.txt + +DISTFILES += \ + LICENSE From 81a5685d1c6c6a72c4a981eb1545fa951596498b Mon Sep 17 00:00:00 2001 From: mattcieslak Date: Sun, 18 Jan 2026 16:57:43 -0500 Subject: [PATCH 04/11] remove gui components --- dsi_studio.pro | 7 ------- 1 file changed, 7 deletions(-) diff --git a/dsi_studio.pro b/dsi_studio.pro index 7486c267..da24d81d 100644 --- a/dsi_studio.pro +++ b/dsi_studio.pro @@ -75,10 +75,8 @@ HEADERS += mainwindow.h \ view_image.h \ manual_alignment.h \ tracking/tract_report.hpp \ - tracking/color_bar_dialog.hpp \ tracking/connectivity_matrix_dialog.h \ tracking/atlasdialog.h \ - filebrowser.h \ qcompletelineedit.h \ libs/mapping/connectometry_db.hpp \ connectometry/createdbdialog.h \ @@ -90,7 +88,6 @@ HEADERS += mainwindow.h \ auto_track.h \ tracking/device.h \ tracking/devicetablewidget.h \ - mac_filesystem.hpp \ libs/dsi/hist_process.hpp \ xnat_dialog.h @@ -102,10 +99,8 @@ FORMS += mainwindow.ui \ view_image.ui \ manual_alignment.ui \ tracking/tract_report.ui \ - tracking/color_bar_dialog.ui \ tracking/connectivity_matrix_dialog.ui \ tracking/atlasdialog.ui \ - filebrowser.ui \ connectometry/createdbdialog.ui \ connectometry/match_db.ui \ connectometry/db_window.ui \ @@ -147,7 +142,6 @@ SOURCES += main.cpp \ view_image.cpp \ manual_alignment.cpp \ tracking/tract_report.cpp \ - tracking/color_bar_dialog.cpp \ cmd/exp.cpp \ tracking/connectivity_matrix_dialog.cpp \ libs/dsi/tessellated_icosahedron.cpp \ @@ -155,7 +149,6 @@ SOURCES += main.cpp \ tracking/atlasdialog.cpp \ cmd/cnt.cpp \ cmd/vis.cpp \ - filebrowser.cpp \ qcompletelineedit.cpp \ libs/tracking/fib_data.cpp \ libs/tracking/tracking_thread.cpp \ From 4f8c26af0a516da28d1c57f50700428c34ca5666 Mon Sep 17 00:00:00 2001 From: mattcieslak Date: Sun, 18 Jan 2026 18:26:27 -0500 Subject: [PATCH 05/11] Update pro --- dsi_studio.pro | 2 ++ 1 file changed, 2 insertions(+) diff --git a/dsi_studio.pro b/dsi_studio.pro index da24d81d..11907839 100644 --- a/dsi_studio.pro +++ b/dsi_studio.pro @@ -117,6 +117,7 @@ SOURCES += main.cpp \ mainwindow.cpp \ dicom/dicom_parser.cpp \ dicom/dwi_header.cpp \ + cmd/img.cpp \ libs/dsi/dsi_interface_imp.cpp \ libs/dsi/gqi_process.cpp \ libs/tracking/tract_cluster.cpp \ @@ -150,6 +151,7 @@ SOURCES += main.cpp \ cmd/cnt.cpp \ cmd/vis.cpp \ qcompletelineedit.cpp \ + libs/tracking/roi.cpp \ libs/tracking/fib_data.cpp \ libs/tracking/tracking_thread.cpp \ cmd/ren.cpp \ From ade501374477b1cedbb47dc593742aad2d077330 Mon Sep 17 00:00:00 2001 From: mattcieslak Date: Tue, 20 Jan 2026 13:31:18 -0500 Subject: [PATCH 06/11] update registratio for chen mode --- auto_track.cpp | 2 + libs/tracking/fib_data.cpp | 209 +++++++++++++++++++++++++++++++++++++ libs/tracking/fib_data.hpp | 2 + reg.hpp | 59 +++++++++++ 4 files changed, 272 insertions(+) diff --git a/auto_track.cpp b/auto_track.cpp index 01f0e06d..8c90a8e6 100644 --- a/auto_track.cpp +++ b/auto_track.cpp @@ -145,6 +145,7 @@ std::string run_auto_track(tipl::program_option& po,const std::vector std::vector tract_name_list; { std::shared_ptr fib(new fib_data); + fib->use_chen_normalization = chen_mode; set_template(fib,po); auto list = fib->get_tractography_all_levels(); { @@ -256,6 +257,7 @@ std::string run_auto_track(tipl::program_option& po,const std::vector handle = std::make_shared(); if(!handle->load_from_file(fib_file_name.c_str())) return handle->error_msg; + handle->use_chen_normalization = chen_mode; set_template(handle,po); } std::shared_ptr tract_model(new TractModel(handle)); diff --git a/libs/tracking/fib_data.cpp b/libs/tracking/fib_data.cpp index 3fcb5b3d..f494203b 100644 --- a/libs/tracking/fib_data.cpp +++ b/libs/tracking/fib_data.cpp @@ -2033,6 +2033,14 @@ bool to_t1wt2w_templates(dual_reg& reg,size_t template_id,bool be) } std::string fib_data::get_mapping_file_name(void) const { + if(use_chen_normalization) + { + std::string output_file_name(fib_file_name); + output_file_name += "."; + output_file_name += QFileInfo(fa_template_list[template_id].c_str()).baseName().toLower().toStdString(); + output_file_name += ".map.gz"; + return output_file_name; + } std::string output_file_name(fib_file_name + "." + std::filesystem::path(fa_template_list[template_id]).filename().stem().stem().stem().string()); if(alternative_mapping_index) @@ -2040,8 +2048,209 @@ std::string fib_data::get_mapping_file_name(void) const output_file_name += ".mz"; return output_file_name; } +namespace +{ +void chen_cdm_common(const tipl::image<3>& It,const tipl::image<3>& It2, + const tipl::image<3>& Is,const tipl::image<3>& Is2, + tipl::image<3,tipl::vector<3> >& dis, + tipl::image<3,tipl::vector<3> >& inv_dis, + bool& terminated) +{ + tipl::reg::cdm_param param; + tipl::reg::cdm_common(tipl::reg::make_list(It,It2), + tipl::reg::make_list(Is,Is2), + dis,terminated,param,has_cuda); + if(terminated) + return; + tipl::reg::cdm_common(tipl::reg::make_list(Is,Is2), + tipl::reg::make_list(It,It2), + inv_dis,terminated,param,has_cuda); +} +} +bool fib_data::map_to_mni_chen(bool background) +{ + if(!load_template()) + return false; + if(is_mni && template_id == matched_template_id) + return true; + if(!s2t.empty() && !t2s.empty()) + return true; + auto output_file_name = get_mapping_file_name(); + + tipl::io::gz_mat_read in; + + // check 1. mapping files was created later than the FIB file + // 2. the recon steps are the same + // 3. has inv_mapping matrix (new format after Sep 2021 + if(QFileInfo(output_file_name.c_str()).lastModified() > QFileInfo(fib_file_name.c_str()).lastModified() && + in.load_from_file(output_file_name.c_str()) && in.has("from2to") && in.has("to2from") + && in.read("steps") == steps) + { + const float* t2s_ptr = nullptr; + unsigned int t2s_row,t2s_col,s2t_row,s2t_col; + const float* s2t_ptr = nullptr; + if(in.read("to2from",t2s_row,t2s_col,t2s_ptr) && + in.read("from2to",s2t_row,s2t_col,s2t_ptr)) + { + if(size_t(t2s_col)*size_t(t2s_row) == template_I.size()*3 && + size_t(s2t_col)*size_t(s2t_row) == dim.size()*3) + { + tipl::out() << "loading mapping fields from " << output_file_name << std::endl; + t2s.clear(); + t2s.resize(template_I.shape()); + s2t.clear(); + s2t.resize(dim); + std::copy(t2s_ptr,t2s_ptr+t2s_col*t2s_row,&t2s[0][0]); + std::copy(s2t_ptr,s2t_ptr+s2t_col*s2t_row,&s2t[0][0]); + prog = 6; + return true; + } + } + } + tipl::progress prog_("running normalization"); + prog = 0; + bool terminated = false; + auto lambda = [this,output_file_name,&terminated]() + { + tipl::transformation_matrix T; + + auto It = template_I; + auto It2 = template_I2; + + if(dir.index_name[0] == "image" && // not FIB file + tipl::io::gz_nifti::load_to_space(t1w_template_file_name.c_str(),It,template_to_mni)) // not FIB file, use t1w as template + { + tipl::out() << "using structure image for normalization" << std::endl; + It2.clear(); + } + + tipl::image<3> Is(dir.fa[0],dim); + tipl::image<3> Is2; + + { + size_t iso_index = get_name_index("iso"); + if(slices.size() == iso_index) + iso_index = get_name_index("md"); + if(slices.size() != iso_index) + Is2 = slices[iso_index]->get_image(); + } + bool no_iso = Is2.empty() || It2.empty(); + + tipl::filter::gaussian(Is); + tipl::filter::gaussian(Is); + if(!no_iso) + { + tipl::filter::gaussian(Is2); + tipl::filter::gaussian(Is2); + } + + prog = 1; + if(!has_manual_atlas) + linear_with_mi(It,template_vs,Is,vs,T,tipl::reg::affine,terminated); + else + T = manual_template_T; + + if(terminated) + return; + prog = 2; + tipl::image<3> Iss(It.shape()); + tipl::resample_mt(Is,Iss,T); + tipl::image<3> Iss2(It2.shape()); + if(!no_iso) + tipl::resample_mt(Is2,Iss2,T); + prog = 3; + + if(dir.index_name[0] == "image") + { + tipl::out() << "matching t1w t2w contrast" << std::endl; + tipl::image<3> t2w(It.shape()); + if(tipl::io::gz_nifti::load_to_space(t2w_template_file_name.c_str(),t2w,template_to_mni)) + { + std::vector X(It.size()*3); + tipl::par_for(It.size(),[&](size_t pos) + { + if(Iss[pos] == 0.0f) + return; + size_t i = pos; + pos = pos+pos+pos; + X[pos] = 1; + X[pos+1] = It[i]; + X[pos+2] = t2w[i]; + }); + tipl::multiple_regression m; + if(m.set_variables(X.begin(),3,It.size())) + { + float b[3] = {0.0f,0.0f,0.0f}; + m.regress(Iss.begin(),b); + tipl::out() << "image=" << b[0] << " + " << b[1] << " × t1w + " << b[2] << " × t2w "; + tipl::par_for(It.size(),[&](size_t pos) + { + if(Iss[pos] == 0.0f) + return; + It[pos] = b[0] + b[1]*It[pos] + b[2]*t2w[pos]; + if(It[pos] < 0.0f) + It[pos] = 0.0f; + }); + } + } + } + tipl::image<3,tipl::vector<3> > dis,inv_dis; + tipl::reg::cdm_pre(It,It2,Iss,Iss2); + + chen_cdm_common(It,It2,Iss,Iss2,dis,inv_dis,terminated); + + tipl::displacement_to_mapping(dis,t2s,T); + tipl::compose_mapping(Is,t2s,Iss); + tipl::out() << "R2:" << tipl::correlation(Iss.begin(),Iss.end(),It.begin()) << std::endl; + + s2t.resize(dim); + tipl::inv_displacement_to_mapping(inv_dis,s2t,T); + + if(terminated) + return; + tipl::io::gz_mat_write out(output_file_name.c_str()); + if(out) + { + out.write("to_dim",dis.shape()); + out.write("to_vs",template_vs); + out.write("from_dim",dim); + out.write("from_vs",vs); + out.write("steps",steps); + prog = 4; + tipl::out() << "calculating template to subject warp field"; + out.write("to2from",&t2s[0][0],3,t2s.size()); + prog = 5; + tipl::out() << "calculating subject to template warp field"; + out.write("from2to",&s2t[0][0],3,s2t.size()); + } + + prog = 6; + }; + + if(background) + { + std::thread t(lambda); + while(prog_(prog,6)) + std::this_thread::yield(); + if(prog_.aborted()) + { + error_msg = "aborted."; + terminated = true; + } + t.join(); + return !prog_.aborted(); + } + else + { + tipl::out() << "Subject normalization to MNI space." << std::endl; + lambda(); + } + return true; +} bool fib_data::map_to_mni(bool background) { + if(use_chen_normalization) + return map_to_mni_chen(background); if(!load_template()) return false; if(is_mni && template_id == matched_template_id) diff --git a/libs/tracking/fib_data.hpp b/libs/tracking/fib_data.hpp index 0d27a339..be2d8fe8 100644 --- a/libs/tracking/fib_data.hpp +++ b/libs/tracking/fib_data.hpp @@ -226,6 +226,7 @@ class fib_data bool is_histology = true; bool is_mni = false; bool is_be = false; + bool use_chen_normalization = false; size_t matched_template_id = 0; bool trackable = true; float min_length(void) const @@ -308,6 +309,7 @@ class fib_data public: std::string get_mapping_file_name(void) const; bool map_to_mni(bool background = true); + bool map_to_mni_chen(bool background = true); void temp2sub(std::vector >&tracts) const; void temp2sub(tipl::vector<3>& pos) const; void sub2temp(tipl::vector<3>& pos); diff --git a/reg.hpp b/reg.hpp index c424d2f5..59c3058b 100644 --- a/reg.hpp +++ b/reg.hpp @@ -30,6 +30,46 @@ inline auto template_image_pre(const tipl::image& I) return template_image_pre(tipl::image(I)); } +inline size_t linear_with_mi(const tipl::image<3,float>& from, + const tipl::vector<3>& from_vs, + const tipl::image<3,float>& to, + const tipl::vector<3>& to_vs, + tipl::affine_transform& arg, + tipl::reg::reg_type reg_type, + bool& terminated, + const float* bound = tipl::reg::reg_bound) +{ + auto new_to_vs = to_vs; + if(reg_type == tipl::reg::affine) + new_to_vs = tipl::reg::adjust_to_vs(from,from_vs,to,to_vs); + if(new_to_vs != to_vs) + tipl::transformation_matrix(arg,from,from_vs,to,to_vs).to_affine_transform(arg,from,from_vs,to,new_to_vs); + float result = tipl::reg::linear_reg(tipl::reg::make_list(from),from_vs, + tipl::reg::make_list(to),new_to_vs, + arg,bound,reg_type,tipl::reg::mutual_info, + true,has_cuda,terminated); + if(new_to_vs != to_vs) + tipl::transformation_matrix(arg,from,from_vs,to,new_to_vs).to_affine_transform(arg,from,from_vs,to,to_vs); + return result; +} + +inline size_t linear_with_mi(const tipl::image<3,float>& from, + const tipl::vector<3>& from_vs, + const tipl::image<3,float>& to, + const tipl::vector<3>& to_vs, + tipl::transformation_matrix& T, + tipl::reg::reg_type reg_type, + bool& terminated, + const float* bound = tipl::reg::reg_bound) +{ + tipl::affine_transform arg; + size_t result = linear_with_mi(from,from_vs,to,to_vs,arg,reg_type,terminated,bound); + tipl::out() << arg << std::endl; + T = tipl::transformation_matrix(arg,from.shape(),from_vs,to.shape(),to_vs); + tipl::out() << T << std::endl; + return result; +} + extern int map_ver; struct dual_reg{ static constexpr int dimension = 3; @@ -137,4 +177,23 @@ struct dual_reg{ +namespace tipl +{ +namespace reg +{ +template +inline void cdm_pre(image_type& It,image_type& It2,image_type& Is,image_type& Is2) +{ + if(!It.empty()) + tipl::normalize_upper_lower2(It,It,255.99f); + if(!Is.empty()) + tipl::normalize_upper_lower2(Is,Is,255.99f); + if(!It2.empty()) + tipl::normalize_upper_lower2(It2,It2,255.99f); + if(!Is2.empty()) + tipl::normalize_upper_lower2(Is2,Is2,255.99f); +} +} +} + #endif//REG_HPP From 13ef858e3481ba86ebcc0054264101efbf507f34 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Wed, 21 Jan 2026 16:54:38 -0500 Subject: [PATCH 07/11] seems like it works? --- auto_track.cpp | 2 -- libs/tracking/fib_data.cpp | 10 +++++----- reg.hpp | 8 ++++---- 3 files changed, 9 insertions(+), 11 deletions(-) diff --git a/auto_track.cpp b/auto_track.cpp index 8c90a8e6..ef4c0b45 100644 --- a/auto_track.cpp +++ b/auto_track.cpp @@ -119,8 +119,6 @@ std::string run_auto_track(tipl::program_option& po,const std::vector po.set("tip_iteration",po.get("tip_iteration",4)); po.set("use_roi",po.get("use_roi",0)); po.set("check_ending",po.get("check_ending",1)); - if(!po.has("tractography_atlas")) - po.set("tractography_atlas",0); } std::vector tolerance; { diff --git a/libs/tracking/fib_data.cpp b/libs/tracking/fib_data.cpp index f494203b..583183c5 100644 --- a/libs/tracking/fib_data.cpp +++ b/libs/tracking/fib_data.cpp @@ -2118,7 +2118,7 @@ bool fib_data::map_to_mni_chen(bool background) auto It2 = template_I2; if(dir.index_name[0] == "image" && // not FIB file - tipl::io::gz_nifti::load_to_space(t1w_template_file_name.c_str(),It,template_to_mni)) // not FIB file, use t1w as template + tipl::io::gz_nifti(t1w_template_file_name.c_str(),std::ios::in).to_space(It,template_to_mni)) // not FIB file, use t1w as template { tipl::out() << "using structure image for normalization" << std::endl; It2.clear(); @@ -2148,23 +2148,23 @@ bool fib_data::map_to_mni_chen(bool background) if(!has_manual_atlas) linear_with_mi(It,template_vs,Is,vs,T,tipl::reg::affine,terminated); else - T = manual_template_T; + T = tipl::transformation_matrix(manual_template_T,It.shape(),template_vs,Is.shape(),vs); if(terminated) return; prog = 2; tipl::image<3> Iss(It.shape()); - tipl::resample_mt(Is,Iss,T); + tipl::resample(Is,Iss,T); tipl::image<3> Iss2(It2.shape()); if(!no_iso) - tipl::resample_mt(Is2,Iss2,T); + tipl::resample(Is2,Iss2,T); prog = 3; if(dir.index_name[0] == "image") { tipl::out() << "matching t1w t2w contrast" << std::endl; tipl::image<3> t2w(It.shape()); - if(tipl::io::gz_nifti::load_to_space(t2w_template_file_name.c_str(),t2w,template_to_mni)) + if(tipl::io::gz_nifti(t2w_template_file_name.c_str(),std::ios::in).to_space(t2w,template_to_mni)) { std::vector X(It.size()*3); tipl::par_for(It.size(),[&](size_t pos) diff --git a/reg.hpp b/reg.hpp index 59c3058b..7b9061fb 100644 --- a/reg.hpp +++ b/reg.hpp @@ -37,19 +37,19 @@ inline size_t linear_with_mi(const tipl::image<3,float>& from, tipl::affine_transform& arg, tipl::reg::reg_type reg_type, bool& terminated, - const float* bound = tipl::reg::reg_bound) + const float (*bound)[8] = tipl::reg::reg_bound) { auto new_to_vs = to_vs; if(reg_type == tipl::reg::affine) new_to_vs = tipl::reg::adjust_to_vs(from,from_vs,to,to_vs); if(new_to_vs != to_vs) - tipl::transformation_matrix(arg,from,from_vs,to,to_vs).to_affine_transform(arg,from,from_vs,to,new_to_vs); + arg = tipl::transformation_matrix(arg,from,from_vs,to,to_vs).to_affine_transform(from,from_vs,to,new_to_vs); float result = tipl::reg::linear_reg(tipl::reg::make_list(from),from_vs, tipl::reg::make_list(to),new_to_vs, arg,bound,reg_type,tipl::reg::mutual_info, true,has_cuda,terminated); if(new_to_vs != to_vs) - tipl::transformation_matrix(arg,from,from_vs,to,new_to_vs).to_affine_transform(arg,from,from_vs,to,to_vs); + arg = tipl::transformation_matrix(arg,from,from_vs,to,new_to_vs).to_affine_transform(from,from_vs,to,to_vs); return result; } @@ -60,7 +60,7 @@ inline size_t linear_with_mi(const tipl::image<3,float>& from, tipl::transformation_matrix& T, tipl::reg::reg_type reg_type, bool& terminated, - const float* bound = tipl::reg::reg_bound) + const float (*bound)[8] = tipl::reg::reg_bound) { tipl::affine_transform arg; size_t result = linear_with_mi(from,from_vs,to,to_vs,arg,reg_type,terminated,bound); From 939362e462b2d12752fe6b73ef3df507e653dbc8 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Wed, 21 Jan 2026 17:01:01 -0500 Subject: [PATCH 08/11] remove .pro --- dsi_studio.pro | 183 ------------------------------------------------- 1 file changed, 183 deletions(-) delete mode 100644 dsi_studio.pro diff --git a/dsi_studio.pro b/dsi_studio.pro deleted file mode 100644 index 11907839..00000000 --- a/dsi_studio.pro +++ /dev/null @@ -1,183 +0,0 @@ -QT += core \ - gui \ - opengl \ - charts \ - network - -greaterThan(QT_MAJOR_VERSION, 5): QT += widgets -equals(QT_MAJOR_VERSION, 6): DEFINES += QT6_PATCH - -CONFIG += c++17 -# CONFIG += console - -TARGET = dsi_studio -TEMPLATE = app - -INCLUDEPATH += ./plot -DEFINES += TIPL_USE_QT - -linux* { -QMAKE_CXXFLAGS += -fpermissive -Wno-sign-compare -LIBS += -lGLU -lz -} - -win32* { -INCLUDEPATH += ../include $$[QT_INSTALL_HEADERS]/QtZlib -QMAKE_CXXFLAGS += -wd4244 -wd4267 -wd4018 -LIBS += -lOpenGL32 -lGlu32 -RC_ICONS = dsi_studio.ico -} - -mac { -QMAKE_MACOSX_DEPLOYMENT_TARGET = 10.15 -INCLUDEPATH += /Users/admin/include -ICON = dsi_studio.icns -LIBS += -lz -} - -INCLUDEPATH += libs \ - libs/dsi \ - libs/tracking \ - libs/mapping \ - dicom - -HEADERS += mainwindow.h \ - console.h \ - dicom/dicom_parser.h \ - dicom/dwi_header.hpp \ - libs/dsi/tessellated_icosahedron.hpp \ - libs/dsi/odf_process.hpp \ - libs/dsi/image_model.hpp \ - libs/dsi/gqi_process.hpp \ - libs/dsi/gqi_mni_reconstruction.hpp \ - libs/dsi/dti_process.hpp \ - libs/dsi/basic_voxel.hpp \ - SliceModel.h \ - tracking/tracking_window.h \ - reconstruction/reconstruction_window.h \ - tracking/slice_view_scene.h \ - opengl/glwidget.h \ - opengl/renderingtablewidget.h \ - opengl/tract_render.hpp \ - opengl/region_render.hpp \ - libs/tracking/tracking_method.hpp \ - libs/tracking/roi.hpp \ - libs/tracking/fib_data.hpp \ - libs/tracking/basic_process.hpp \ - libs/tracking/tract_cluster.hpp \ - tracking/region/regiontablewidget.h \ - tracking/region/Regions.h \ - libs/tracking/tract_model.hpp \ - tracking/tract/tracttablewidget.h \ - qcolorcombobox.h \ - libs/tracking/tracking_thread.hpp \ - libs/mapping/atlas.hpp \ - view_image.h \ - manual_alignment.h \ - tracking/tract_report.hpp \ - tracking/connectivity_matrix_dialog.h \ - tracking/atlasdialog.h \ - qcompletelineedit.h \ - libs/mapping/connectometry_db.hpp \ - connectometry/createdbdialog.h \ - connectometry/match_db.h \ - connectometry/db_window.h \ - connectometry/group_connectometry.hpp \ - connectometry/group_connectometry_analysis.h \ - regtoolbox.h \ - auto_track.h \ - tracking/device.h \ - tracking/devicetablewidget.h \ - libs/dsi/hist_process.hpp \ - xnat_dialog.h - -FORMS += mainwindow.ui \ - console.ui \ - tracking/tracking_window.ui \ - reconstruction/reconstruction_window.ui \ - dicom/dicom_parser.ui \ - view_image.ui \ - manual_alignment.ui \ - tracking/tract_report.ui \ - tracking/connectivity_matrix_dialog.ui \ - tracking/atlasdialog.ui \ - connectometry/createdbdialog.ui \ - connectometry/match_db.ui \ - connectometry/db_window.ui \ - connectometry/group_connectometry.ui \ - regtoolbox.ui \ - auto_track.ui \ - xnat_dialog.ui - -RESOURCES += \ - icons.qrc - -SOURCES += main.cpp \ - console.cpp \ - mainwindow.cpp \ - dicom/dicom_parser.cpp \ - dicom/dwi_header.cpp \ - cmd/img.cpp \ - libs/dsi/dsi_interface_imp.cpp \ - libs/dsi/gqi_process.cpp \ - libs/tracking/tract_cluster.cpp \ - SliceModel.cpp \ - tracking/tracking_window.cpp \ - tracking/tracking_window_action.cpp \ - reconstruction/reconstruction_window.cpp \ - tracking/slice_view_scene.cpp \ - opengl/glwidget.cpp \ - opengl/tract_render.cpp \ - opengl/region_render.cpp \ - opengl/renderingtablewidget.cpp \ - tracking/region/regiontablewidget.cpp \ - tracking/region/Regions.cpp \ - libs/tracking/tract_model.cpp \ - tracking/tract/tracttablewidget.cpp \ - qcolorcombobox.cpp \ - cmd/trk.cpp \ - cmd/rec.cpp \ - cmd/src.cpp \ - libs/mapping/atlas.cpp \ - cmd/ana.cpp \ - view_image.cpp \ - manual_alignment.cpp \ - tracking/tract_report.cpp \ - cmd/exp.cpp \ - tracking/connectivity_matrix_dialog.cpp \ - libs/dsi/tessellated_icosahedron.cpp \ - cmd/atl.cpp \ - tracking/atlasdialog.cpp \ - cmd/cnt.cpp \ - cmd/vis.cpp \ - qcompletelineedit.cpp \ - libs/tracking/roi.cpp \ - libs/tracking/fib_data.cpp \ - libs/tracking/tracking_thread.cpp \ - cmd/ren.cpp \ - libs/mapping/connectometry_db.cpp \ - connectometry/createdbdialog.cpp \ - connectometry/match_db.cpp \ - connectometry/db_window.cpp \ - connectometry/group_connectometry.cpp \ - connectometry/group_connectometry_analysis.cpp \ - regtoolbox.cpp \ - cmd/cnn.cpp \ - cmd/qc.cpp \ - libs/dsi/basic_voxel.cpp \ - libs/dsi/image_model.cpp \ - cmd/reg.cpp \ - auto_track.cpp \ - cmd/atk.cpp \ - tracking/device.cpp \ - tracking/devicetablewidget.cpp \ - cmd/xnat.cpp \ - xnat_dialog.cpp - -OTHER_FILES += \ - options.txt \ - dicom_tag.txt \ - FreeSurferColorLUT.txt - -DISTFILES += \ - LICENSE From 4d907cdcdf37dac536bcf237926409a8d5589d8a Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Thu, 22 Jan 2026 12:50:40 -0500 Subject: [PATCH 09/11] remove chen registration code --- auto_track.cpp | 42 +++++++- libs/tracking/fib_data.cpp | 209 ------------------------------------- libs/tracking/fib_data.hpp | 2 - reg.hpp | 12 --- 4 files changed, 37 insertions(+), 228 deletions(-) diff --git a/auto_track.cpp b/auto_track.cpp index ef4c0b45..bbf34bb2 100644 --- a/auto_track.cpp +++ b/auto_track.cpp @@ -114,7 +114,6 @@ std::string run_auto_track(tipl::program_option& po,const std::vector uint32_t thread_count = tipl::max_thread_count; if(chen_mode) { - po.set("template",po.get("template",0)); po.set("track_voxel_ratio",po.get("track_voxel_ratio",2.0f)); po.set("tip_iteration",po.get("tip_iteration",4)); po.set("use_roi",po.get("use_roi",0)); @@ -143,7 +142,6 @@ std::string run_auto_track(tipl::program_option& po,const std::vector std::vector tract_name_list; { std::shared_ptr fib(new fib_data); - fib->use_chen_normalization = chen_mode; set_template(fib,po); auto list = fib->get_tractography_all_levels(); { @@ -214,6 +212,42 @@ std::string run_auto_track(tipl::program_option& po,const std::vector scan_names.push_back(cur_file_base_name); tipl::out() << "processing " << cur_file_base_name << std::endl; std::shared_ptr handle; + std::string fib_base = QFileInfo(fib_file_name.c_str()).baseName().toStdString(); + bool save_atlas_track = po.has("template_tracks_in_native_space"); + std::string atlas_track_option; + if(save_atlas_track) + atlas_track_option = po.get("template_tracks_in_native_space",std::string()); + + if(save_atlas_track) + { + if(!handle.get()) + { + handle = std::make_shared(); + if(!handle->load_from_file(fib_file_name.c_str())) + return handle->error_msg; + set_template(handle,po); + } + if(!handle->load_track_atlas(true/*symmetric*/)) + return handle->error_msg + " at " + fib_file_name; + std::string atlas_track_name = atlas_track_option; + if(atlas_track_name.empty()) + atlas_track_name = dir + "/" + fib_base + ".atlas.tt.gz"; + else if(std::filesystem::is_directory(atlas_track_name)) + atlas_track_name = (std::filesystem::path(atlas_track_name) / (fib_base + ".atlas.tt.gz")).string(); + tipl::out() << "saving template tracks in native space to " << atlas_track_name; + if(!handle->track_atlas || !handle->track_atlas->save_tracts_to_file(atlas_track_name)) + return std::string("failed to save atlas track to ") + atlas_track_name; + auto atlas_label_file = handle->tractography_atlas_file_name + ".txt"; + auto output_label_file = atlas_track_name + ".txt"; + if(std::filesystem::exists(atlas_label_file)) + { + std::error_code ec; + std::filesystem::copy_file(atlas_label_file,output_label_file, + std::filesystem::copy_options::overwrite_existing,ec); + if(ec) + return std::string("failed to copy atlas label file to ") + output_label_file; + } + } tipl::progress prog1("tracking pathways"); for(size_t j = 0;prog1(j,tract_name_list.size());++j) @@ -228,7 +262,6 @@ std::string run_auto_track(tipl::program_option& po,const std::vector if (!dir.exists() && !dir.mkpath(".")) tipl::out() << std::string("cannot create directory: ") + output_path << std::endl; } - std::string fib_base = QFileInfo(fib_file_name.c_str()).baseName().toStdString(); std::string trk_base = output_path + "/" + fib_base+"."+tract_name; std::string no_result_file_name = trk_base+".no_result.txt"; std::string trk_file_name = trk_base + "." + trk_format; @@ -255,7 +288,6 @@ std::string run_auto_track(tipl::program_option& po,const std::vector handle = std::make_shared(); if(!handle->load_from_file(fib_file_name.c_str())) return handle->error_msg; - handle->use_chen_normalization = chen_mode; set_template(handle,po); } std::shared_ptr tract_model(new TractModel(handle)); @@ -268,7 +300,7 @@ std::string run_auto_track(tipl::program_option& po,const std::vector { ThreadData thread(handle); { - if(!handle->load_track_atlas(!chen_mode/*symmetric*/)) + if(!handle->load_track_atlas(true/*symmetric*/)) return handle->error_msg + " at " + fib_file_name; if (po.has("threshold_index") && !handle->dir.set_tracking_index(po.get("threshold_index"))) diff --git a/libs/tracking/fib_data.cpp b/libs/tracking/fib_data.cpp index 583183c5..3fcb5b3d 100644 --- a/libs/tracking/fib_data.cpp +++ b/libs/tracking/fib_data.cpp @@ -2033,14 +2033,6 @@ bool to_t1wt2w_templates(dual_reg& reg,size_t template_id,bool be) } std::string fib_data::get_mapping_file_name(void) const { - if(use_chen_normalization) - { - std::string output_file_name(fib_file_name); - output_file_name += "."; - output_file_name += QFileInfo(fa_template_list[template_id].c_str()).baseName().toLower().toStdString(); - output_file_name += ".map.gz"; - return output_file_name; - } std::string output_file_name(fib_file_name + "." + std::filesystem::path(fa_template_list[template_id]).filename().stem().stem().stem().string()); if(alternative_mapping_index) @@ -2048,209 +2040,8 @@ std::string fib_data::get_mapping_file_name(void) const output_file_name += ".mz"; return output_file_name; } -namespace -{ -void chen_cdm_common(const tipl::image<3>& It,const tipl::image<3>& It2, - const tipl::image<3>& Is,const tipl::image<3>& Is2, - tipl::image<3,tipl::vector<3> >& dis, - tipl::image<3,tipl::vector<3> >& inv_dis, - bool& terminated) -{ - tipl::reg::cdm_param param; - tipl::reg::cdm_common(tipl::reg::make_list(It,It2), - tipl::reg::make_list(Is,Is2), - dis,terminated,param,has_cuda); - if(terminated) - return; - tipl::reg::cdm_common(tipl::reg::make_list(Is,Is2), - tipl::reg::make_list(It,It2), - inv_dis,terminated,param,has_cuda); -} -} -bool fib_data::map_to_mni_chen(bool background) -{ - if(!load_template()) - return false; - if(is_mni && template_id == matched_template_id) - return true; - if(!s2t.empty() && !t2s.empty()) - return true; - auto output_file_name = get_mapping_file_name(); - - tipl::io::gz_mat_read in; - - // check 1. mapping files was created later than the FIB file - // 2. the recon steps are the same - // 3. has inv_mapping matrix (new format after Sep 2021 - if(QFileInfo(output_file_name.c_str()).lastModified() > QFileInfo(fib_file_name.c_str()).lastModified() && - in.load_from_file(output_file_name.c_str()) && in.has("from2to") && in.has("to2from") - && in.read("steps") == steps) - { - const float* t2s_ptr = nullptr; - unsigned int t2s_row,t2s_col,s2t_row,s2t_col; - const float* s2t_ptr = nullptr; - if(in.read("to2from",t2s_row,t2s_col,t2s_ptr) && - in.read("from2to",s2t_row,s2t_col,s2t_ptr)) - { - if(size_t(t2s_col)*size_t(t2s_row) == template_I.size()*3 && - size_t(s2t_col)*size_t(s2t_row) == dim.size()*3) - { - tipl::out() << "loading mapping fields from " << output_file_name << std::endl; - t2s.clear(); - t2s.resize(template_I.shape()); - s2t.clear(); - s2t.resize(dim); - std::copy(t2s_ptr,t2s_ptr+t2s_col*t2s_row,&t2s[0][0]); - std::copy(s2t_ptr,s2t_ptr+s2t_col*s2t_row,&s2t[0][0]); - prog = 6; - return true; - } - } - } - tipl::progress prog_("running normalization"); - prog = 0; - bool terminated = false; - auto lambda = [this,output_file_name,&terminated]() - { - tipl::transformation_matrix T; - - auto It = template_I; - auto It2 = template_I2; - - if(dir.index_name[0] == "image" && // not FIB file - tipl::io::gz_nifti(t1w_template_file_name.c_str(),std::ios::in).to_space(It,template_to_mni)) // not FIB file, use t1w as template - { - tipl::out() << "using structure image for normalization" << std::endl; - It2.clear(); - } - - tipl::image<3> Is(dir.fa[0],dim); - tipl::image<3> Is2; - - { - size_t iso_index = get_name_index("iso"); - if(slices.size() == iso_index) - iso_index = get_name_index("md"); - if(slices.size() != iso_index) - Is2 = slices[iso_index]->get_image(); - } - bool no_iso = Is2.empty() || It2.empty(); - - tipl::filter::gaussian(Is); - tipl::filter::gaussian(Is); - if(!no_iso) - { - tipl::filter::gaussian(Is2); - tipl::filter::gaussian(Is2); - } - - prog = 1; - if(!has_manual_atlas) - linear_with_mi(It,template_vs,Is,vs,T,tipl::reg::affine,terminated); - else - T = tipl::transformation_matrix(manual_template_T,It.shape(),template_vs,Is.shape(),vs); - - if(terminated) - return; - prog = 2; - tipl::image<3> Iss(It.shape()); - tipl::resample(Is,Iss,T); - tipl::image<3> Iss2(It2.shape()); - if(!no_iso) - tipl::resample(Is2,Iss2,T); - prog = 3; - - if(dir.index_name[0] == "image") - { - tipl::out() << "matching t1w t2w contrast" << std::endl; - tipl::image<3> t2w(It.shape()); - if(tipl::io::gz_nifti(t2w_template_file_name.c_str(),std::ios::in).to_space(t2w,template_to_mni)) - { - std::vector X(It.size()*3); - tipl::par_for(It.size(),[&](size_t pos) - { - if(Iss[pos] == 0.0f) - return; - size_t i = pos; - pos = pos+pos+pos; - X[pos] = 1; - X[pos+1] = It[i]; - X[pos+2] = t2w[i]; - }); - tipl::multiple_regression m; - if(m.set_variables(X.begin(),3,It.size())) - { - float b[3] = {0.0f,0.0f,0.0f}; - m.regress(Iss.begin(),b); - tipl::out() << "image=" << b[0] << " + " << b[1] << " × t1w + " << b[2] << " × t2w "; - tipl::par_for(It.size(),[&](size_t pos) - { - if(Iss[pos] == 0.0f) - return; - It[pos] = b[0] + b[1]*It[pos] + b[2]*t2w[pos]; - if(It[pos] < 0.0f) - It[pos] = 0.0f; - }); - } - } - } - tipl::image<3,tipl::vector<3> > dis,inv_dis; - tipl::reg::cdm_pre(It,It2,Iss,Iss2); - - chen_cdm_common(It,It2,Iss,Iss2,dis,inv_dis,terminated); - - tipl::displacement_to_mapping(dis,t2s,T); - tipl::compose_mapping(Is,t2s,Iss); - tipl::out() << "R2:" << tipl::correlation(Iss.begin(),Iss.end(),It.begin()) << std::endl; - - s2t.resize(dim); - tipl::inv_displacement_to_mapping(inv_dis,s2t,T); - - if(terminated) - return; - tipl::io::gz_mat_write out(output_file_name.c_str()); - if(out) - { - out.write("to_dim",dis.shape()); - out.write("to_vs",template_vs); - out.write("from_dim",dim); - out.write("from_vs",vs); - out.write("steps",steps); - prog = 4; - tipl::out() << "calculating template to subject warp field"; - out.write("to2from",&t2s[0][0],3,t2s.size()); - prog = 5; - tipl::out() << "calculating subject to template warp field"; - out.write("from2to",&s2t[0][0],3,s2t.size()); - } - - prog = 6; - }; - - if(background) - { - std::thread t(lambda); - while(prog_(prog,6)) - std::this_thread::yield(); - if(prog_.aborted()) - { - error_msg = "aborted."; - terminated = true; - } - t.join(); - return !prog_.aborted(); - } - else - { - tipl::out() << "Subject normalization to MNI space." << std::endl; - lambda(); - } - return true; -} bool fib_data::map_to_mni(bool background) { - if(use_chen_normalization) - return map_to_mni_chen(background); if(!load_template()) return false; if(is_mni && template_id == matched_template_id) diff --git a/libs/tracking/fib_data.hpp b/libs/tracking/fib_data.hpp index be2d8fe8..0d27a339 100644 --- a/libs/tracking/fib_data.hpp +++ b/libs/tracking/fib_data.hpp @@ -226,7 +226,6 @@ class fib_data bool is_histology = true; bool is_mni = false; bool is_be = false; - bool use_chen_normalization = false; size_t matched_template_id = 0; bool trackable = true; float min_length(void) const @@ -309,7 +308,6 @@ class fib_data public: std::string get_mapping_file_name(void) const; bool map_to_mni(bool background = true); - bool map_to_mni_chen(bool background = true); void temp2sub(std::vector >&tracts) const; void temp2sub(tipl::vector<3>& pos) const; void sub2temp(tipl::vector<3>& pos); diff --git a/reg.hpp b/reg.hpp index 7b9061fb..31b0c479 100644 --- a/reg.hpp +++ b/reg.hpp @@ -181,18 +181,6 @@ namespace tipl { namespace reg { -template -inline void cdm_pre(image_type& It,image_type& It2,image_type& Is,image_type& Is2) -{ - if(!It.empty()) - tipl::normalize_upper_lower2(It,It,255.99f); - if(!Is.empty()) - tipl::normalize_upper_lower2(Is,Is,255.99f); - if(!It2.empty()) - tipl::normalize_upper_lower2(It2,It2,255.99f); - if(!Is2.empty()) - tipl::normalize_upper_lower2(Is2,Is2,255.99f); -} } } From 0711ca5762d11d28d81cc7ba8e1bdfd0e7971b40 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Wed, 4 Feb 2026 16:01:41 -0500 Subject: [PATCH 10/11] add too much debugging and add random restarts to linear reg --- auto_track.cpp | 277 +++++++++++++++++++++++++++- cmd/reg.cpp | 66 ++++++- libs/dsi/gqi_mni_reconstruction.hpp | 45 +++++ libs/tracking/fib_data.cpp | 65 +++++++ libs/tracking/tract_model.cpp | 131 ++++++++++++- reg.hpp | 1 + 6 files changed, 564 insertions(+), 21 deletions(-) diff --git a/auto_track.cpp b/auto_track.cpp index bbf34bb2..28bcb3c6 100644 --- a/auto_track.cpp +++ b/auto_track.cpp @@ -7,6 +7,9 @@ #include "fib_data.hpp" #include "libs/tracking/tracking_thread.hpp" #include +#include +#include +#include extern std::vector fa_template_list; auto_track::auto_track(QWidget *parent) : QMainWindow(parent), @@ -106,6 +109,10 @@ int trk_post(tipl::program_option& po, std::string run_auto_track(tipl::program_option& po,const std::vector& file_list,int& prog) { const bool chen_mode = po.get("chen_mode",0); + const bool debug_shape = po.get("debug_shape",0); + const bool debug_streamlines = po.get("debug_streamlines",0) || debug_shape; + const bool delete_repeats = po.has("delete_repeat") || po.get("action") == "atk" || po.has("track_id"); + const float delete_repeat_distance = po.get("delete_repeat",float(0.5f)); std::string trk_format = po.get("trk_format","tt.gz"); std::string tolerance_string = po.get("tolerance","22,26,30"); float yield_rate = po.get("yield_rate",0.00001f); @@ -214,9 +221,45 @@ std::string run_auto_track(tipl::program_option& po,const std::vector std::shared_ptr handle; std::string fib_base = QFileInfo(fib_file_name.c_str()).baseName().toStdString(); bool save_atlas_track = po.has("template_tracks_in_native_space"); + bool save_template_images = po.has("templates_images_in_native_space"); std::string atlas_track_option; if(save_atlas_track) atlas_track_option = po.get("template_tracks_in_native_space",std::string()); + std::string template_images_option; + if(save_template_images) + template_images_option = po.get("templates_images_in_native_space",std::string()); + + if(save_template_images) + { + if(!handle.get()) + { + handle = std::make_shared(); + if(!handle->load_from_file(fib_file_name.c_str())) + return handle->error_msg; + set_template(handle,po); + } + if(!handle->map_to_mni(tipl::show_prog)) + return handle->error_msg + " at " + fib_file_name; + if(handle->s2t.empty()) + return std::string("no subject-to-template mapping at ") + fib_file_name; + std::string output_base = template_images_option; + if(output_base.empty()) + output_base = dir + "/" + fib_base + ".template"; + else if(std::filesystem::is_directory(output_base)) + output_base = (std::filesystem::path(output_base) / (fib_base + ".template")).string(); + tipl::out() << "saving template images in native space to " << output_base; + { + auto warped_qa = tipl::compose_mapping(handle->template_I,handle->s2t); + if(!(tipl::io::gz_nifti(output_base + ".qa.nii.gz",std::ios::out) << handle->bind(warped_qa))) + return std::string("failed to save warped template qa to ") + output_base + ".qa.nii.gz"; + } + if(!handle->template_I2.empty()) + { + auto warped_iso = tipl::compose_mapping(handle->template_I2,handle->s2t); + if(!(tipl::io::gz_nifti(output_base + ".iso.nii.gz",std::ios::out) << handle->bind(warped_iso))) + return std::string("failed to save warped template iso to ") + output_base + ".iso.nii.gz"; + } + } if(save_atlas_track) { @@ -234,6 +277,20 @@ std::string run_auto_track(tipl::program_option& po,const std::vector atlas_track_name = dir + "/" + fib_base + ".atlas.tt.gz"; else if(std::filesystem::is_directory(atlas_track_name)) atlas_track_name = (std::filesystem::path(atlas_track_name) / (fib_base + ".atlas.tt.gz")).string(); + if(!tipl::ends_with(atlas_track_name,"tt.gz") && + !tipl::ends_with(atlas_track_name,".trk") && + !tipl::ends_with(atlas_track_name,".trk.gz") && + !tipl::ends_with(atlas_track_name,".tck") && + !tipl::ends_with(atlas_track_name,".txt") && + !tipl::ends_with(atlas_track_name,".mat") && + !tipl::ends_with(atlas_track_name,".nii") && + !tipl::ends_with(atlas_track_name,".nii.gz")) + { + std::string track_ext = trk_format; + if(!track_ext.empty() && track_ext[0] != '.') + track_ext = "." + track_ext; + atlas_track_name += track_ext.empty() ? ".tt.gz" : track_ext; + } tipl::out() << "saving template tracks in native space to " << atlas_track_name; if(!handle->track_atlas || !handle->track_atlas->save_tracts_to_file(atlas_track_name)) return std::string("failed to save atlas track to ") + atlas_track_name; @@ -294,9 +351,143 @@ std::string run_auto_track(tipl::program_option& po,const std::vector if(!overwrite && std::filesystem::exists(trk_file_name)) tract_model->load_tracts_from_file(trk_file_name,handle.get()); + auto format_bids_value = [](float value) + { + std::ostringstream out; + out << std::fixed << std::setprecision(3) << value; + std::string text = out.str(); + while(text.size() > 1 && text.back() == '0') + text.pop_back(); + if(!text.empty() && text.back() == '.') + text.pop_back(); + for(char& c : text) + if(c == '.') + c = 'p'; + if(!text.empty() && text.front() == '-') + text.front() = 'm'; + return text; + }; + + auto format_numeric_value = [](float value) + { + std::ostringstream out; + out << std::fixed << std::setprecision(3) << value; + std::string text = out.str(); + while(text.size() > 1 && text.back() == '0') + text.pop_back(); + if(!text.empty() && text.back() == '.') + text.pop_back(); + return text; + }; + + auto make_debug_track_name = [&](const std::string& debug_base,const std::string& entity)->std::string + { + std::string name = debug_base + "_desc-debug"; + if(!entity.empty()) + name += "_" + entity; + if(!tipl::ends_with(name,{".tt.gz",".trk.gz",".trk",".tck",".txt",".mat"})) + { + std::string track_ext = trk_format; + if(!track_ext.empty() && track_ext[0] != '.') + track_ext = "." + track_ext; + name += track_ext.empty() ? ".tt.gz" : track_ext; + } + return name; + }; + auto add_prefix_to_filename = [](const std::filesystem::path& path, + const std::string& prefix)->std::string + { + std::filesystem::path output(path); + output.replace_filename(prefix + output.filename().string()); + return output.string(); + }; + + auto save_debug_tracks = [&](TractModel* model, + const std::string& debug_base, + const std::string& entity, + std::shared_ptr& output_model)->std::string + { + if(delete_repeats) + { + output_model = std::make_shared(*model); + output_model->delete_repeated(delete_repeat_distance); + } + else + { + output_model = std::shared_ptr(model,[](TractModel*){}); + } + std::string output_name = make_debug_track_name(debug_base,entity); + if(!output_model->save_tracts_to_file(output_name)) + { + tipl::warning() << "failed to save debug streamlines to " << output_name; + return std::string(); + } + if(!output_model->save_tracts_in_template_space(handle, + add_prefix_to_filename(output_name,"T_"))) + tipl::warning() << "failed to save debug streamlines in template space for " << output_name; + return output_name; + }; + + std::ofstream debug_shape_out; + std::vector debug_shape_titles; + bool debug_shape_header_written = false; + const std::string debug_shape_table_name = trk_base + "_desc-debugshape.tsv"; + auto write_debug_shape_row = [&](TractModel* model, + const std::string& streamline_file, + float tolerance_value, + int tip_iteration, + float min_length_value, + float max_length_value, + int official_flag) + { + if(!debug_shape || streamline_file.empty()) + return; + std::vector debug_shape_data; + std::vector current_titles; + model->get_quantitative_info(handle,current_titles,debug_shape_data); + if(!debug_shape_header_written) + { + debug_shape_out.open(debug_shape_table_name); + if(!debug_shape_out) + { + tipl::warning() << "failed to save debug shape stats to " << debug_shape_table_name; + return; + } + debug_shape_out << "streamline_file\ttolerance\ttip_iteration\tmin_length_mm\tmax_length_mm\tofficial"; + debug_shape_titles = current_titles; + for(const auto& title : debug_shape_titles) + debug_shape_out << "\t" << title; + debug_shape_out << "\n"; + debug_shape_header_written = true; + } + if(!debug_shape_out) + return; + debug_shape_out << streamline_file << "\t" + << format_numeric_value(tolerance_value) << "\t" + << tip_iteration << "\t" + << format_numeric_value(min_length_value) << "\t" + << format_numeric_value(max_length_value) << "\t" + << official_flag; + size_t count = std::min(debug_shape_titles.size(),debug_shape_data.size()); + for(size_t index = 0;index < count;++index) + { + if(std::isnan(debug_shape_data[index])) + debug_shape_out << "\t"; + else + debug_shape_out << "\t" << std::to_string(debug_shape_data[index]); + } + debug_shape_out << "\n"; + }; + + float official_tolerance = std::numeric_limits::quiet_NaN(); + int official_tip_iteration = 0; + float official_min_length = std::numeric_limits::quiet_NaN(); + float official_max_length = std::numeric_limits::quiet_NaN(); + bool official_set = false; + // each iteration increases tolerance - for(size_t tracking_iteration = 0;tracking_iteration < tolerance.size() && - !tract_model->get_visible_track_count();++tracking_iteration) + for(size_t tracking_iteration = 0;tracking_iteration < tolerance.size() && + (debug_streamlines || !tract_model->get_visible_track_count());++tracking_iteration) { ThreadData thread(handle); { @@ -398,17 +589,82 @@ std::string run_auto_track(tipl::program_option& po,const std::vector } - if(no_result) + if(no_result && !debug_streamlines) continue; + + std::shared_ptr iteration_model = tract_model; + if(debug_streamlines) + iteration_model = std::make_shared(handle); + // fetch both front and back buffer - thread.fetchTracks(tract_model.get()); - thread.fetchTracks(tract_model.get()); + thread.fetchTracks(iteration_model.get()); + thread.fetchTracks(iteration_model.get()); + if(debug_streamlines) + { + std::string debug_base = trk_base; + std::string tolerance_entity = "tolerance-" + format_bids_value(tolerance[tracking_iteration]); + std::shared_ptr output_model; + save_debug_tracks(iteration_model.get(),debug_base,tolerance_entity,output_model); + } + if(thread.param.step_size != 0.0f) - tract_model->resample(1.0f); - tract_model->trim(thread.param.tip_iteration); + { + iteration_model->resample(1.0f); + if(debug_streamlines) + { + std::string debug_base = trk_base; + std::string tolerance_entity = "stage-resample_tolerance-" + + format_bids_value(tolerance[tracking_iteration]); + std::shared_ptr output_model; + std::string saved_name = save_debug_tracks(iteration_model.get(),debug_base,tolerance_entity,output_model); + write_debug_shape_row(output_model.get(),saved_name,tolerance[tracking_iteration],0, + thread.param.min_length,thread.param.max_length,0); + } + } + if(debug_streamlines) + { + std::string debug_base = trk_base; + std::string tolerance_entity = "tip-0_tolerance-" + + format_bids_value(tolerance[tracking_iteration]); + std::shared_ptr output_model; + std::string saved_name = save_debug_tracks(iteration_model.get(),debug_base,tolerance_entity,output_model); + write_debug_shape_row(output_model.get(),saved_name,tolerance[tracking_iteration],0, + thread.param.min_length,thread.param.max_length,0); + } + iteration_model->trim(thread.param.tip_iteration); // if trim removes too many tract, undo to at least get the smallest possible bundle. - if(thread.param.tip_iteration && tract_model->get_visible_track_count() == 0) - tract_model->undo(); + if(thread.param.tip_iteration && iteration_model->get_visible_track_count() == 0) + iteration_model->undo(); + if(debug_streamlines) + { + std::string debug_base = trk_base; + std::string tolerance_entity = "tip-" + std::to_string(thread.param.tip_iteration) + + "_tolerance-" + format_bids_value(tolerance[tracking_iteration]); + std::shared_ptr output_model; + std::string saved_name = save_debug_tracks(iteration_model.get(),debug_base,tolerance_entity,output_model); + write_debug_shape_row(output_model.get(),saved_name,tolerance[tracking_iteration],thread.param.tip_iteration, + thread.param.min_length,thread.param.max_length,0); + } + + if(debug_streamlines && + tract_model->get_visible_track_count() == 0 && + iteration_model->get_visible_track_count() != 0) + { + tract_model = iteration_model; + official_tolerance = tolerance[tracking_iteration]; + official_tip_iteration = thread.param.tip_iteration; + official_min_length = thread.param.min_length; + official_max_length = thread.param.max_length; + official_set = true; + } + if(!official_set && iteration_model->get_visible_track_count() != 0) + { + official_tolerance = tolerance[tracking_iteration]; + official_tip_iteration = thread.param.tip_iteration; + official_min_length = thread.param.min_length; + official_max_length = thread.param.max_length; + official_set = true; + } } if(tract_model->get_visible_track_count() == 0) @@ -422,6 +678,9 @@ std::string run_auto_track(tipl::program_option& po,const std::vector po.set("output",trk_file_name); if(trk_post(po,handle,tract_model,trk_file_name,true)) return std::string("terminated due to error"); + if(debug_shape) + write_debug_shape_row(tract_model.get(),trk_file_name,official_tolerance,official_tip_iteration, + official_min_length,official_max_length,1); } } } diff --git a/cmd/reg.cpp b/cmd/reg.cpp index 40e599c6..f6ee1cdb 100644 --- a/cmd/reg.cpp +++ b/cmd/reg.cpp @@ -1,5 +1,7 @@ #include #include +#include +#include #include "img.hpp" #include "reg.hpp" #include "tract_model.hpp" @@ -280,8 +282,68 @@ float dual_reg::linear_reg(bool& terminated) tipl::progress prog("linear registration"); float cost = 0.0f; if(!skip_linear) - cost = tipl::reg::linear(tipl::reg::make_list(It),Itvs,tipl::reg::make_list(I),Ivs, - arg,reg_type,terminated,bound,cost_type,use_cuda && has_cuda); + { + auto best_arg = arg; + float best_cost = std::numeric_limits::max(); + auto run_linear = [&](tipl::affine_transform& start_arg) + { + auto local_arg = start_arg; + float local_cost = tipl::reg::linear( + tipl::reg::make_list(It),Itvs, + tipl::reg::make_list(I),Ivs, + local_arg,reg_type,terminated,bound,cost_type,use_cuda && has_cuda); + if(!terminated && local_cost < best_cost) + { + best_cost = local_cost; + best_arg = local_arg; + } + }; + + tipl::out() << "linear registration restarts planned: " << linear_restarts; + run_linear(arg); + if(linear_restarts > 1) + { + std::mt19937 rng(42); + std::uniform_real_distribution unit01(0.0f,1.0f); + auto image_extent = [&](unsigned int axis) + { + float from_mm = It[0].shape()[axis]*Itvs[axis]; + float to_mm = I[0].shape()[axis]*Ivs[axis]; + return std::max(from_mm,to_mm)*0.5f; + }; + for(size_t restart = 1;restart < linear_restarts;++restart) + { + if(terminated) + { + tipl::out() << "linear registration restarts aborted"; + break; + } + tipl::out() << "linear registration restart " << (restart + 1) + << "/" << linear_restarts; + tipl::affine_transform seed; + seed.clear(); + for(unsigned int axis = 0;axis < dimension;++axis) + { + float range = image_extent(axis); + float tmin = range*bound[axis][1]; + float tmax = range*bound[axis][0]; + seed.translocation[axis] = tmin + (tmax-tmin)*unit01(rng); + float smin = bound[axis][5]; + float smax = bound[axis][4]; + seed.scaling[axis] = smin + (smax-smin)*unit01(rng); + float rmin = bound[axis][3]*3.14159265358979323846f; + float rmax = bound[axis][2]*3.14159265358979323846f; + seed.rotation[axis] = rmin + (rmax-rmin)*unit01(rng); + float amin = bound[axis][7]; + float amax = bound[axis][6]; + seed.affine[axis] = amin + (amax-amin)*unit01(rng); + } + run_linear(seed); + } + } + arg = best_arg; + cost = best_cost; + } calculate_linear_r(); diff --git a/libs/dsi/gqi_mni_reconstruction.hpp b/libs/dsi/gqi_mni_reconstruction.hpp index fe73f10f..22d406d9 100644 --- a/libs/dsi/gqi_mni_reconstruction.hpp +++ b/libs/dsi/gqi_mni_reconstruction.hpp @@ -2,6 +2,7 @@ #define MNI_RECONSTRUCTION_HPP #include #include +#include #include "basic_voxel.hpp" #include "basic_process.hpp" #include "gqi_process.hpp" @@ -72,6 +73,32 @@ class DWINormalization : public BaseProcess } else { + auto nonzero_voxels = [](const auto& image) + { + size_t count = 0; + for(size_t i = 0;i < image.size();++i) + if(image[i]) + ++count; + return count; + }; + if(!reg.I[0].empty() && !reg.It[0].empty()) + { + auto subject_voxels = nonzero_voxels(reg.I[0]); + auto template_voxels = nonzero_voxels(reg.It[0]); + double subject_volume = double(subject_voxels) * reg.Ivs[0] * reg.Ivs[1] * reg.Ivs[2]; + double template_volume = double(template_voxels) * reg.Itvs[0] * reg.Itvs[1] * reg.Itvs[2]; + if(template_volume > 0.0) + tipl::out() << "brain volume ratio (subject/template): " + << (subject_volume / template_volume); + } + reg.linear_restarts = 6; + // Use a wider affine search (especially scaling) for lifespan differences. + static const float lifespan_affine_bound[3][8] = { + {1.0f,-1.0f, 0.5f,-0.5f, 3.0f,0.3f, 0.25f,-0.25f}, + {1.0f,-1.0f, 0.4f,-0.4f, 3.0f,0.3f, 0.25f,-0.25f}, + {1.0f,-1.0f, 0.4f,-0.4f, 3.0f,0.3f, 0.25f,-0.25f} + }; + reg.bound = lifespan_affine_bound; tipl::run("linear registration",[&](void) { reg.linear_reg(tipl::prog_aborted); @@ -98,6 +125,24 @@ class DWINormalization : public BaseProcess tipl::out() << "nonlinear R2: " << voxel.R2 << std::endl; if(voxel.R2 < 0.3f) tipl::warning() << "poor registration found in nonlinear registration. Please check image quality or image orientation"; + if(!reg.t2f_dis.empty()) + { + float max_displacement_mm = 0.0f; + for(size_t i = 0;i < reg.t2f_dis.size();++i) + { + const auto& v = reg.t2f_dis[i]; + float dx = v[0]*reg.Itvs[0]; + float dy = v[1]*reg.Itvs[1]; + float dz = v[2]*reg.Itvs[2]; + float dis = std::sqrt(dx*dx + dy*dy + dz*dz); + if(dis > max_displacement_mm) + max_displacement_mm = dis; + } + if(max_displacement_mm > 10.0f) + tipl::warning() << "large nonlinear displacement detected (max " + << max_displacement_mm + << " mm). Please verify template-to-native alignment."; + } auto new_ItR = reg.ItR; auto new_Its = reg.Its; diff --git a/libs/tracking/fib_data.cpp b/libs/tracking/fib_data.cpp index 3fcb5b3d..3ffc96fb 100644 --- a/libs/tracking/fib_data.cpp +++ b/libs/tracking/fib_data.cpp @@ -1,5 +1,7 @@ #include #include +#include +#include #include #include #include @@ -1788,6 +1790,67 @@ bool fib_data::load_track_atlas(bool symmetric) cluster.insert(cluster.end(),new_cluster.begin(),new_cluster.end()); + auto log_tract_stats = [&](const char* label) + { + const auto& cur_tracts = track_atlas->get_tracts(); + size_t empty_count = 0; + size_t invalid_count = 0; + size_t non_finite_count = 0; + double min_len = std::numeric_limits::max(); + double max_len = 0.0; + size_t len_count = 0; + size_t sample_count = 0; + for(size_t i = 0;i < cur_tracts.size();++i) + { + const auto& tract = cur_tracts[i]; + if(tract.empty()) + ++empty_count; + if(tract.size() < 3 || (tract.size() % 3) != 0) + { + ++invalid_count; + } + else + { + bool has_non_finite = false; + for(float v : tract) + if(!std::isfinite(v)) + { + has_non_finite = true; + break; + } + if(has_non_finite) + ++non_finite_count; + double length = track_atlas->get_tract_length_in_mm(i); + if(length > 0.0) + { + min_len = std::min(min_len,length); + max_len = std::max(max_len,length); + ++len_count; + } + } + if((tract.empty() || (tract.size() < 3 || (tract.size() % 3) != 0)) && sample_count < 3) + { + auto c = (i < cluster.size() ? cluster[i] : tractography_name_list.size()); + std::string name = (c < tractography_name_list.size() ? tractography_name_list[c] : "unknown"); + tipl::warning() << label << " invalid tract index=" << i + << " size=" << tract.size() + << " name=" << name; + ++sample_count; + } + } + if(min_len == std::numeric_limits::max()) + min_len = 0.0; + tipl::out() << label << " tracts total=" << cur_tracts.size() + << " empty=" << empty_count + << " invalid=" << invalid_count + << " non_finite=" << non_finite_count + << " length_min=" << min_len + << " length_max=" << max_len + << " length_count=" << len_count; + }; + + log_tract_stats("pre-warp"); + // get distance scaling auto& s2t = get_sub2temp_mapping(); if(s2t.empty()) @@ -1796,6 +1859,8 @@ bool fib_data::load_track_atlas(bool symmetric) // warp tractography atlas to subject space temp2sub(track_atlas->get_tracts()); + log_tract_stats("post-warp"); + auto& tract_data = track_atlas->get_tracts(); // get min max length std::vector min_length(tractography_name_list.size()),max_length(tractography_name_list.size()); diff --git a/libs/tracking/tract_model.cpp b/libs/tracking/tract_model.cpp index 5f27763b..0dd8f537 100644 --- a/libs/tracking/tract_model.cpp +++ b/libs/tracking/tract_model.cpp @@ -20,6 +20,9 @@ #include "tracking_method.hpp" #include "reg.hpp" #include +#include +#include +#include void prepare_idx(const std::string& file_name,std::shared_ptr in); void save_idx(const std::string& file_name,std::shared_ptr in); @@ -78,9 +81,18 @@ class TinyTrack{ const std::vector& color) { tipl::progress prog0("saving " + file_name); + tipl::out() << "saving tt: file=" << file_name + << " tracts=" << tract_data.size() + << " cluster=" << cluster.size() + << " color=" << color.size(); + if(tract_data.empty()) + tipl::warning() << "no tracts to save for " << file_name; tipl::io::gz_mat_write out(file_name); if (!out) + { + tipl::warning() << "failed to open " << file_name << " for writing"; return false; + } out.write("dimension",geo); out.write("voxel_size",vs); out.write("trans_to_mni",trans_to_mni); @@ -92,16 +104,19 @@ class TinyTrack{ std::vector buf_size(track32.size()); { - tipl::progress prog1("compressing trajectories"); - size_t p = 0; - tipl::adaptive_par_for(track32.size(),[&](size_t i) + auto compress_track = [&](size_t i)->bool { - prog1(p++,tract_data.size()); auto& t32 = track32[i]; + if(tract_data[i].size() < 3 || (tract_data[i].size() % 3) != 0) + return false; t32.resize(tract_data[i].size()); // all coordinates multiply by 32 and convert to integer for(size_t j = 0;j < t32.size();j++) + { + if(!std::isfinite(tract_data[i][j])) + return false; t32[j] = int(std::round(std::ldexp(tract_data[i][j],5))); + } // Calculate coordinate displacement, skipping the first coordinate for(size_t j = t32.size()-1;j >= 3;j--) t32[j] -= t32[j-3]; @@ -146,9 +161,60 @@ class TinyTrack{ new_t32.swap(t32); } buf_size[i] = sizeof(tract_header)+t32.size()-3; + return true; + }; + + tipl::progress prog1("compressing trajectories"); + std::atomic invalid_tract(false); + std::atomic invalid_index(std::numeric_limits::max()); + std::atomic invalid_size(0); + size_t p = 0; + tipl::adaptive_par_for(track32.size(),[&](size_t i) + { + if(invalid_tract.load(std::memory_order_relaxed)) + return; + prog1(p++,tract_data.size()); + if(!compress_track(i)) + { + invalid_tract.store(true,std::memory_order_relaxed); + invalid_index.store(i,std::memory_order_relaxed); + invalid_size.store(tract_data[i].size(),std::memory_order_relaxed); + } }); - if(prog1.aborted()) + if(invalid_tract.load(std::memory_order_relaxed)) + { + size_t bad_index = invalid_index.load(std::memory_order_relaxed); + size_t bad_size = invalid_size.load(std::memory_order_relaxed); + tipl::warning() << "invalid tract data while saving " << file_name + << " at index " << bad_index + << " size " << bad_size + << " (expect >=3 and multiple of 3, finite values)"; return false; + } + if(prog1.aborted()) + { + tipl::warning() << "compression aborted while saving " << file_name + << ", retrying single-threaded"; + tipl::progress prog1_single("compressing trajectories (single-thread)"); + for(size_t i = 0;i < track32.size();++i) + { + prog1_single(i,track32.size()); + if(!compress_track(i)) + { + tipl::warning() << "invalid tract data while saving " << file_name + << " at index " << i + << " size " << tract_data[i].size() + << " (expect >=3 and multiple of 3, finite values)"; + return false; + } + } + if(prog1_single.aborted()) + { + tipl::warning() << "compression aborted while saving " << file_name + << " (single-thread)"; + return false; + } + } } { for(size_t block = 0,cur_track_block = 0;prog0(cur_track_block,track32.size());++block) @@ -180,6 +246,9 @@ class TinyTrack{ out[j] = char(t32[j]); }); + tipl::out() << "writing track block " << block + << " count=" << pos.size() + << " bytes=" << total_size; if(block == 0) out.write("track",&out_buf[0],total_size,1); else @@ -187,7 +256,12 @@ class TinyTrack{ cur_track_block += pos.size(); } } - return !prog0.aborted(); + if(prog0.aborted()) + { + tipl::warning() << "save aborted while writing " << file_name; + return false; + } + return true; } static bool load_from_file(const std::string& file_name, std::vector >& tract_data, @@ -900,9 +974,46 @@ bool TractModel::save_tracts_to_file(const std::string& file_name) tipl::out() << "trans:" << trans_to_mni; if(tipl::ends_with(file_name,"tt.gz")) { - return TinyTrack::save_to_file(file_name,geo,vs,trans_to_mni, + const int save_retries = 3; + const int save_timeout_sec = 10; + bool saved_ok = false; + auto save_start = std::chrono::steady_clock::now(); + for(int attempt = 1;attempt <= save_retries;++attempt) + { + if(TinyTrack::save_to_file(file_name,geo,vs,trans_to_mni, tract_data,std::vector(tract_cluster.begin(),tract_cluster.end()),report,parameter_id, - std::vector{get_cluster_color(tract_color)}); + std::vector{get_cluster_color(tract_color)})) + { + saved_ok = true; + break; + } + if(save_timeout_sec > 0) + { + auto elapsed = std::chrono::duration_cast( + std::chrono::steady_clock::now()-save_start).count(); + if(elapsed >= save_timeout_sec) + break; + } + if(attempt < save_retries) + { + tipl::warning() << "failed to save " << file_name << ", retrying (" + << (attempt+1) << "/" << save_retries << ")"; + std::this_thread::sleep_for(std::chrono::seconds(1)); + } + } + if(saved_ok) + return true; + std::filesystem::path fallback_path(file_name); + if(fallback_path.extension() == ".gz") + { + fallback_path.replace_extension(); // drop .gz -> .tt + tipl::warning() << "failed to save " << file_name + << ", falling back to uncompressed " << fallback_path.string(); + return TinyTrack::save_to_file(fallback_path.string(),geo,vs,trans_to_mni, + tract_data,std::vector(tract_cluster.begin(),tract_cluster.end()),report,parameter_id, + std::vector{get_cluster_color(tract_color)}); + } + return false; } if(tipl::ends_with(file_name,".trk") || tipl::ends_with(file_name,".trk.gz")) { @@ -1153,7 +1264,7 @@ std::string TractModel::get_obj(unsigned int& coordinate_count, //--------------------------------------------------------------------------- bool TractModel::save_all(const std::string& file_name, const std::vector >& all) -{ +{ if(all.empty()) return false; tipl::progress prog("save " + file_name); @@ -1594,7 +1705,7 @@ bool TractModel::select_tracts(const std::vector& tracts_to_select } //--------------------------------------------------------------------------- bool TractModel::delete_repeated(float d) -{ +{ std::vector > x_reg; std::vector track_location1,track_location2; { diff --git a/reg.hpp b/reg.hpp index 31b0c479..747b5440 100644 --- a/reg.hpp +++ b/reg.hpp @@ -90,6 +90,7 @@ struct dual_reg{ bool skip_nonlinear = false; bool match_fov = true; bool masked_r = true; + size_t linear_restarts = 6; public: dual_reg(void):modality_names(max_modality),I(max_modality),J(max_modality),It(max_modality),r(max_modality) { From 04dbbafdb4cfced496ea842d2cce5f0f7a2b8098 Mon Sep 17 00:00:00 2001 From: Matt Cieslak Date: Wed, 4 Feb 2026 17:43:45 -0500 Subject: [PATCH 11/11] remove linear_mi, refactor autotrack --debug --- auto_track.cpp | 144 +++++++++++++++++++++++-------------------------- reg.hpp | 47 ---------------- 2 files changed, 66 insertions(+), 125 deletions(-) diff --git a/auto_track.cpp b/auto_track.cpp index 28bcb3c6..0888f928 100644 --- a/auto_track.cpp +++ b/auto_track.cpp @@ -10,6 +10,8 @@ #include #include #include +#include +#include extern std::vector fa_template_list; auto_track::auto_track(QWidget *parent) : QMainWindow(parent), @@ -108,9 +110,58 @@ int trk_post(tipl::program_option& po, std::string tract_file_name,bool output_track); std::string run_auto_track(tipl::program_option& po,const std::vector& file_list,int& prog) { - const bool chen_mode = po.get("chen_mode",0); - const bool debug_shape = po.get("debug_shape",0); - const bool debug_streamlines = po.get("debug_streamlines",0) || debug_shape; + bool debug_streamlines = false; + bool debug_shape = false; + bool debug_template_tracks = false; + bool debug_template_images = false; + auto enable_all_debug = [&]() + { + debug_streamlines = true; + debug_shape = true; + debug_template_tracks = true; + debug_template_images = true; + }; + if(po.has("debug")) + { + std::string debug_option = po.get("debug",std::string()); + auto items = tipl::split(debug_option,','); + bool has_item = false; + for(auto& raw_item : items) + { + auto trim_item = [](std::string text) + { + auto is_space = [](unsigned char c){return std::isspace(c);}; + text.erase(text.begin(),std::find_if_not(text.begin(),text.end(),is_space)); + text.erase(std::find_if_not(text.rbegin(),text.rend(),is_space).base(),text.end()); + std::transform(text.begin(),text.end(),text.begin(), + [](unsigned char c){return std::tolower(c);}); + return text; + }; + std::string item = trim_item(raw_item); + if(item.empty()) + continue; + has_item = true; + if(item == "all" || item == "true" || item == "1") + { + enable_all_debug(); + break; + } + if(item == "streamlines") + debug_streamlines = true; + else if(item == "shape") + debug_shape = true; + else if(item == "template_images") + debug_template_images = true; + else if(item == "template_streamlines") + debug_template_tracks = true; + else + tipl::warning() << "unknown debug option: " << item; + } + if(!has_item) + enable_all_debug(); + } + if(debug_shape) + debug_streamlines = true; const bool delete_repeats = po.has("delete_repeat") || po.get("action") == "atk" || po.has("track_id"); const float delete_repeat_distance = po.get("delete_repeat",float(0.5f)); std::string trk_format = po.get("trk_format","tt.gz"); @@ -119,13 +170,6 @@ std::string run_auto_track(tipl::program_option& po,const std::vector size_t yield_check_count = 10.0f/yield_rate; bool overwrite = po.get("overwrite",0); uint32_t thread_count = tipl::max_thread_count; - if(chen_mode) - { - po.set("track_voxel_ratio",po.get("track_voxel_ratio",2.0f)); - po.set("tip_iteration",po.get("tip_iteration",4)); - po.set("use_roi",po.get("use_roi",0)); - po.set("check_ending",po.get("check_ending",1)); - } std::vector tolerance; { if(!po.has("export")) @@ -220,14 +264,8 @@ std::string run_auto_track(tipl::program_option& po,const std::vector tipl::out() << "processing " << cur_file_base_name << std::endl; std::shared_ptr handle; std::string fib_base = QFileInfo(fib_file_name.c_str()).baseName().toStdString(); - bool save_atlas_track = po.has("template_tracks_in_native_space"); - bool save_template_images = po.has("templates_images_in_native_space"); - std::string atlas_track_option; - if(save_atlas_track) - atlas_track_option = po.get("template_tracks_in_native_space",std::string()); - std::string template_images_option; - if(save_template_images) - template_images_option = po.get("templates_images_in_native_space",std::string()); + bool save_atlas_track = debug_template_tracks; + bool save_template_images = debug_template_images; if(save_template_images) { @@ -242,11 +280,7 @@ std::string run_auto_track(tipl::program_option& po,const std::vector return handle->error_msg + " at " + fib_file_name; if(handle->s2t.empty()) return std::string("no subject-to-template mapping at ") + fib_file_name; - std::string output_base = template_images_option; - if(output_base.empty()) - output_base = dir + "/" + fib_base + ".template"; - else if(std::filesystem::is_directory(output_base)) - output_base = (std::filesystem::path(output_base) / (fib_base + ".template")).string(); + std::string output_base = dir + "/" + fib_base + ".template"; tipl::out() << "saving template images in native space to " << output_base; { auto warped_qa = tipl::compose_mapping(handle->template_I,handle->s2t); @@ -272,25 +306,7 @@ std::string run_auto_track(tipl::program_option& po,const std::vector } if(!handle->load_track_atlas(true/*symmetric*/)) return handle->error_msg + " at " + fib_file_name; - std::string atlas_track_name = atlas_track_option; - if(atlas_track_name.empty()) - atlas_track_name = dir + "/" + fib_base + ".atlas.tt.gz"; - else if(std::filesystem::is_directory(atlas_track_name)) - atlas_track_name = (std::filesystem::path(atlas_track_name) / (fib_base + ".atlas.tt.gz")).string(); - if(!tipl::ends_with(atlas_track_name,"tt.gz") && - !tipl::ends_with(atlas_track_name,".trk") && - !tipl::ends_with(atlas_track_name,".trk.gz") && - !tipl::ends_with(atlas_track_name,".tck") && - !tipl::ends_with(atlas_track_name,".txt") && - !tipl::ends_with(atlas_track_name,".mat") && - !tipl::ends_with(atlas_track_name,".nii") && - !tipl::ends_with(atlas_track_name,".nii.gz")) - { - std::string track_ext = trk_format; - if(!track_ext.empty() && track_ext[0] != '.') - track_ext = "." + track_ext; - atlas_track_name += track_ext.empty() ? ".tt.gz" : track_ext; - } + std::string atlas_track_name = dir + "/" + fib_base + ".atlas.tt.gz"; tipl::out() << "saving template tracks in native space to " << atlas_track_name; if(!handle->track_atlas || !handle->track_atlas->save_tracts_to_file(atlas_track_name)) return std::string("failed to save atlas track to ") + atlas_track_name; @@ -503,47 +519,19 @@ std::string run_auto_track(tipl::program_option& po,const std::vector thread.param.step_size = po.get("step_size",thread.param.step_size); thread.param.smooth_fraction = po.get("smoothing",thread.param.smooth_fraction); - if(chen_mode) - { - auto track_ids = handle->get_track_ids(tract_name); - float min_l = 0.0f,max_l = 0.0f; - for(size_t idx = 0;idx < track_ids.size();++idx) - { - auto id = track_ids[idx]; - if(idx == 0) - { - min_l = handle->tract_atlas_min_length[id]; - max_l = handle->tract_atlas_max_length[id]; - } - else - { - min_l = std::min(min_l,handle->tract_atlas_min_length[id]); - max_l = std::max(max_l,handle->tract_atlas_max_length[id]); - } - } - thread.param.min_length = handle->vs[0]*std::max(tolerance[tracking_iteration], - min_l-2.0f*tolerance[tracking_iteration])/handle->tract_atlas_jacobian; - thread.param.max_length = handle->vs[0]*(max_l+2.0f*tolerance[tracking_iteration])/handle->tract_atlas_jacobian; - } - else - { - auto minmax = handle->get_track_minmax_length(tract_name); - thread.param.min_length = handle->vs[0]*std::max(tolerance[tracking_iteration], - minmax.first-2.0f*tolerance[tracking_iteration])/handle->tract_atlas_jacobian; - thread.param.max_length = handle->vs[0]*(minmax.second+2.0f*tolerance[tracking_iteration])/handle->tract_atlas_jacobian; - } + auto minmax = handle->get_track_minmax_length(tract_name); + thread.param.min_length = handle->vs[0]*std::max(tolerance[tracking_iteration], + minmax.first-2.0f*tolerance[tracking_iteration])/handle->tract_atlas_jacobian; + thread.param.max_length = handle->vs[0]*(minmax.second+2.0f*tolerance[tracking_iteration])/handle->tract_atlas_jacobian; tipl::out() << "min_length(mm): " << thread.param.min_length << std::endl; tipl::out() << "max_length(mm): " << thread.param.max_length << std::endl; - thread.param.tip_iteration = po.get("tip_iteration",chen_mode ? 4 : 32); - if(chen_mode) - thread.param.check_ending = po.get("check_ending",1) && !tipl::contains_case_insensitive(tract_name,"Cingulum"); - else - thread.param.check_ending = po.get("check_ending",1); - thread.param.track_voxel_ratio = po.get("track_voxel_ratio",chen_mode ? 2.0f : thread.param.track_voxel_ratio); + thread.param.tip_iteration = po.get("tip_iteration",32); + thread.param.check_ending = po.get("check_ending",1); + thread.param.track_voxel_ratio = po.get("track_voxel_ratio",thread.param.track_voxel_ratio); } { thread.roi_mgr->use_auto_track = true; - thread.roi_mgr->use_roi = po.get("use_roi",chen_mode ? 0 : 1); + thread.roi_mgr->use_roi = po.get("use_roi",1); thread.roi_mgr->tract_name = tract_name; thread.roi_mgr->tolerance_dis_in_icbm152_mm = tolerance[tracking_iteration]; } diff --git a/reg.hpp b/reg.hpp index 747b5440..23f78dd2 100644 --- a/reg.hpp +++ b/reg.hpp @@ -30,46 +30,6 @@ inline auto template_image_pre(const tipl::image& I) return template_image_pre(tipl::image(I)); } -inline size_t linear_with_mi(const tipl::image<3,float>& from, - const tipl::vector<3>& from_vs, - const tipl::image<3,float>& to, - const tipl::vector<3>& to_vs, - tipl::affine_transform& arg, - tipl::reg::reg_type reg_type, - bool& terminated, - const float (*bound)[8] = tipl::reg::reg_bound) -{ - auto new_to_vs = to_vs; - if(reg_type == tipl::reg::affine) - new_to_vs = tipl::reg::adjust_to_vs(from,from_vs,to,to_vs); - if(new_to_vs != to_vs) - arg = tipl::transformation_matrix(arg,from,from_vs,to,to_vs).to_affine_transform(from,from_vs,to,new_to_vs); - float result = tipl::reg::linear_reg(tipl::reg::make_list(from),from_vs, - tipl::reg::make_list(to),new_to_vs, - arg,bound,reg_type,tipl::reg::mutual_info, - true,has_cuda,terminated); - if(new_to_vs != to_vs) - arg = tipl::transformation_matrix(arg,from,from_vs,to,new_to_vs).to_affine_transform(from,from_vs,to,to_vs); - return result; -} - -inline size_t linear_with_mi(const tipl::image<3,float>& from, - const tipl::vector<3>& from_vs, - const tipl::image<3,float>& to, - const tipl::vector<3>& to_vs, - tipl::transformation_matrix& T, - tipl::reg::reg_type reg_type, - bool& terminated, - const float (*bound)[8] = tipl::reg::reg_bound) -{ - tipl::affine_transform arg; - size_t result = linear_with_mi(from,from_vs,to,to_vs,arg,reg_type,terminated,bound); - tipl::out() << arg << std::endl; - T = tipl::transformation_matrix(arg,from.shape(),from_vs,to.shape(),to_vs); - tipl::out() << T << std::endl; - return result; -} - extern int map_ver; struct dual_reg{ static constexpr int dimension = 3; @@ -178,11 +138,4 @@ struct dual_reg{ -namespace tipl -{ -namespace reg -{ -} -} - #endif//REG_HPP