Add debugging options for AutoTrack and add random starts to linear reg

cmd/reg.cpp

@@ -1,5 +1,7 @@
 #include <QString>
 #include <QImage>
+#include <limits>
+#include <random>
 #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::out>(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<float>::max();
+        auto run_linear = [&](tipl::affine_transform<float,dimension>& start_arg)
+        {
+            auto local_arg = start_arg;
+            float local_cost = tipl::reg::linear<tipl::out>(
+                        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<float> 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<float,dimension> 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();
 

libs/dsi/gqi_mni_reconstruction.hpp

@@ -2,6 +2,7 @@
 #define MNI_RECONSTRUCTION_HPP
 #include <QFileInfo>
 #include <chrono>
+#include <cmath>
 #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;

libs/tracking/fib_data.cpp

@@ -1,5 +1,7 @@
 #include <filesystem>
 #include <unordered_set>
+#include <limits>
+#include <cmath>
 #include <QCoreApplication>
 #include <QFileInfo>
 #include <QDateTime>
@@ -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<double>::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<double>::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<float> min_length(tractography_name_list.size()),max_length(tractography_name_list.size());