<commit_before>/*
 * Copyright (c) 2008-2018 the MRtrix3 contributors.
 *
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, you can obtain one at http://mozilla.org/MPL/2.0/
 *
 * MRtrix3 is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty
 * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 *
 * For more details, see http://www.mrtrix.org/
 */


#ifndef __dwi_tractography_sift_output_h__
#define __dwi_tractography_sift_output_h__


#include "image.h"
#include "header.h"

#include "algo/loop.h"

#include "dwi/fixel_map.h"
#include "dwi/tractography/SIFT/model_base.h"

#include "file/ofstream.h"

#include "math/SH.h"

#include "fixel/legacy/fixel_metric.h"
#include "fixel/legacy/image.h"
#include "fixel/legacy/keys.h"


namespace MR
{
  namespace DWI
  {
    namespace Tractography
    {
      namespace SIFT
      {

      // Output functions - non-essential, mostly debugging outputs
      template <class Fixel>
      void ModelBase<Fixel>::output_target_image (const std::string& path) const
      {
        auto out = Image<float>::create (path, Fixel_map<Fixel>::header());
        VoxelAccessor v (accessor());
        for (auto l = Loop(out) (out, v); l; ++l) {
          if (v.value()) {
            default_type value = 0.0;
            for (typename Fixel_map<Fixel>::ConstIterator i = begin (v); i; ++i)
              value += i().get_FOD();
            out.value() = value;
          } else {
            out.value() = NAN;
          }
        }
      }

      template <class Fixel>
      void ModelBase<Fixel>::output_target_image_sh (const std::string& path) const
      {
        const size_t L = 8;
        const size_t N = Math::SH::NforL (L);
        Math::SH::aPSF<default_type> aPSF (L);
        Header H_sh (Fixel_map<Fixel>::header());
        H_sh.ndim() = 4;
        H_sh.size(3) = N;
        H_sh.stride (3) = 0;
        auto out = Image<float>::create (path, H_sh);
        VoxelAccessor v (accessor());
        for (auto l = Loop (0, 3) (out, v); l; ++l) {
          if (v.value()) {
            Eigen::Matrix<default_type, Eigen::Dynamic, 1> sum = Eigen::Matrix<default_type, Eigen::Dynamic, 1>::Zero (N);
            for (typename Fixel_map<Fixel>::ConstIterator i = begin (v); i; ++i) {
              if (i().get_FOD()) {
                Eigen::Matrix<default_type, Eigen::Dynamic, 1> this_lobe;
                aPSF (this_lobe, i().get_dir());
                for (size_t c = 0; c != N; ++c)
                  sum[c] += i().get_FOD() * this_lobe[c];
              }
            }
            for (auto l = Loop (3) (out); l; ++l)
              out.value() = sum[out.index(3)];
          } else {
            for (auto l = Loop (3) (out); l; ++l)
              out.value() = NAN;
          }
        }
      }

      template <class Fixel>
      void ModelBase<Fixel>::output_target_image_fixel (const std::string& path) const
      {
        using MR::Fixel::Legacy::FixelMetric;
        Header H_fixel (Fixel_map<Fixel>::header());
        H_fixel.datatype() = DataType::UInt64;
        H_fixel.datatype().set_byte_order_native();
        H_fixel.keyval()[MR::Fixel::Legacy::name_key] = str(typeid(FixelMetric).name());
        H_fixel.keyval()[MR::Fixel::Legacy::size_key] = str(sizeof(FixelMetric));
        MR::Fixel::Legacy::Image<FixelMetric> out (path, H_fixel);
        VoxelAccessor v (accessor());
        for (auto l = Loop (out) (out, v); l; ++l) {
          if (v.value()) {
            out.value().set_size ((*v.value()).num_fixels());
            size_t index = 0;
            for (typename Fixel_map<Fixel>::ConstIterator iter = begin (v); iter; ++iter, ++index) {
              FixelMetric fixel (iter().get_dir().template cast<float>(), iter().get_FOD(), iter().get_FOD());
              out.value()[index] = fixel;
            }
          }
        }
      }

      template <class Fixel>
      void ModelBase<Fixel>::output_tdi (const std::string& path) const
      {
        const default_type current_mu = mu();
        auto out = Image<float>::create (path, Fixel_map<Fixel>::header());
        VoxelAccessor v (accessor());
        for (auto l = Loop (out) (out, v); l; ++l) {
          if (v.value()) {
            default_type value = 0.0;
            for (typename Fixel_map<Fixel>::ConstIterator i = begin (v); i; ++i)
              value += i().get_TD();
            out.value() = value * current_mu;
          } else {
            out.value() = NaN;
          }
        }
      }

