Skip to content

Instantly share code, notes, and snippets.

What would you like to do?
MRtrix3 stand-alone command for calculating the centre of mass of an image
/* Copyright (c) 2008-2017 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
* MRtrix is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty
* For more details, see
#include "command.h"
#include "image.h"
#include "image_helpers.h"
#include "transform.h"
#include "types.h"
using namespace MR;
using namespace App;
const char* space_options[] = { "scanner", "voxel", nullptr };
void usage ()
AUTHOR = "Robert E. Smith (";
SYNOPSIS = "Determine the centre of mass / centre of gravity of an image";
+ Argument ("input", "the input image").type_image_in()
+ Argument ("space", "the coordinate space in which to provide the centre (options are: " + join (space_options, ",") + ")").type_choice (space_options);
+ Option ("mask", "only include voxels within a mask in the calculation")
+ Argument ("image").type_image_in();
typedef float value_type;
void run ()
Image<value_type> image = Image<value_type>::open (argument[0]);
if (image.ndim() > 3)
throw Exception ("Command does not accept images with more than 3 dimensions");
Image<bool> mask;
auto opt = get_options ("mask");
if (opt.size()) {
mask = Image<bool>::open (opt[0][0]);
check_dimensions (image, mask);
Eigen::Vector3 com (0.0, 0.0, 0.0);
default_type mass = 0.0;
if (mask.valid()) {
for (auto l = Loop(image) (image, mask); l; ++l) {
if (mask.value()) {
com += Eigen::Vector3 (image.index(0), image.index(1), image.index(2)) * image.value();
mass += image.value();
} else {
for (auto l = Loop(image) (image); l; ++l) {
com += Eigen::Vector3 (image.index(0), image.index(1), image.index(2)) * image.value();
mass += image.value();
com /= mass;
if (!int(argument[1]))
com = image.transform() * com;
std::cout << com.transpose();
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.