reorient vector field after transformation

vectorreorient.cpp

/* Copyright (c) 2008-2020 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/.
 *
 * Covered Software is provided under this License on an "as is"
 * basis, without warranty of any kind, either expressed, implied, or
 * statutory, including, without limitation, warranties that the
 * Covered Software is free of defects, merchantable, fit for a
 * particular purpose or non-infringing.
 * See the Mozilla Public License v. 2.0 for more details.
 *
 * For more details, see http://www.mrtrix.org/.
 */

#include "command.h"
#include "progressbar.h"
#include "algo/loop.h"
#include "image.h"
#include "adapter/jacobian.h"
#include "registration/warp/helpers.h"

using namespace MR;
using namespace App;


void usage ()
{
  AUTHOR = "Max Pietsch, Dave Raffelt";

  SYNOPSIS = "Reorient a vector field using warp";

  DESCRIPTION
  + "Reorientation is performed by transforming the vector using "
    "the Jacobian (local affine transform) computed at each voxel in the warp.";

  ARGUMENTS
  + Argument ("vector_in", "the input vector file").type_image_in()
  + Argument ("warp", "a 4D deformation field used to perform reorientation.").type_image_in ()
  + Argument ("vector_out", "the output vector.").type_image_out();

  OPTIONS
  + Option ("modulate", "modulate by the Jacobian of the warp.");
}


void run ()
{
  auto input_vec = Image<float>::open(argument[0]);
  if (input_vec.ndim() != 4) {
    throw Exception("Input image needs to be 4D");
  if (input_vec.size(3) != 3) {
    throw Exception("Input vector field needs 3 volumes");
  }

  bool modulate = get_options("modulate").size();

  Header warp_header = Header::open (argument[1]);
  Registration::Warp::check_warp (warp_header);
  check_dimensions (input_vec, warp_header, 0, 3);


  Adapter::Jacobian<Image<float> > jacobian (warp_header.get_image<float>());
  auto output_vec = Image<float>::create(argument[2], input_vec).with_direct_io();

  for (auto l = Loop ("reorienting vector", input_vec, 0, 3) (input_vec, output_vec, jacobian); l; ++l){
    Eigen::Matrix<float, 3, 3> transform = jacobian.value().inverse();
    output_vec.row(3) = (transform * Eigen::Vector3f (input_vec.row(3)));
    if (!modulate)
      output_vec.row(3) /= transform.determinant();
  }
}