      template <class Fixel>
      void ModelBase<Fixel>::output_tdi_null_lobes (const std::string& path) const
      {
        const default_type current_mu = mu();
        auto out = Image<float>::create (path, Fixel_map<Fixel>::header());
        VoxelAccessor v (accessor());
        for (auto l = Loop (out) (out, v); l; ++l) {
          if (v.value()) {
            default_type value = 0.0;
            for (typename Fixel_map<Fixel>::ConstIterator i = begin (v); i; ++i) {
              if (!i().get_FOD())
                value += i().get_TD();
            }
            out.value() = value * current_mu;
          } else {
            out.value() = NaN;
          }
        }
      }

      template <class Fixel>
      void ModelBase<Fixel>::output_tdi_sh (const std::string& path) const
      {
        const default_type current_mu = mu();
        const size_t L = 8;
        const size_t N = Math::SH::NforL (L);
        Math::SH::aPSF<default_type> aPSF (L);
        Header H_sh (Fixel_map<Fixel>::header());
        H_sh.ndim() = 4;
        H_sh.size(3) = N;
        H_sh.stride (3) = 0;
        auto out = Image<float>::create (path, H_sh);
        VoxelAccessor v (accessor());
        for (auto l = Loop (out) (out, v); l; ++l) {
          if (v.value()) {
            Eigen::Matrix<default_type, Eigen::Dynamic, 1> sum = Eigen::Matrix<default_type, Eigen::Dynamic, 1>::Zero (N);
            for (typename Fixel_map<Fixel>::ConstIterator i = begin (v); i; ++i) {
              if (i().get_FOD()) {
                Eigen::Matrix<default_type, Eigen::Dynamic, 1> this_lobe;
                aPSF (this_lobe, i().get_dir());
                for (size_t c = 0; c != N; ++c)
                  sum[c] += i().get_TD() * this_lobe[c];
              }
              for (auto l = Loop (3) (out); l; ++l)
                out.value() = sum[out.index(3)] * current_mu;
            }
          } else {
            for (auto l = Loop (3) (out); l; ++l)
              out.value() = NaN;
          }
        }
      }

      template <class Fixel><commit_msg>SIFT: Fix TDI SH image output<commit_after>/*
 * Copyright (c) 2008-2018 the MRtrix3 contributors.
 *
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, you can obtain one at http://mozilla.org/MPL/2.0/
 *
 * MRtrix3 is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty
 * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 *
 * For more details, see http://www.mrtrix.org/
 */


#ifndef __dwi_tractography_sift_output_h__
#define __dwi_tractography_sift_output_h__


#include "image.h"
#include "header.h"

#include "algo/loop.h"

#include "dwi/fixel_map.h"
#include "dwi/tractography/SIFT/model_base.h"

#include "file/ofstream.h"

#include "math/SH.h"

#include "fixel/legacy/fixel_metric.h"
#include "fixel/legacy/image.h"
#include "fixel/legacy/keys.h"


namespace MR
{
  namespace DWI
  {
    namespace Tractography
    {
      namespace SIFT
      {

      // Output functions - non-essential, mostly debugging outputs
      template <class Fixel>
      void ModelBase<Fixel>::output_target_image (const std::string& path) const
      {
        auto out = Image<float>::create (path, Fixel_map<Fixel>::header());
        VoxelAccessor v (accessor());
        for (auto l = Loop(out) (out, v); l; ++l) {
          if (v.value()) {
            default_type value = 0.0;
            for (typename Fixel_map<Fixel>::ConstIterator i = begin (v); i; ++i)
              value += i().get_FOD();
            out.value() = value;
          } else {
            out.value() = NAN;
          }
        }
      }

      template <class Fixel>
      void ModelBase<Fixel>::output_target_image_sh (const std::string& path) const
      {
        const size_t L = 8;
        const size_t N = Math::SH::NforL (L);
        Math::SH::aPSF<default_type> aPSF (L);
        Header H_sh (Fixel_map<Fixel>::header());
        H_sh.ndim() = 4;
        H_sh.size(3) = N;
        H_sh.stride (3) = 0;
        auto out = Image<float>::create (path, H_sh);
        VoxelAccessor v (accessor());
        for (auto l = Loop (0, 3) (out, v); l; ++l) {
          if (v.value()) {
            Eigen::Matrix<default_type, Eigen::Dynamic, 1> sum = Eigen::Matrix<default_type, Eigen::Dynamic, 1>::Zero (N);
            for (typename Fixel_map<Fixel>::ConstIterator i = begin (v); i; ++i) {
              if (i().get_FOD()) {
                Eigen::Matrix<default_type, Eigen::Dynamic, 1> this_lobe;
                aPSF (this_lobe, i().get_dir());
                for (size_t c = 0; c != N; ++c)
                  sum[c] += i().get_FOD() * this_lobe[c];
              }
            }
            for (auto l = Loop (3) (out); l; ++l)
              out.value() = sum[out.index(3)];
          } else {
            for (auto l = Loop (3) (out); l; ++l)
              out.value() = NAN;
          }
        }
      }

