1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
|
/* Copyright (c) 2008-2025 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 "header.h"
#include "image.h"
#include "image_helpers.h"
#include "transform.h"
#include "types.h"
#include "connectome/connectome.h"
using namespace MR;
using namespace App;
const char * field_choices[] = { "mass", "centre", nullptr };
void usage ()
{
AUTHOR = "Robert E. Smith (robert.smith@florey.edu.au)";
SYNOPSIS = "Compute statistics of parcels within a label image";
ARGUMENTS
+ Argument ("input", "the input label image").type_image_in();
OPTIONS
+ Option ("output", "output only the field specified; "
"options are: " + join(field_choices, ","))
+ Argument ("choice").type_choice (field_choices)
+ Option ("voxelspace", "report parcel centres of mass in voxel space rather than scanner space");
}
using Connectome::node_t;
using vector_type = Eigen::Array<default_type, Eigen::Dynamic, 1>;
using matrix_type = Eigen::Matrix<default_type, Eigen::Dynamic, 3>;
void run ()
{
Header H = Header::open (argument[0]);
if (H.ndim() > 3)
throw Exception ("Command does not accept images with more than 3 dimensions");
Connectome::check (H);
Image<node_t> image = H.get_image<node_t>();
matrix_type coms;
vector_type masses;
for (auto l = Loop(image) (image); l; ++l) {
const node_t value = image.value();
if (value) {
if (value > coms.rows()) {
coms.conservativeResizeLike (matrix_type::Zero (value, 3));
masses.conservativeResizeLike (vector_type::Zero (value));
}
coms.row(value-1) += Eigen::Vector3d (image.index(0), image.index(1), image.index(2));
masses[value-1]++;
}
}
coms = coms.array().colwise() / masses;
if (!get_options("voxelspace").size())
coms = (image.transform() * coms.transpose()).transpose();
auto opt = get_options ("output");
if (opt.size()) {
switch (int(opt[0][0])) {
case 0: std::cout << masses; break;
case 1: std::cout << coms; break;
default: assert (0);
}
return;
}
// Ensure that the printout of centres-of-mass is nicely formatted
Eigen::IOFormat fmt (Eigen::StreamPrecision, 0, ", ", "\n", "[ ", " ]", "", "");
std::stringstream com_stringstream;
com_stringstream << coms.format (fmt);
const vector<std::string> com_strings = split(com_stringstream.str(), "\n");
assert (com_strings.size() == size_t(masses.size()));
// Find width of first non-empty string, in order to centralise header label
size_t com_width = 0;
for (const auto& i : com_strings) {
if (i.size()) {
com_width = i.size();
break;
}
}
std::cout << std::setw(8) << std::right << "index" << " "
<< std::setw(8) << std::right << "mass" << " "
<< std::setw(4) << std::right << " "
<< std::setw(com_width/2 + std::string("centre of mass").size()/2) << std::right << "centre of mass" << "\n";
for (ssize_t i = 0; i != masses.size(); ++i) {
if (masses[i]) {
std::cout << std::setw(8) << std::right << i+1 << " "
<< std::setw(8) << std::right << masses[i] << " "
<< std::setw(4+com_width) << std::right << com_strings[i] << "\n";
}
}
}
|