      template <class Fixel>
      void ModelBase<Fixel>::output_target_image_fixel (const std::string& path) const
      {
        using MR::Fixel::Legacy::FixelMetric;
        Header H_fixel (Fixel_map<Fixel>::header());
        H_fixel.datatype() = DataType::UInt64;
        H_fixel.datatype().set_byte_order_native();
        H_fixel.keyval()[MR::Fixel::Legacy::name_key] = str(typeid(FixelMetric).name());
        H_fixel.keyval()[MR::Fixel::Legacy::size_key] = str(sizeof(FixelMetric));
        MR::Fixel::Legacy::Image<FixelMetric> out (path, H_fixel);
        VoxelAccessor v (accessor());
        for (auto l = Loop (out) (out, v); l; ++l) {
          if (v.value()) {
            out.value().set_size ((*v.value()).num_fixels());
            size_t index = 0;
            for (typename Fixel_map<Fixel>::ConstIterator iter = begin (v); iter; ++iter, ++index) {
              FixelMetric fixel (iter().get_dir().template cast<float>(), iter().get_FOD(), iter().get_FOD());
              out.value()[index] = fixel;
            }
          }
        }
      }

      template <class Fixel>
      void ModelBase<Fixel>::output_tdi (const std::string& path) const
      {
        const default_type current_mu = mu();
        auto out = Image<float>::create (path, Fixel_map<Fixel>::header());
        VoxelAccessor v (accessor());
        for (auto l = Loop (out) (out, v); l; ++l) {
          if (v.value()) {
            default_type value = 0.0;
            for (typename Fixel_map<Fixel>::ConstIterator i = begin (v); i; ++i)
              value += i().get_TD();
            out.value() = value * current_mu;
          } else {
            out.value() = NaN;
          }
        }
      }

      template <class Fixel>
      void ModelBase<Fixel>::output_tdi_null_lobes (const std::string& path) const
      {
        const default_type current_mu = mu();
        auto out = Image<float>::create (path, Fixel_map<Fixel>::header());
        VoxelAccessor v (accessor());
        for (auto l = Loop (out) (out, v); l; ++l) {
          if (v.value()) {
            default_type value = 0.0;
            for (typename Fixel_map<Fixel>::ConstIterator i = begin (v); i; ++i) {
              if (!i().get_FOD())
                value += i().get_TD();
            }
            out.value() = value * current_mu;
          } else {
            out.value() = NaN;
          }
        }
      }

      template <class Fixel>
      void ModelBase<Fixel>::output_tdi_sh (const std::string& path) const
      {
        const default_type current_mu = mu();
        const size_t L = 8;
        const size_t N = Math::SH::NforL (L);
        Math::SH::aPSF<default_type> aPSF (L);
        Header H_sh (Fixel_map<Fixel>::header());
        H_sh.ndim() = 4;
        H_sh.size(3) = N;
        H_sh.stride (3) = 0;
        auto out = Image<float>::create (path, H_sh);
        VoxelAccessor v (accessor());
        for (auto l = Loop (v) (out, v); l; ++l) {
          if (v.value()) {
            Eigen::Matrix<default_type, Eigen::Dynamic, 1> sum = Eigen::Matrix<default_type, Eigen::Dynamic, 1>::Zero (N);
            for (typename Fixel_map<Fixel>::ConstIterator i = begin (v); i; ++i) {
              if (i().get_FOD()) {
                Eigen::Matrix<default_type, Eigen::Dynamic, 1> this_lobe;
                aPSF (this_lobe, i().get_dir());
                for (size_t c = 0; c != N; ++c)
                  sum[c] += i().get_TD() * this_lobe[c];
              }
              for (auto l = Loop (3) (out); l; ++l)
                out.value() = sum[out.index(3)] * current_mu;
            }
          } else {
            for (auto l = Loop (3) (out); l; ++l)
              out.value() = NaN;
          }
        }
      }

      template <class Fixel>
{'fraction_non_alphanumeric': 0.1532426778242678, 'fraction_numerical': 0.007322175732217573, 'mean_word_length': 2.0392824704813806, 'pattern_counts': {'":': 0, '<': 83, '<?xml version=': 0, '>': 83, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2118277', 'n_tokens_mistral': 4550, 'n_tokens_neox': 4399, 'n_words': 1090}