diff --git a/Applications/CMakeLists.txt b/Applications/CMakeLists.txt new file mode 100644 index 000000000..5ce41433b --- /dev/null +++ b/Applications/CMakeLists.txt @@ -0,0 +1,148 @@ +# ============================================================================ +# Medical Image Registration ToolKit (MIRTK) +# +# Copyright 2013-2015 Imperial College London +# Copyright 2013-2015 Andreas Schuh +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ============================================================================ + +############################################################################## +# @file CMakeLists.txt +# @brief Build configuration of (core) MIRTK Applications. +############################################################################## + +# ---------------------------------------------------------------------------- +# Include optional dependencies (cf. Config/Depends.cmake) +if (MODULE_PointSet) + basis_use_package(VTK) + if (FLANN_FOUND) + basis_use_package(FLANN) + else () + set(FLANN_LIBRARY) + endif () + if (MATLAB_FOUND) + basis_use_package(MATLAB) + else () + set(MATLAB_mwmclmcrrt_LIBRARY) + endif () +endif () + +## NOTE: OPTIONAL dependencies on MIRTK IO modules added by mirtk_add_executable +## automatically when list of dependencies includes LibImage. + +# Generic file info tool +mirtk_add_executable(info DEPENDS LibImage OPTIONAL LibTransformation LibPointSet ${VTK_LIBRARIES}) + +# ---------------------------------------------------------------------------- +# Image tools +if (MODULE_Image) + + macro(add_image_command cmd) + mirtk_add_executable(${cmd} DEPENDS LibImage ${ARGN}) + endmacro() + + add_image_command(average-images OPTIONAL LibTransformation) + add_image_command(calculate OPTIONAL ${VTK_LIBRARIES}) + add_image_command(calculate-distance-map) + add_image_command(calculate-exponential-map) + add_image_command(calculate-lie-bracket) + add_image_command(calculate-logarithmic-map) + add_image_command(close-image) + add_image_command(convert-image) + add_image_command(cut-brain) + add_image_command(detect-edges) + add_image_command(dilate-image) + add_image_command(downsample-image) + add_image_command(edit-image OPTIONAL LibTransformation) + add_image_command(erode-image) + add_image_command(evaluate-overlap) + add_image_command(extract-connected-components) + add_image_command(extract-image-region) + add_image_command(flip-image) + add_image_command(concatenate-images) + add_image_command(open-image) + add_image_command(reflect-image) + add_image_command(resample-image) + add_image_command(smooth-image) + add_image_command(transform-image LibTransformation) + +endif () + +# ---------------------------------------------------------------------------- +# Registration tools +if (MODULE_Registration) + + macro(add_registration_command cmd) + mirtk_add_executable(${cmd} DEPENDS LibRegistration ${ARGN}) + endmacro() + + add_registration_command(register LibImage OPTIONAL LibLBFGS LibPointSet LibDeformable) + +endif () + +# ---------------------------------------------------------------------------- +# Transformation tools +if (MODULE_Transformation) + + macro(add_transformation_command cmd) + mirtk_add_executable(${cmd} DEPENDS LibTransformation ${ARGN}) + endmacro() + + add_transformation_command(average-dofs LibImage) + add_transformation_command(concatenate-dofs LibImage) + add_transformation_command(convert-dof LibImage) + add_transformation_command(convert-dof2csv LibImage) + add_transformation_command(convert-dof2velo LibImage) + add_transformation_command(edit-dof) + add_transformation_command(evaluate-dof LibImage) + add_transformation_command(evaluate-jacobian LibImage) + add_transformation_command(init-dof LibImage OPTIONAL LibPointSet) + add_transformation_command(invert-dof) + +endif () + +# ---------------------------------------------------------------------------- +# Point set tools +if (MODULE_PointSet) + + macro(add_pointset_command cmd) + mirtk_add_executable(${cmd} DEPENDS LibPointSet ${VTK_LIBRARIES} ${ARGN}) + endmacro() + + add_pointset_command(calculate-surface-attributes) + add_pointset_command(calculate-surface-spectrum) + add_pointset_command(convert-mris) # requires Python VTK modules and FreeSurfer installation + add_pointset_command(convert-pointset LibImage) + add_pointset_command(copy-pointset-attributes) + add_pointset_command(decimate-surface) + add_pointset_command(delete-pointset-attributes) + add_pointset_command(evaluate-distance) + add_pointset_command(evaluate-distortion) + add_pointset_command(evaluate-surface-overlap ${FLANN_LIBRARY}) + add_pointset_command(extract-connected-points) + add_pointset_command(extract-pointset-surface LibImage) + add_pointset_command(extract-surface LibImage) + add_pointset_command(match-points LibTransformation) + add_pointset_command(offset-surface) + add_pointset_command("project-onto-surface") + add_pointset_command(transform-points LibImage LibTransformation) + add_pointset_command(remesh) + add_pointset_command(smooth-surface) + + if (MATLAB_FOUND) + add_pointset_command(convert-pointset-to-mat ${MATLAB_mwmclmcrrt_LIBRARY}) + add_pointset_command(copy-pointset-attributes-from-mat ${MATLAB_mwmclmcrrt_LIBRARY}) + endif () + +endif () diff --git a/Applications/average-dofs.cc b/Applications/average-dofs.cc new file mode 100644 index 000000000..38ee770fd --- /dev/null +++ b/Applications/average-dofs.cc @@ -0,0 +1,753 @@ +/* + * Medical Image Registration ToolKit (MIRTK) + * + * Copyright 2008-2015 Imperial College London + * Copyright 2008-2013 Daniel Rueckert, Julia Schnabel + * Copyright 2013-2015 Andreas Schuh + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +#include +#include +#include +#include +#include + +using namespace mirtk; + + +// =========================================================================== +// Constants +// =========================================================================== + +const double NaN = numeric_limits::quiet_NaN(); +const double EPSILON = 0.001; // default -epsilon value + +// =========================================================================== +// Help +// =========================================================================== + +// --------------------------------------------------------------------------- +void PrintHelp(const char *name) +{ + cout << endl; + cout << "Usage: " << name << " ... [options]" << endl; + cout << " " << name << " -dofnames [options]" << endl; + cout << endl; + cout << "Description:" << endl; + cout << " This command averages a number of input transformations." << endl; + cout << " It can be used to construct a brain image atlas and has been" << endl; + cout << " utilized for the construction of the spatio-temporal fetal/neonatal" << endl; + cout << " brain atlases available at http://brain-development.org/brain-atlases/." << endl; + cout << endl; + cout << " The transformations which can be averaged by this program are listed below," << endl; + cout << " commonly expressed as \"sum of global and local transformations\", i.e.," << endl; + cout << endl; + cout << " .. math::" << endl; + cout << endl; + cout << " T(x) = T_{global}(x) + T_{local}(x)" << endl; + cout << endl; + cout << " ====== ===============================================================================" << endl; + cout << " Linear :math:`T(x) = Ax, T_{global}(x) = A x, T_{local}(x) = 0`" << endl; + cout << " FFD :math:`T(x) = x + d(x), T_{global}(x) = I x, T_{local}(x) = d(x)`" << endl; + cout << " MFFD :math:`T(x) = Ax + d(x), T_{global}(x) = A x, T_{local}(x) = d(x)`" << endl; + cout << " Fluid :math:`T(x) = x + d(Ax), T_{global}(x) = A x, T_{local}(x) = d(Ax) - (Ax - x)`" << endl; + cout << " ====== ===============================================================================" << endl; + cout << endl; + cout << " where A is a linear 4x4 transformation matrix, and d(x) is a non-linear displacement." << endl; + cout << " Note that in case of the fluid transformation, d(x) is the total displacement" << endl; + cout << " instead of the local displacement only. By subtracting the displacement induced" << endl; + cout << " by the global transformation, the displacement corresponding to the local" << endl; + cout << " deformation only is obtained. Note that in case of MFFDs, the total displacement" << endl; + cout << " introduced by all local transformation levels is used to compute the average local" << endl; + cout << " transformation and thus multiple levels are allowed. Representing all transformations" << endl; + cout << " as sum of global and local transformations enables the computation of an average" << endl; + cout << " transformation for a given target domain from an arbitrary mix of input transformations." << endl; + // + // Note: This computation is in fact implemented by the overridden virtual member + // function FluidFreeFormTransformation::LocalDisplacement and hence + // the virtual function Transformation::LocalDisplacement can be used + // in all cases to get the local component of the total displacement field. + // + cout << endl; + cout << " Accordingly, the resulting average transformation is set to be equal the sum of the" << endl; + cout << " Log-Euclidean means of the global transformation matrices, A, of the input transformations" << endl; + cout << " and the corresponding local transformations. The bi-invariant global mean is computed" << endl; + cout << " instead using the barycentric fixed point iteration when :option:`-bi-invariant` is given." << endl; + cout << " Before the computation of the local mean transformation, any dependency of the local" << endl; + cout << " displacement on the global component is removed by pre-multiplying the displacement" << endl; + cout << " vectors by the inverse matrix of the global transformation [1]_. Moreover, the (optional)" << endl; + cout << " local component of the final average transformation has to be taken into consideration" << endl; + cout << " as well. Therefore, all local input displacements are made relative to the same position" << endl; + cout << " after global transformation, i.e.," << endl; + cout << endl; + cout << " .. math::" << endl; + cout << endl; + cout << " d(x) = \\bar{A} \\circ A^{-1} \\circ T_{local}(x)" << endl; + cout << endl; + cout << " The average of all transformation parameters is computed in the \"log-space\" by default." << endl; + cout << " Thus, in case of the local displacement fields, d(x), the corresponding stationary" << endl; + cout << " velocity fields, :math:`v(x) = log(d(x))`, are averaged instead. If the input transformation" << endl; + cout << " itself is parameterized by stationary velocities, these velocities are used directly" << endl; + cout << " when possible to avoid the redundant :math:`log(exp)` computation and its numerical error." << endl; + cout << " The resulting average velocity field is then either combined with a suitable" << endl; + cout << " velocity-based transformation model or exponentiated again to obtain the final" << endl; + cout << " average displacement field, :math:`\\bar{d}(x)`." << endl; + cout << endl; + cout << " The output transformation will be an affine transformation if only the global transformation" << endl; + cout << " parameters were averaged, or a MFFD with one level corresponding to the average FFD otherwise." << endl; + cout << " The type of the output FFD depends on the type of the input transformations. If all transformations" << endl; + cout << " are either a linear transformation or a FFD with cubic B-spline or linear interpolation kernel" << endl; + cout << " (may also be contained in a one-level MFFD), the ouput FFD will contain the average parameters" << endl; + cout << " of the input FFDs. If the average is computed in the log-space, this only applies if all input" << endl; + cout << " FFDs are parameterized by stationary velocity fields. In all other cases, the output FFD will" << endl; + cout << " be either defined in the image domain of the input target image or the first encountered input" << endl; + cout << " FFD with a linear interpolation kernel and a control point spacing corresponding to the voxel size." << endl; + cout << " In other words, the output FFD is in general a dense (\"non-parametric\") displacement field." << endl; + cout << endl; + cout << " .. [1] Rueckert et al., Automatic construction of 3-D statistical deformation" << endl; + cout << " models of the brain using nonrigid registration, IEEE TMI, 22(8), 1014–25 (2003)" << endl; + cout << endl; + cout << "Input options:" << endl; + cout << " -target Common reference image space of input transformations." << endl; + cout << " If not specified and local transformations are to be averaged," << endl; + cout << " the attributes of the first local transformation are used." << endl; + cout << " -dofnames Text file listing input transformations and associated values." << endl; + cout << " -dofdir Directory used to make relative paths in :option:`-dofnames` text file absolute. (default: cwd)" << endl; + cout << " -prefix Prefix for transformation name entries in :option:`-dofnames` text file. (default: none)" << endl; + cout << " -suffix Suffix for transformation name entries in :option:`-dofnames` text file. (default: none)" << endl; + cout << " -gaussian Use Gaussian kernel weights. Requires :option:`-dofnames` to specify transformation values." << endl; + cout << " By default, if :option:`-dofnames` is used or is 0, the values specified in the" << endl; + cout << " text file are directly used as kernel weights for the averaging instead of" << endl; + cout << " using these values as arguments for the Gaussian kernel function." << endl; + cout << " -epsilon Weight threshold. (default: " << EPSILON << ")" << endl; + cout << " -add-identity Assume additional identity transformation as part of input transformations. (default: none)" << endl; + cout << " -add-identity-with-weight Assume additional identity transformation with given value/weight. (default: none)" << endl; + cout << " -add-identity-for-dofname Assume identity transformation for named input transformation." << endl; + cout << " Note that if this option is used, the named input transformation file" << endl; + cout << " which is listed in the :option:`-dofnames` list does not need to exist." << endl; + cout << " (default: read input transformation from file)" << endl; + cout << endl; + cout << "Average transformation options:" << endl; + cout << " -[no]rotation Average rotation or assume none to be present. (default: off)" << endl; + cout << " -[no]translation Average translation or assume none to be present. (default: off)" << endl; + cout << " -[no]scaling Average scaling or assume none to be present. (default: on)" << endl; + cout << " -[no]shearing Average shearing or assume none to be present. (default: on)" << endl; + cout << " -[no]deformation Average deformation or assume none to be present. (default: on)" << endl; + cout << " -[no]log Whether to average local transformations in log-space." << endl; + cout << " (default: yes, unless the :option:`-dofs` are averaged directly)" << endl; + cout << " -log-euclidean Compute Log-Euclidean means. (default: off)" << endl; + cout << " -bi-invariant Compute global bi-invariant mean, i.e., exponential barycenter. (default: on)" << endl; + cout << " -dofs Average the local transformation parameters directly." << endl; + cout << " (default: corresponding (dense) displacement fields are averaged)" << endl; + cout << " -inverse-dofs Average inverse input transformations and invert the average again. (default: off)" << endl; + cout << endl; + cout << " -rigid Enables averaging of :option:`-translation` and :option:`-rotation` components," << endl; + cout << " and disables averaging of :option:`-scaling`, :option:`-shearing`, and :option:`-deformation`." << endl; + cout << " -norigid Disables averaging of :option:`-translation` and :option:`-rotation` components." << endl; + cout << " -affine Enables averaging of :option:`-translation`, :option:`-rotation`," << endl; + cout << " :option:`-scaling`, and :option:`-shearing`, components and disables output of" << endl; + cout << " average :option:`-deformation`." << endl; + cout << " -noaffine Disables averaging of :option:`-rotation`, :option:`-translation`," << endl; + cout << " :option:`-scaling`, and :option:`-shearing` components." << endl; + cout << " -all Enables averaging of :option:`-rotation`, :option:`-translation`," << endl; + cout << " :option:`-scaling`, :option:`-shearing`, and :option:`-deformation`." << endl; + cout << endl; + cout << " -linear Force average output transformation to use linear interpolation." << endl; + cout << " -cubic Force average output transformation to use cubic B-spline interpolation." << endl; + PrintStandardOptions(cout); + cout << endl; +} + +// =========================================================================== +// Auxiliaries +// =========================================================================== + +// --------------------------------------------------------------------------- +bool isRelativePath(const string &path) +{ + return path.size() > 0 && path[0] != '/'; +} + +// =========================================================================== +// Main +// =========================================================================== + +// --------------------------------------------------------------------------- +int main(int argc, char **argv) +{ + REQUIRES_POSARGS(1); + + string dofout; + Array dofin; + Array value; + + // Parse arguments + dofout = POSARG(1); + for (OPTIONAL_POSARGS) dofin.push_back(ARGUMENT); + + // Parse options + const char *target_name = NULL; + const char *doflist = NULL; + const char *dofdir = NULL; + const char *prefix = NULL; + const char *suffix = NULL; + bool translation = false; + bool rotation = false; + bool scaling = true; + bool shearing = true; + bool deformation = true; + bool invert = false; // invert input transformations + bool invavg = false; // invert output transformation + int logspace = -1; // i.e., not set explicitly to "true" + bool avgdofs = false; // average DoFs directly if possible + bool approxglobal = true; // approximate vs. truncate average + // global transformation if #DoFs < 12 + bool bsplineffd = false; // output B-spline FFD when avgdofs == false + bool biinvariant = true; // bi-invariant mean if possible + const char *identity_name = "identity"; // name of (input) transformation which is to be replaced by identity + double identity_value = numeric_limits::infinity(); // i.e., no implicit identity transformation added by default + double mean = .0; + double sigma = .0; + double epsilon = EPSILON; + int frechet_iter = 20; + + for (ALL_OPTIONS) { + if (OPTION("-target")) target_name = ARGUMENT; + else if (OPTION("-dofnames")) doflist = ARGUMENT; + else if (OPTION("-dofdir")) dofdir = ARGUMENT; + else if (OPTION("-prefix")) prefix = ARGUMENT; + else if (OPTION("-suffix")) suffix = ARGUMENT; + else if (OPTION("-invert")) invert = true; + else if (OPTION("-inverse")) invavg = true; + else if (OPTION("-translation")) translation = true; + else if (OPTION("-rotation")) rotation = true; + else if (OPTION("-scaling")) scaling = true; + else if (OPTION("-shearing")) shearing = true; + else if (OPTION("-deformation")) deformation = true; + else if (OPTION("-notranslation")) translation = false; + else if (OPTION("-norotation")) rotation = false; + else if (OPTION("-noscaling")) scaling = false; + else if (OPTION("-noshearing")) shearing = false; + else if (OPTION("-nodeformation")) deformation = false; + else if (OPTION("-bi-invariant")) biinvariant = true; + else if (OPTION("-log-euclidean")) biinvariant = false, logspace = true; + else if (OPTION("-log")) logspace = true; + else if (OPTION("-nolog")) logspace = false; + else if (OPTION("-linear")) bsplineffd = false; + else if (OPTION("-cubic")) bsplineffd = true; + else if (OPTION("-dofs")) avgdofs = true; + else if (OPTION("-inverse-dofs")) avgdofs = true, invert = invavg = true; + else if (OPTION("-epsilon")) epsilon = atof(ARGUMENT); + else if (OPTION("-gaussian")) { mean = atof(ARGUMENT); sigma = atof(ARGUMENT); } + else if (OPTION("-rigid")) { translation = rotation = true; scaling = shearing = deformation = false; } + else if (OPTION("-norigid")) { translation = rotation = false; } + else if (OPTION("-affine")) { translation = rotation = scaling = shearing = true; deformation = false; } + else if (OPTION("-noaffine")) { translation = rotation = scaling = shearing = false; } + else if (OPTION("-all")) { translation = rotation = scaling = shearing = deformation = true; } + else if (OPTION("-add-identity")) { identity_value = NaN; } + else if (OPTION("-add-identity-with-weight")) { identity_value = atof(ARGUMENT); } + else if (OPTION("-add-identity-for-dofname")) { identity_name = ARGUMENT; } + else if (OPTION("-max-frechet-iterations")) frechet_iter = atoi(ARGUMENT); + else HANDLE_STANDARD_OR_UNKNOWN_OPTION(); + } + + if (sigma < .0) { + cerr << EXECNAME << ": Gaussian sigma value must be positive" << endl; + exit(1); + } + + // Parse input text file + if (doflist) { + if (dofin.size() > 0) { + cerr << EXECNAME << ": Text file specified for transformations but " << dofin.size() << " have already been read." << endl; + exit(1); + } + ifstream in(doflist); + if (!in) { + cout << EXECNAME << ": Unable to open file " << doflist << endl; + exit(1); + } + size_t lineno = 0; + while (in) { + lineno++; + dofin.push_back(""); + value.push_back(NaN); + in >> dofin.back() >> value.back(); + if (in.eof()) { + dofin.pop_back(); + value.pop_back(); + break; + } + if (dofin.back().empty() || dofin.back()[0] == '#') { + dofin.pop_back(); + value.pop_back(); + } else if (dofin.back() != identity_name) { + if (prefix) dofin.back().insert(0, prefix); + if (suffix) dofin.back() += suffix; + } + if (in.fail()) { + cout << EXECNAME << ": Failed to parse line " << lineno << " of file " << doflist << endl; + exit(1); + } + } + in.close(); + if (dofin.empty()) { + cerr << EXECNAME << ": No transformations listed in file " << doflist << endl; + exit(1); + } + if (verbose) { + cout << "Read " << dofin.size() << " transformations from file " << doflist << endl; + } + } else if (sigma != .0) { + cerr << EXECNAME << ": Input transformations must be specified using -dofnames when using -gaussian kernel." << endl; + exit(1); + } + + // Number of input transformations, including optional additional implicit identity transformation + const int N = static_cast(IsInf(identity_value) ? dofin.size() : dofin.size()+1u); + + // Compute actual weights + double *w = NULL; // actual weights + double W = .0; // sum of all weights + + if (value.size() > 0) { + // Add value for implicit identity transformation + if (!IsInf(identity_value)) value.push_back(identity_value); + w = new double[value.size()]; + for (size_t i = 0; i < value.size(); ++i) { + if (sigma > .0) { + w[i] = 1.0 / sqrt(2.0) / sigma; + if (!IsNaN(value[i])) w[i] *= exp(-pow((mean - value[i])/sigma, 2.0)/2.0); + } else { + if (IsNaN(value[i])) w[i] = 1.0; + else w[i] = value[i]; + } + if (w[i] < epsilon) w[i] = .0; + W += w[i]; + } + } + + // If no weights specified, all transformations contribute equally + if (W == .0) W = N; + + // Some verbose reporting + if (verbose > 1) { + int nzero = 0; + size_t width = 0; + double minw = (w ? w[0] : 1.0); + double maxw = (w ? w[0] : 1.0); + double dw = 1.0; + for (size_t i = 0; i < value.size(); ++i) { + if (w[i] < minw) minw = w[i]; + if (w[i] > maxw) maxw = w[i]; + if (w[i] == .0) { ++nzero; if (verbose < 2) continue; } + if (dofin[i].length() > width) width = dofin[i].length(); + for (size_t j = i+1; j < value.size(); ++j) { + double d = fabs(w[i] - w[j]); + if (d > 1e-6 && d < dw) dw = d; + } + } + if (value.size() > 0) { + bool *reported = new bool[value.size()]; + for (size_t i = 0; i < value.size(); ++i) reported[i] = (verbose < 2 && w[i] == .0); + for (double level = maxw; level >= minw; level -= dw) { + for (size_t i = 0; i < value.size(); ++i) { + if (reported[i] || w[i] < level) continue; + printf("%3d. Weight of %-*s = %.5f\n", int(i+1), int(width), dofin[i].c_str(), w[i]); + reported[i] = true; + } + } + delete[] reported; + printf("\n"); + } + printf("Number of transformations with zero weight: %d\n", nzero); + printf("Number of transformations with non-zero weight: %d\n", N - nzero); + printf("\n"); + printf("Total number of transformations: %d\n", N); + printf("Total weight of transformations: %.5f\n", W); + } + + // Make input file paths absolute + if (dofdir) { + for (size_t i = 0; i < dofin.size(); ++i) { + if (dofin[i] != identity_name && isRelativePath(dofin[i])) { + dofin[i] = string(dofdir) + '/' + dofin[i]; + } + } + } + + // (Common) type and attributes of local input/output transformations + string mtype; // mulit-level transformation type + string type; // common transformation type + ImageAttributes attr; // common FFD lattice attributes + + // Average displacement fields sampled in specified target image domain + if (!avgdofs && target_name) { + unique_ptr target(BaseImage::New(target_name)); + attr = target->GetImageAttributes(); + // Otherwise, determine common type of local input transformations + // if a mix of transformations is given, use attributes of first FFD + // to define the image domain of the average displacement field + } else { + if (verbose) cout << "Checking type of input transformations...", cout.flush(); + for (size_t i = 0; i < dofin.size(); ++i) { + if (dofin[i] == identity_name) continue; + unique_ptr t(Transformation::New(dofin[i].c_str())); + Transformation *p = t.get(); + MultiLevelTransformation *mffd = dynamic_cast(t.get()); + if (mffd) { + if (mffd->NumberOfLevels() == 0) p = mffd->GetGlobalTransformation(); + else if (mffd->NumberOfLevels() == 1) p = mffd->GetLocalTransformation(0); + else { + attr = mffd->GetLocalTransformation(-1)->Attributes(); + type.clear(); + break; + } + if (mtype.empty()) mtype = mffd->NameOfClass(); + else if (mtype != mffd->NameOfClass()) mtype.clear(); + } + FreeFormTransformation *ffd = dynamic_cast(p); + if (ffd) { + if (type.empty()) { + type = ffd->NameOfClass(); + attr = ffd->Attributes(); + } else if (type != ffd->NameOfClass() || ffd->Attributes() != attr) { + mtype.clear(); + type .clear(); + break; + } + } + } + if (verbose) cout << " done" << endl; + } + + // Cases in which local transformation parameters are averaged directly and stored + // using the very same transformation model as given as input which does not require + // a conversion of the input transformations to dense displacement fields + if (avgdofs && deformation) { + if (attr.NumberOfLatticePoints() > 0 && // input contains FFDs to be averaged + type != "BSplineFreeFormTransformationSV" && + type != "LinearFreeFormTransformation3D" && + type != "BSplineFreeFormTransformation3D") { + if (verbose) cout << endl; + cerr << EXECNAME << ": Cannot average local input transformation parameters directly (-dofs option)." << endl; + cerr << endl; + cerr << " Transformations must therefore either all be linear transformations or transformations with local" << endl; + cerr << " component of type \"Linear FFD\", \"B-spline FFD\", or \"SV B-spline FFD\" (see dofprint output)." << endl; + exit(1); + } + if (logspace == 1) { + if (type != "BSplineFreeFormTransformationSV") { + if (verbose) cout << endl; + cerr << EXECNAME << ": Combining options -log and -dofs only possible for SV B-spline FFDs! Use only one of them." << endl; + exit(1); + } + } else if (logspace == 0) { + if (type != "LinearFreeFormTransformation" && type != "BSplineFreeFormTransformation3D") { + if (verbose) cout << endl; + cerr << EXECNAME << ": Combining options -nolog and -dofs only possible for Linear/B-spline FFDs! Use only one of them." << endl; + exit(1); + } + } + } + + if (verbose) { + cout << "Average deformation: " << ToString(deformation) << endl; + cout << "Average DoFs directly: " << ToString(avgdofs) << endl; + cout << "Average in log-space: " << ToString(logspace) << endl; + } + + // Initialize intermediate data structures + Matrix *A = NULL; + GenericImage d, avgD; + + if (translation || rotation || scaling || shearing) { + A = new Matrix[N]; + for (int i = 0; i < N; ++i) { + A[i].Initialize(4, 4); + A[i].Ident(); + } + } + if (deformation && attr.NumberOfLatticePoints() > 0) { + d .Initialize(attr, 3); + avgD.Initialize(attr, 3); + } + + // Read parameters of input transformations + for (size_t i = 0; i < dofin.size(); ++i) { + if (dofin[i] == identity_name) continue; + // Read transformation from file + unique_ptr t(Transformation::New(dofin[i].c_str())); + // Determine actual type of transformation + HomogeneousTransformation *global = NULL; + RigidTransformation *rigid = NULL; + SimilarityTransformation *similarity = NULL; + AffineTransformation *affine = NULL; + FreeFormTransformation3D *ffd = NULL; + MultiLevelTransformation *mffd = NULL; + FluidFreeFormTransformation *fluid = NULL; + if (!((fluid = dynamic_cast(t.get())) || + (mffd = dynamic_cast(t.get())) || + (ffd = dynamic_cast(t.get())) || + (affine = dynamic_cast(t.get())) || + (similarity = dynamic_cast(t.get())) || + (rigid = dynamic_cast(t.get())))) { + cerr << EXECNAME << ": Cannot process transformation \"" << dofin[i] << "\" of type " << t->NameOfClass() << endl; + exit(1); + } + // Set base class pointers for common interface access, + // e.g., rigid components can be modified using the common + // RigidTransformation interface even if the actual transformation + // is of type SimilarityTransformation, AffineTransformation, + // or even a subclass of MultiLevelTransformation. + if (fluid) mffd = fluid; + if (mffd) { + affine = mffd->GetGlobalTransformation(); + if (mffd->NumberOfLevels() > 0) { + FreeFormTransformation *affd = mffd->GetLocalTransformation(mffd->NumberOfLevels() - 1); + ffd = dynamic_cast(affd); + if (!ffd) { + cerr << EXECNAME << ": Cannot process MFFD \"" << dofin[i] << "\" with level of type " << affd->NameOfClass() << endl; + exit(1); + } + } + } + if (affine ) similarity = affine; + if (similarity) rigid = similarity; + if (rigid ) global = rigid; + // Get local transformation parameters + if (deformation && ffd) { + if (invert) { + cerr << EXECNAME << ": -invert option only supported for rigid/affine transformations" << endl; + exit(1); + } + // Get inverse of the linear transformation matrix + Matrix invA(4, 4); + if (global) { + invA = global->GetMatrix(); + invA.Invert(); + } else invA.Ident(); + // Average parameters at control points directly when possible + double x, y, z; + if (avgdofs) { + for (int k = 0; k < attr._z; ++k) { + for (int j = 0; j < attr._y; ++j) { + for (int i = 0; i < attr._x; ++i) { + // Get control point parameters + ffd->Get(i, j, k, x, y, z); + // Remove dependency on global transformation + // Note: Applies also when the FFD parameters are stationary velocities, see below. + d(i, j, k, 0) = invA(0, 0) * x + invA(0, 1) * y + invA(0, 2) * z; + d(i, j, k, 1) = invA(1, 0) * x + invA(1, 1) * y + invA(1, 2) * z; + d(i, j, k, 2) = invA(2, 0) * x + invA(2, 1) * y + invA(2, 2) * z; + } + } + } + // Otherwise,... + } else { + // Get local displacement field + Matrix i2w = attr.GetImageToWorldMatrix(); + for (int k = 0; k < attr._z; ++k) { + for (int j = 0; j < attr._y; ++j) { + for (int i = 0; i < attr._x; ++i) { + // Convert voxel indices to world coordinates + x = i2w(0, 0) * i + i2w(0, 1) * j + i2w(0, 2) * k + i2w(0, 3); + y = i2w(1, 0) * i + i2w(1, 1) * j + i2w(1, 2) * k + i2w(1, 3); + z = i2w(2, 0) * i + i2w(2, 1) * j + i2w(2, 2) * k + i2w(2, 3); + // Evaluate (total) local voxel displacement + t->LocalDisplacement(x, y, z); + // Remove dependency on global transformation + d(i, j, k, 0) = invA(0, 0) * x + invA(0, 1) * y + invA(0, 2) * z; + d(i, j, k, 1) = invA(1, 0) * x + invA(1, 1) * y + invA(1, 2) * z; + d(i, j, k, 2) = invA(2, 0) * x + invA(2, 1) * y + invA(2, 2) * z; + } + } + } + // Compute stationary velocity field + if (logspace) { + DisplacementToVelocityFieldBCH dtov; + dtov.Input (&d); + dtov.Output(&d); + dtov.Run(); + // Smooth velocities if only few transformations are being averaged, + // otherwise rely on the average velocity field to be sufficiently smooth + if (N < 5) { + GaussianBlurring blur(max(attr._dx, max(attr._dy, attr._dz))); + blur.Input (&d); + blur.Output(&d); + blur.Run(); + } + } + } + // Add to sum of local transformations + if (w) d *= w[i]; + avgD += d; + } + // Get global transformation parameters + // (**after** local parameters as the following modifies the transformation) + if (A && global) { + if (invert) global->Invert(); + if (!translation) { + rigid->PutTranslationX(0); + rigid->PutTranslationY(0); + rigid->PutTranslationZ(0); + } + if (!rotation) { + rigid->PutRotationX(0); + rigid->PutRotationY(0); + rigid->PutRotationZ(0); + } + if (!scaling) { + if (affine) { + affine->PutScaleX(100); + affine->PutScaleY(100); + affine->PutScaleZ(100); + } else if (similarity) { + similarity->PutScale(100); + } + } + if (!shearing && affine) { + affine->PutShearXY(0); + affine->PutShearXZ(0); + affine->PutShearYZ(0); + } + A[i] = global->GetMatrix(); + } + } + + // Compute average global transformation + AffineTransformation globalAvg; + globalAvg.AllowTranslations(translation); + globalAvg.AllowRotations(rotation); + globalAvg.AllowScaling(scaling); + globalAvg.AllowShearing(shearing); + + if (A) { + if (verbose) cout << "\nComputing average global transformation..." << endl; + Matrix avgA; + if (biinvariant) avgA = BiInvariantMean (N, A, w, frechet_iter); + else avgA = LogEuclideanMean(N, A, w); + if (invavg) avgA.Invert(); + if (approxglobal) globalAvg.ApproximateAsNew(avgA); + else globalAvg.PutMatrix(avgA); + if (verbose) { + Indent indent(1); + if (verbose > 1) { + cout << "\n"; + avgA.Print(indent); + } + const Matrix avgAInv = avgA.Inverse(); + Matrix residual(4, 4); + for (int i = 0; i < N; ++i) { + if (w[i] == .0) continue; + residual += (avgAInv * A[i]).Log() * w[i]; + } + residual /= W; + cout << "\n" << indent << "Deviation from barycentric equation = " << residual.Norm() << endl; + residual.Ident(); + residual -= avgA * globalAvg.GetMatrix().Inverse(); + cout << indent << "Norm of (A * inv(global._matrix) - I) = " << residual.Norm() << "\n" << endl; + globalAvg.Print(indent + 1); + } + if (verbose) cout << "\nComputing average global transformation... done" << endl; + } + + // If global average only... + if (avgD.IsEmpty()) { + + // Write average affine transformation + globalAvg.Write(dofout.c_str()); + + // Otherwise,... + } else { + + if (invavg) { + cerr << EXECNAME << ": -inverse[-dofs] option only supported for global transformations" << endl; + exit(1); + } + + // Divide sum of local transformations by total weight to get average + avgD /= W; + // Convert velocities back to displacements, unless the output transformation + // is parameterized by a stationary velocity field + if (logspace && !avgdofs) { + VelocityToDisplacementFieldSS vtod; + vtod.Input (&avgD); + vtod.Output(&avgD); + vtod.Run(); + } + // Add dependency on average global transformation such that + // avgT = avgA (x + avgD) = avgA x + avgA avgD = avgT_global + avgT_local + // + // Note that also in case of avgD being a stationary velocity field (avgdofs == true), + // we can pre-multiply it the same way with the global transformation matrix as in + // case of a displacement field. This can be illustrated by looking at the formula + // of the family of explicit Runge-Kutta methods for the exponentiation step which + // consists of a sum of the velocities at different time steps, i.e., + // avgT(x) = avgA x + avgA sum_i b_i h v(x_i) = avgA x + sum_i b_i h (avgA v(x_i)). + double x, y, z; + const Matrix &avgA = globalAvg.GetMatrix(); + for (int k = 0; k < attr._z; ++k) { + for (int j = 0; j < attr._y; ++j) { + for (int i = 0; i < attr._x; ++i) { + x = avgD(i, j, k, 0), y = avgD(i, j, k, 1), z = avgD(i, j, k, 2); + avgD(i, j, k, 0) = avgA(0, 0) * x + avgA(0, 1) * y + avgA(0, 2) * z; + avgD(i, j, k, 1) = avgA(1, 0) * x + avgA(1, 1) * y + avgA(1, 2) * z; + avgD(i, j, k, 2) = avgA(2, 0) * x + avgA(2, 1) * y + avgA(2, 2) * z; + } + } + } + // Construct FFD from average deformation + FreeFormTransformation *ffd = NULL; + if (avgdofs) { + if (type == "LinearFreeFormTransformation3D") { + ffd = new LinearFreeFormTransformation3D(avgD, true); + } else if (type == "BSplineFreeFormTransformation3D") { + ffd = new BSplineFreeFormTransformation3D(avgD, true); + } else if (type == "BSplineFreeFormTransformationSV") { + ffd = new BSplineFreeFormTransformationSV(avgD, true); + } + } + if (!ffd) { + if (bsplineffd) { + ConvertToCubicBSplineCoefficients(avgD, 0); + ConvertToCubicBSplineCoefficients(avgD, 1); + ConvertToCubicBSplineCoefficients(avgD, 2); + ffd = new BSplineFreeFormTransformation3D(avgD, true); + } else { + ffd = new LinearFreeFormTransformation3D(avgD, true); + } + } + // Write average non-linear transformation + MultiLevelTransformation *mffd = NULL; + if (mtype == "MultiLevelStationaryVelocityTransformation") { + mffd = new MultiLevelStationaryVelocityTransformation(globalAvg); + } else { + mffd = new MultiLevelFreeFormTransformation(globalAvg); + } + mffd->PushLocalTransformation(ffd); + mffd->Write(dofout.c_str()); + delete mffd; + } + + // Clean up + delete[] w; + delete[] A; + + return 0; +} diff --git a/Applications/average-images.cc b/Applications/average-images.cc new file mode 100644 index 000000000..c2dcf2701 --- /dev/null +++ b/Applications/average-images.cc @@ -0,0 +1,537 @@ +/* + * Medical Image Registration ToolKit (MIRTK) + * + * Copyright 2013-2015 Imperial College London + * Copyright 2013-2015 Andreas Schuh + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include + +#include +#include +#include + +#ifdef HAVE_MIRTK_Transformation +# include +# include +# include +#endif + +using namespace mirtk; + + +// ============================================================================= +// Help +// ============================================================================= + +/// Print program usage information +void PrintHelp(const char* name) +{ + cout << endl; + cout << "Usage: " << name << " -images [options]" << endl; + cout << " " << name << " -image [-dof ] -image [-dof ]... [options]" << endl; + cout << " " << name << " [] []... [options]" << endl; + cout << " " << name << " [options]" << endl; + cout << endl; + cout << "Description:" << endl; + cout << " Computes voxelwise average image of the given (transformed) input images." << endl; +#ifdef HAVE_MIRTK_Transformation + cout << endl; + cout << " The input images are optionally transformed by an image-to-average space transformation" << endl; + cout << " (:option:`-dof`) or its inverse (:option:`-dof_i`)." << endl; +#endif + cout << endl; + cout << "Required arguments:" << endl; + cout << " Voxel-wise average image." << endl; + cout << " -images Text file with N lines containing the file name of each input image," << endl; + cout << " optionally a transformation file name (see :option:`-dof` and :option:`-dof_i`)," << endl; + cout << " optionally followed by a weight for a weighted output average (default weight is 1)." << endl; + cout << " The first line of the text file must specify the common base directory" << endl; + cout << " of all relative image and transformation file paths occurring on the" << endl; + cout << " subsequent N lines. A path starting with './' must be relative to the" << endl; + cout << " directory containing the input text file itself." << endl; + cout << " -image A single input image." << endl; +#ifdef HAVE_MIRTK_Transformation + cout << " -dof Specifies a transformation to be applied to the preceeding image. (default: none)" << endl; + cout << " -dof_i Specifies a transformation whose inverse is to be applied to the" << endl; + cout << " preceeding image. (default: none)" << endl; +#endif + cout << endl; + cout << "Optional arguments:" << endl; + cout << " -size [ []] Voxel size of intensity average image. (default: average of input images)" << endl; + cout << " -padding Input padding and output background value. No input padding if not specified." << endl; + cout << " Output background value zero by default or minimum average intensity minus 1. (default: 0)" << endl; + cout << " -interp Interpolation mode, e.g., NN, Linear, BSpline, Cubic, Sinc. (default: Linear)" << endl; + cout << " -label Segmentation label of which to create an average probability map." << endl; + PrintCommonOptions(cout); + cout << endl; +} + +// ============================================================================= +// Types +// ============================================================================= + +// Type of input images +typedef RealPixel InputType; +typedef GenericImage InputImage; +typedef GenericInterpolateImageFunction InputImageFunction; + +// Type of output image +typedef float OutputType; +typedef GenericImage OutputImage; + +// ============================================================================= +// Auxiliary functions +// ============================================================================= + +/// Read image file names and times from input list file +int read_image_list_file(const char *image_list_name, + Array &names, + Array &dofs, + Array &weights) +{ + string base_dir = Directory(image_list_name); // Base directory containing image files + ifstream iff(image_list_name); // List input file stream + string line; // Input line + string part; // Line part + string entry; // Column entry + bool inquote = false; // Parsing parts of a quoted string + int l = 0; // Line number + + names .clear(); + dofs .clear(); + weights.clear(); + + // Read base directory for image files + if (!getline(iff, line)) { + cerr << "Error: Cannot parse list file " << image_list_name << endl; + return 0; + } + + if (!base_dir.empty() && line[0] != PATHSEP) { + if (line != ".") base_dir += PATHSEP + line; + } else { + base_dir = line; + } + l++; + + while (getline(iff, line)) { + l++; + // Ignore comment lines starting with # character + if (line[0] == '#') continue; + // Split line at space + istringstream lss(line); + int c = 0; + while (getline(lss, part, ' ')) { + if (inquote) entry += ' '; + if (part.empty()) continue; + if (part[0] == '"') { + part.erase(0, 1); + inquote = !inquote; + } + if (part[part.size()-1] == '"') { + part.erase(part.size()-1, 1); + inquote = !inquote; + } + entry += part; + if (inquote) continue; + // Add column entry to the proper output container + if (c == 0) { + if (!base_dir.empty() && entry[0] != PATHSEP) entry = base_dir + PATHSEP + entry; + names.push_back(entry); + } else if (c == 1 && !isdigit(entry[0])) { + if (!base_dir.empty() && entry[0] != PATHSEP) entry = base_dir + PATHSEP + entry; + dofs.push_back(entry); + } else if (c == 2 || (c == 1 && isdigit(entry[0]))) { + if (c == 1) dofs.push_back(""); + istringstream iss(entry); + OutputType weight; + iss >> weight; + if (!iss) { + cerr << "Error: " << image_list_name << ":" << l << ": Failed to parse weight!" << endl; + names .clear(); + weights.clear(); + return 0; + } + weights.push_back(weight); + } else { + cerr << "Error: " << image_list_name << ":" << l << ": Too many entries!"<< endl; + names .clear(); + weights.clear(); + return 0; + } + // Parse next column + entry.clear(); + c++; + } + if (c == 0) continue; // Skip empty lines + // Add default for non-specified columns + dofs .resize(names.size()); + weights.resize(names.size(), 1.0); + } + + return static_cast(names.size()); +} + +// ----------------------------------------------------------------------------- +#ifdef HAVE_MIRTK_Transformation +void Read(const string &image_name, const string &imdof_name, + InputImage &image, Transformation *&imdof, + bool imdof_invert = false) +{ + image.Read(image_name.c_str()); + if (imdof_name.empty()) { + imdof = NULL; + } else { + imdof = Transformation::New(imdof_name.c_str()); + HomogeneousTransformation *lin = dynamic_cast(imdof); + FluidFreeFormTransformation *fluid = dynamic_cast(imdof); + if (fluid && imdof_invert) lin = fluid->GetGlobalTransformation(); + if (lin) { + Matrix mat = lin->GetMatrix(); + if (imdof_invert) mat.Invert(); + image.PutAffineMatrix(mat, true); + Delete(imdof); + } + if (fluid && lin) fluid->GetGlobalTransformation()->Reset(); + } +} +#endif // HAVE_MIRTK_Transformation + +// ----------------------------------------------------------------------------- +struct AddVoxelValueToAverage : public VoxelFunction +{ + OutputImage *_Average; + InputImageFunction *_Image; + double _Weight; + int _Label; + + #ifdef HAVE_MIRTK_Transformation + Transformation *_Transformation; + GenericImage *_Displacement; + bool _Invert; + #endif + + void operator()(int i, int j, int k, int, OutputType *avg) + { + double x = i, y = j, z = k; + _Average->ImageToWorld(x, y, z); + #ifdef HAVE_MIRTK_Transformation + if (_Displacement) { + x += _Displacement->Get(i, j, k, 0); + y += _Displacement->Get(i, j, k, 1); + z += _Displacement->Get(i, j, k, 2); + } else if (_Transformation) { + // Note: Input transformation is from image to average! + if (_Invert) _Transformation->Transform(x, y, z); + else _Transformation->Inverse (x, y, z); + } + #endif // HAVE_MIRTK_Transformation + _Image->Input()->WorldToImage(x, y, z); + if (_Label > 0) { + if (static_cast(_Image->Evaluate(x, y, z)) == _Label) { + *avg += static_cast(_Weight); + } + } else { + const OutputType bg = static_cast(_Average->GetBackgroundValueAsDouble()); + if ((IsNaN(bg) && IsNaN(*avg)) || (*avg == bg)) { + *avg = static_cast(_Weight * _Image->Evaluate(x, y, z)); + } else { + *avg += static_cast(_Weight * _Image->Evaluate(x, y, z)); + } + } + } +}; + +// ----------------------------------------------------------------------------- +void Add(OutputImage &average, InputImage &image, int label = -1, + #ifdef HAVE_MIRTK_Transformation + Transformation *transformation = NULL, + bool invert = false, + #endif // HAVE_MIRTK_Transformation + OutputType weight = 1.0, + InterpolationMode interpolation = Interpolation_Linear) +{ + if (label > 0) interpolation = Interpolation_NN; + AddVoxelValueToAverage add; + #ifdef HAVE_MIRTK_Transformation + GenericImage disp; + if (transformation && transformation->RequiresCachingOfDisplacements()) { + disp.Initialize(average.Attributes(), 3); + // Note: Input transformation is from image to average! + if (invert) transformation->Displacement(disp); + else transformation->InverseDisplacement(disp); + } + add._Transformation = transformation; + add._Displacement = (disp.IsEmpty() ? NULL : &disp); + add._Invert = invert; + #endif // HAVE_MIRTK_Transformation + add._Average = &average; + add._Weight = weight; + add._Label = label; + add._Image = InputImageFunction::New(interpolation, &image); + add._Image->Initialize(); + ParallelForEachVoxel(average.Attributes(), average, add); + delete add._Image; +} + +// ============================================================================= +// Main +// ============================================================================= + +// ----------------------------------------------------------------------------- +int main(int argc, char **argv) +{ + InputImage sequence, image; + Array image_name; + Array image_weight; + int nimages; + + #ifdef HAVE_MIRTK_Transformation + Transformation *imdof = NULL; + Array imdof_name; + Array imdof_invert; + #endif // HAVE_MIRTK_Transformation + + // --------------------------------------------------------------------------- + // Command-line parsing + + // Required positional argument(s) + REQUIRES_POSARGS(1); + const char *output_name = POSARG(1); + + // Optional positional argument(s) + int nposarg; + for (nposarg = NUM_POSARGS + 1; nposarg < argc; ++nposarg) { + bool img = (strcmp(argv[nposarg], "-image") == 0); + bool inv = (strcmp(argv[nposarg], "-dof_i") == 0); + bool dof = inv || (strcmp(argv[nposarg], "-dof") == 0); + if (!img && !dof) { + if (argv[nposarg][0] == '-') break; + int l = strlen(argv[nposarg]); + dof = ((l > 4 && strcmp(&argv[nposarg][l-4], ".dof") == 0) || + (l > 7 && strcmp(&argv[nposarg][l-7], ".dof.gz") == 0)); + } else { + ++nposarg; + if (nposarg == argc) { + PrintHelp(EXECNAME); + cout << endl; + FatalError("Option " << argv[nposarg-1] << " requires an argument"); + } + } + if (dof) { + #ifdef HAVE_MIRTK_Transformation + imdof_name .resize(image_name.size()); + imdof_invert.resize(image_name.size()); + imdof_name .back() = argv[nposarg]; + imdof_invert.back() = inv; + #else // HAVE_MIRTK_Transformation + FatalError("Cannot apply an image transformation, rebuild with Transformation module enabled!"); + #endif // HAVE_MIRTK_Transformation + } else { + image_name.push_back(argv[nposarg]); + } + } + --nposarg; + + nimages = static_cast(image_name.size()); + image_weight.resize(nimages, 1.0); + #ifdef HAVE_MIRTK_Transformation + imdof_name .resize(nimages); + imdof_invert.resize(nimages); + #endif // HAVE_MIRTK_Transformation + + // Parse arguments + const char *image_list_name = NULL; + const char *reference_name = NULL; + double padding = numeric_limits::quiet_NaN(); + InterpolationMode interpolation = Interpolation_Linear; + bool voxelwise = false; + int label = -1; + int margin = -1; + double dx = .0, dy = .0, dz = .0; + + for (ARGUMENTS_AFTER(nposarg)) { + if (OPTION("-images" )) image_list_name = ARGUMENT; + else if (OPTION("-reference")) reference_name = ARGUMENT; + else if (OPTION("-size")) { + PARSE_ARGUMENT(dx); + dy = dz = dx; + if (HAS_ARGUMENT) { + PARSE_ARGUMENT(dy); + if (HAS_ARGUMENT) PARSE_ARGUMENT(dz); + else dz = .0; + } + } + else if (OPTION("-padding")) PARSE_ARGUMENT(padding); + else if (OPTION("-voxelwise")) voxelwise = true; + else if (OPTION("-margin")) PARSE_ARGUMENT(margin); + else if (OPTION("-label")) PARSE_ARGUMENT(label); + else if (OPTION("-interp")) PARSE_ARGUMENT(interpolation); + else HANDLE_COMMON_OR_UNKNOWN_OPTION(); + } + + if (label > 0 && IsNaN(padding)) padding = .0; + + // --------------------------------------------------------------------------- + // Collect (further) input image meta-data... + + // ...from files specified as positional arguments + if (image_name.size() == 1) { + if (verbose) cout << "Reading sequence ... ", cout.flush(); + sequence.Read(image_name.front().c_str()); + nimages = sequence.T(); + if (nimages < 2) { + if (verbose) cout << endl; + FatalError("Input sequence contains only one temporal frame!"); + } + if (verbose) cout << " done\n" << endl; + } + + // ...from files specified in list file + if (image_list_name) { + Array names; + Array imdofs; + Array weights; + int n = read_image_list_file(image_list_name, names, imdofs, weights); + if (n == 0) { + if (verbose) cout << endl; + FatalError("Failed to parse input list file " << image_list_name << "!"); + } + if (verbose) cout << "Found " << n << " images in list file " << image_list_name << endl; + nimages += n; + image_name .insert(image_name .end(), names .begin(), names .end()); + image_weight.insert(image_weight.end(), weights.begin(), weights.end()); + imdof_name .insert(imdof_name .end(), imdofs .begin(), imdofs .end()); + imdof_invert.resize(nimages, false); + } + + // Check that any input image is given + if (nimages < 1) { + PrintHelp(EXECNAME); + FatalError("No input image(s) specified!"); + } + + // Compute sum of weights + OutputType wsum = .0; + for (int n = 0; n < nimages; ++n) wsum += image_weight[n]; + + // Compute voxel-wise average of (interpolated) input intensities on + // common (mapped) image domain of all input images; for example, + // used for anatomical intensity atlas construction + + // Compute field-of-view which contains all images + if (verbose) cout << "Determine field-of-view which contains all images...", cout.flush(); + ImageAttributes fov; + if (reference_name) { + GreyImage reference(reference_name); + fov = reference.Attributes(); + } else { + Array attr; + for (int n = 0; n < nimages; ++n) { + #ifdef HAVE_MIRTK_Transformation + Read(image_name[n], imdof_name[n], image, imdof, imdof_invert[n]); + Delete(imdof); + #else // HAVE_MIRTK_Transformation + image.Read(image_name[n].c_str()); + #endif // HAVE_MIRTK_Transformation + if (image.T() > 1) { + if (verbose) cout << " failed" << endl; + FatalError("Image " << (n+1) << " has four dimensions!"); + } + attr.push_back(OrthogonalFieldOfView(image.Attributes())); + } + fov = OverallFieldOfView(attr); + } + if (verbose) cout << " done" << endl; + + // Set desired output voxel size + if (dx > .0) { + fov._x = ceil(fov._x * fov._dx / dx); + fov._dx = dx; + } + if (dy > .0) { + fov._y = ceil(fov._y * fov._dy / dy); + fov._dy = dy; + } + if (dz > .0) { + fov._z = ceil(fov._z * fov._dz / dz); + fov._dz = dz; + } + + // Compute average image + OutputImage average(fov); + average = padding; + average.PutBackgroundValueAsDouble(padding); + for (int n = 0; n < nimages; ++n) { + if (sequence.IsEmpty()) { + if (verbose) { + cout << "Add image " << setw(3) << (n+1) << " out of " << nimages << "... "; + cout.flush(); + } + #ifdef HAVE_MIRTK_Transformation + Read(image_name[n], imdof_name[n], image, imdof, imdof_invert[n]); + #else // HAVE_MIRTK_Transformation + image.Read(image_name[n].c_str()); + #endif // HAVE_MIRTK_Transformation + } else { + if (verbose) { + cout << "Add frame " << setw(3) << (n+1) << " out of " << nimages << "... "; + cout.flush(); + } + sequence.GetFrame(image, n); + #ifdef HAVE_MIRTK_Transformation + imdof = NULL; + #endif + } + + Add(average, image, label, + #ifdef HAVE_MIRTK_Transformation + imdof, imdof_invert[n], + #endif + image_weight[n], interpolation); + + #ifdef HAVE_MIRTK_Transformation + Delete(imdof); + #endif + + if (verbose) cout << " done" << endl; + } + if (wsum > .0) average /= wsum; + + // Crop/pad average image + if (margin >= 0) average.CropPad(margin); + + // Replace NaN background values + if (IsNaN(average.GetBackgroundValueAsDouble())) { + OutputImage::VoxelType min_value, max_value; + average.GetMinMax(min_value, max_value); + padding = min(.0, double(min_value) - 1.0); + for (int idx = 0; idx < average.NumberOfVoxels(); ++idx) { + if (IsNaN(average(idx))) { + average(idx) = voxel_cast(padding); + } + } + average.PutBackgroundValueAsDouble(padding); + } + + // Write average image + average.Write(output_name); + + + return 0; +} diff --git a/Applications/calculate-distance-map.cc b/Applications/calculate-distance-map.cc new file mode 100644 index 000000000..bc567ac79 --- /dev/null +++ b/Applications/calculate-distance-map.cc @@ -0,0 +1,266 @@ +/* + * Medical Image Registration ToolKit (MIRTK) + * + * Copyright 2008-2015 Imperial College London + * Copyright 2008-2013 Daniel Rueckert, Julia Schnabel + * Copyright 2013-2015 Andreas Schuh + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +#include +#include + +#include +#include + +using namespace mirtk; + + +// ============================================================================= +// Help +// ============================================================================= + +// ----------------------------------------------------------------------------- +void PrintHelp(const char *name) +{ + cout << "\n"; + cout << "Usage: " << name << " [options]\n"; + cout << "\n"; + cout << "Description:\n"; + cout << " Computes distance transformation of binary object mask.\n"; + cout << "\n"; + cout << "Arguments:\n"; + cout << " input Binary object mask.\n"; + cout << " output Distance transform.\n"; + cout << "\n"; + cout << "Optional arguments:\n"; + cout << " -distance Name of the distance transform. (default: euclidean)\n"; + cout << " - ``cityblock``: Manhatten city block (L1) distance.\n"; + cout << " - ``euclidean``: Euclidean distance.\n"; + cout << "\n"; + cout << "Euclidean distance transform options:\n"; + cout << " -2D Compute distance transform in 2D.\n"; + cout << " -3D Compute distance transform in 3D. (default)\n"; + cout << " -radial Compute radial distance transform. (default: off)\n"; + cout << " -tradial Compute radial distance transform in t dimension only.\n"; + cout << " -isotropic [] Resample distance transform to isotropic voxel size,\n"; + cout << " where the voxel size is s times the minimum voxel size.\n"; + cout << " When no argument given, s defaults to 1. (default: 0/off)\n"; + PrintStandardOptions(cout); + cout << "\n"; +} + +// ============================================================================= +// Auxiliaries +// ============================================================================= + +namespace mirtk { + + +/// Type of distance transform +enum DistanceTransformType +{ + DT_Unknown, + DT_Euclidean, + DT_CityBlock, + DT_Last +}; + +// ----------------------------------------------------------------------------- +/// Convert distance transform type enumeration value to string +template <> +inline string ToString(const DistanceTransformType &value, int w, char c, bool left) +{ + const char *str; + switch (value) { + case DT_Euclidean: str = "euclidean"; break; + case DT_CityBlock: str = "cityblock"; break; + default: str = "unknown"; break; + } + return ToString(str, w, c, left); +} + +// ----------------------------------------------------------------------------- +/// Convert string to distance transform type enumeration value +template <> +inline bool FromString(const char *str, DistanceTransformType &value) +{ + const string lstr = ToLower(str); + value = static_cast(DT_Last - 1); + while (value != DT_Unknown) { + if (ToString(value) == lstr) break; + value = static_cast(value - 1); + } + return value != DT_Unknown; +} + + +} // namespace mirtk + +// ============================================================================= +// Main +// ============================================================================= + +// ----------------------------------------------------------------------------- +int main(int argc, char *argv[]) +{ + typedef EuclideanDistanceTransform EuclideanDistanceTransformType; + typedef CityBlockDistanceTransform CityBlockDistanceTransformType; + + EXPECTS_POSARGS(2); + + const char *input_name = POSARG(1); + const char *output_name = POSARG(2); + + DistanceTransformType type = DT_Euclidean; + + EuclideanDistanceTransformType::Mode euclidean_mode; + euclidean_mode = EuclideanDistanceTransformType::DT_3D; + + int radial = 0; + double isotropic = .0; + + for (ALL_OPTIONS) { + if (OPTION("-distance") || OPTION("-mode")) { + PARSE_ARGUMENT(type); + } + else if (OPTION("-2D")) { + euclidean_mode = EuclideanDistanceTransformType::DT_2D; + } + else if (OPTION("-3D")) { + euclidean_mode = EuclideanDistanceTransformType::DT_2D; + } + else if (OPTION("-radial")) { + radial = 1; + } + else if (OPTION("-tradial")) { + radial = 2; + } + else if (OPTION("-isotropic")) { + if (HAS_ARGUMENT) PARSE_ARGUMENT(isotropic); + else isotropic = 1.0; + } + } + + // No linear isotropic resampling for city block distance + if (type == DT_CityBlock) isotropic = .0; + + // Read object image + if (verbose) cout << "Reading image " << input_name << "...", cout.flush(); + RealImage image(input_name); + if (verbose) cout << " done" << endl; + + // Compute distance transform + if (verbose) cout << "Computing " << ToString(type) << " distance transform..." << endl; + RealImage dmap; + + switch (type) { + + // City Block distance transform + case DT_CityBlock: { + CityBlockDistanceTransformType cbdt; + RealImage &outputA = dmap, outputB; + + if (verbose) cout << " Computing distance transform of interior..." << endl; + + cbdt.Input(&image); + cbdt.Output(&outputA); + cbdt.Run(); + + if (verbose) cout << " done\n Computing distance transform of exterior..." << endl; + + for (int idx = 0; idx < image.NumberOfVoxels(); ++idx) { + image(idx) = (image(idx) > .0f ? 0 : 1); + } + + cbdt.Output(&outputB); + cbdt.Run(); + + if (verbose) cout << " done" << endl; + + for (int idx = 0; idx < image.NumberOfVoxels(); ++idx) { + dmap(idx) = outputB(idx) - outputA(idx); + } + + } break; + + // Euclidean distance transform + case DT_Euclidean: { + EuclideanDistanceTransformType edt(euclidean_mode); + RealImage inputA, inputB, &outputA = dmap, outputB; + + for (int idx = 0; idx < image.NumberOfVoxels(); ++idx) { + if (image(idx) > .5) { + inputA(idx) = 1.0; + inputB(idx) = 0.0; + } else { + inputA(idx) = 0.0; + inputB(idx) = 1.0; + } + } + + if (verbose) cout << " Computing distance transform of interior..." << endl; + + edt.Input (& inputA); + edt.Output(&outputA); + edt.Run(); + + if (verbose) cout << " done\n Computing distance transform of exterior..." << endl; + + edt.Input (& inputB); + edt.Output(&outputB); + edt.Run(); + + if (verbose) cout << " done" << endl; + + for (int idx = 0; idx < image.NumberOfVoxels(); ++idx) { + dmap(idx) = sqrt(outputA(idx)) - sqrt(outputA(idx)); + } + + if (radial > 0) { + edt.Input (&dmap); + edt.Output(&dmap); + switch (radial) { + case 1: edt.Radial(); break; + case 2: edt.TRadial(); break; + } + } + } + + default: + FatalError("Unknown distance transform type: " << type); + } + + // Resample image to isotropic voxels (smallest voxel dimension) + if (isotropic > .0) { + double ds = min(min(dmap.GetXSize(), dmap.GetYSize()), dmap.GetZSize()) * isotropic; + if (verbose) cout << "Resampling distance map to isotropic voxel size: " << ds << " [mm]...", cout.flush(); + Resampling resampling(ds, ds, ds); + GenericLinearInterpolateImageFunction > interpolator; + resampling.Interpolator(&interpolator); + resampling.Input (&dmap); + resampling.Output(&dmap); + resampling.Run(); + if (verbose) cout << " done" << endl; + } + + if (verbose) cout << "Writing distance transform to " << output_name << "..."; + dmap.Write(output_name); + if (verbose) cout << " done" << endl; + + return 0; +} diff --git a/Applications/calculate-exponential-map.cc b/Applications/calculate-exponential-map.cc new file mode 100644 index 000000000..998f34a4a --- /dev/null +++ b/Applications/calculate-exponential-map.cc @@ -0,0 +1,195 @@ +/* + * Medical Image Registration ToolKit (MIRTK) + * + * Copyright 2013-2015 Imperial College London + * Copyright 2013-2015 Andreas Schuh + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +#include +#include +#include + +using namespace mirtk; + + +// ============================================================================= +// Global variables with default arguments +// ============================================================================= + +char IntegrationMethod[32] = {"SS"}; +int NumberOfSteps = 128; + +// ============================================================================= +// Help +// ============================================================================= + +// ----------------------------------------------------------------------------- +void PrintHelp(const char *name) +{ + cout << endl; + cout << "Usage: " << name << " [options]" << endl; + cout << " " << name << " [] [] [options]" << endl; + cout << endl; + cout << "Description:" << endl; + cout << " Reads a dense 3D (2D) stationary velocity field, computes the corresponding" << endl; + cout << " dense 3D (2D) displacement field, and writes the resulting vector field." << endl; + cout << endl; + cout << "Options:" << endl; + cout << " -ss Use scaling and squaring (SS) method. (default)" << endl; + cout << " -euler Use forward Euler method." << endl; + cout << " -steps Number of integration steps (2^n in case of SS). " + << "(default: " << NumberOfSteps << ")" << endl; + PrintCommonOptions(cout); + cout << endl; +} + +// ============================================================================= +// Main +// ============================================================================= + +// ----------------------------------------------------------------------------- +int main(int argc, char **argv) +{ + // Print help or version and exit if requested + HANDLE_HELP_OR_VERSION(); + + // Parse positional arguments + const char *velo_fname = NULL; + const char *disp_fname = NULL; + const char *vx_fname = NULL; + const char *vy_fname = NULL; + const char *vz_fname = NULL; + const char *dx_fname = NULL; + const char *dy_fname = NULL; + const char *dz_fname = NULL; + + int posargc = 0; + while (posargc+1 < argc && argv[posargc+1][0] != '-') posargc++; + + switch (posargc) { + case 2: + velo_fname = argv[1]; + disp_fname = argv[2]; + break; + case 4: + vx_fname = argv[1]; + vy_fname = argv[2]; + dx_fname = argv[3]; + dy_fname = argv[4]; + break; + case 6: + vx_fname = argv[1]; + vy_fname = argv[2]; + vz_fname = argv[3]; + dx_fname = argv[4]; + dy_fname = argv[5]; + dz_fname = argv[6]; + break; + default: + PrintHelp(EXECNAME); + exit(1); + } + + // Parse optional arguments + for (ARGUMENTS_AFTER(posargc)) { + if (OPTION("-int")) { + const char *arg = ARGUMENT; + unsigned int len = strlen(arg); + if (len > 99) len = 99; + for (unsigned int i = 0; i < len; i++) { + if ('a' <= arg[i] && arg[i] <= 'z') IntegrationMethod[i] = 'A' + (arg[i] - 'a'); + else IntegrationMethod[i] = arg[i]; + } + IntegrationMethod[len] = '\0'; + } else if (OPTION("-ss")) { + IntegrationMethod[0] = 'S'; + IntegrationMethod[1] = 'S'; + IntegrationMethod[2] = '\0'; + } else if (OPTION("-euler")) { + IntegrationMethod[0] = 'R'; + IntegrationMethod[1] = 'K'; + IntegrationMethod[2] = 'E'; + IntegrationMethod[3] = '1'; + IntegrationMethod[4] = '\0'; + } else if (OPTION("-steps")) { + NumberOfSteps = atoi(ARGUMENT); + } else { + HANDLE_COMMON_OR_UNKNOWN_OPTION(); + } + } + + // Read input velocities + unique_ptr v; + + if (velo_fname) { + v.reset(new RealImage(velo_fname)); + } else { + RealImage vx(vx_fname); + RealImage vy(vy_fname); + unique_ptr vz(vz_fname ? new RealImage(vz_fname) : nullptr); + + ImageAttributes attr = vx.Attributes(); + if (attr._t > 1) { + FatalError("Separate input component images must be scalar velocity fields.\n" + "Try calling this program with single vector field input and output images instead."); + } + if (vy.Attributes() != attr || (vz.get() && vz->Attributes() != attr)) { + FatalError("Image attributes of input velocity component images do not match!"); + } + attr._t = (vz.get() ? 3 : 2); + + v.reset(new RealImage(attr)); + for (int k = 0; k < attr._z; ++k) + for (int j = 0; j < attr._y; ++j) + for (int i = 0; i < attr._x; ++i) { + v->Put(i, j, k, 0, vx.Get(i, j, k)); + v->Put(i, j, k, 1, vy.Get(i, j, k)); + if (vz.get()) v->Put(i, j, k, 2, vz->Get(i, j, k)); + } + } + + // Instantiate filter which implements the desired integration method + unique_ptr > vtod; + + if (strcmp(IntegrationMethod, "SS") == 0) { + vtod.reset(new VelocityToDisplacementFieldSS()); + } else if (strcmp(IntegrationMethod, "RKE1") == 0 || strcmp(IntegrationMethod, "FORWARDEULER") == 0 || strcmp(IntegrationMethod, "EULER") == 0) { + vtod.reset(new VelocityToDisplacementFieldEuler()); + } else { + FatalError("Unknown integration method: " << IntegrationMethod); + } + + // Compute displacement field + RealImage d; + + vtod->Input(v.get()); + vtod->Output(&d); + vtod->NumberOfSteps(NumberOfSteps); + vtod->Run(); + + // Save output velocities + if (disp_fname) { + d.Write(disp_fname); + } else { + d.GetFrame(0).Write(dx_fname); + d.GetFrame(1).Write(dy_fname); + if (dz_fname) d.GetFrame(2).Write(dz_fname); + } + + return 0; +} diff --git a/Applications/calculate-lie-bracket.cc b/Applications/calculate-lie-bracket.cc new file mode 100644 index 000000000..35d363d24 --- /dev/null +++ b/Applications/calculate-lie-bracket.cc @@ -0,0 +1,207 @@ +/* + * Medical Image Registration ToolKit (MIRTK) + * + * Copyright 2013-2015 Imperial College London + * Copyright 2013-2015 Andreas Schuh + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +#include +#include + +using namespace mirtk; + + +// =========================================================================== +// Help +// =========================================================================== + +// --------------------------------------------------------------------------- +void PrintHelp(const char *name) +{ + cout << endl; + cout << "Usage: " << name << " [options]" << endl; + cout << " " << name << " [] [] [] [options]" << endl; + cout << endl; + cout << "Description:" << endl; + cout << " Computes the Lie bracket of two vector fiels :math:`z = [x, y]`." << endl; + cout << endl; + cout << "Options:" << endl; + cout << " -jac Use Jacobian definition of the Lie bracket. (default)" << endl; + cout << " -dif Use difference of compositions: :math:`z = x(y) - y(x)`." << endl; + PrintCommonOptions(cout); + cout << endl; +} + +// =========================================================================== +// Image type +// =========================================================================== + +typedef GenericImage ImageType; + +// =========================================================================== +// Main +// =========================================================================== + +// --------------------------------------------------------------------------- +int main(int argc, char **argv) +{ + // Print help or version and exit if requested + HANDLE_HELP_OR_VERSION(); + + // Parse positional arguments + const char *x_fname = NULL; + const char *y_fname = NULL; + const char *z_fname = NULL; + const char *xx_fname = NULL; + const char *xy_fname = NULL; + const char *xz_fname = NULL; + const char *yx_fname = NULL; + const char *yy_fname = NULL; + const char *yz_fname = NULL; + const char *zx_fname = NULL; + const char *zy_fname = NULL; + const char *zz_fname = NULL; + + int posargc = 0; + while (posargc+1 < argc && argv[posargc+1][0] != '-') posargc++; + + switch (posargc) { + case 3: + x_fname = argv[1]; + y_fname = argv[2]; + z_fname = argv[3]; + break; + case 6: + xx_fname = argv[1]; + xy_fname = argv[2]; + yx_fname = argv[3]; + yy_fname = argv[4]; + zx_fname = argv[5]; + zy_fname = argv[6]; + break; + case 9: + xx_fname = argv[1]; + xy_fname = argv[2]; + xz_fname = argv[3]; + yx_fname = argv[4]; + yy_fname = argv[5]; + yz_fname = argv[6]; + zx_fname = argv[7]; + zy_fname = argv[8]; + zz_fname = argv[9]; + break; + default: + PrintHelp(EXECNAME); + exit(1); + } + + // Parse optional arguments + bool usejac = true; + + for (ARGUMENTS_AFTER(posargc)) { + if (OPTION("-jac")) usejac = true; + else if (OPTION("-dif")) usejac = false; + else HANDLE_COMMON_OR_UNKNOWN_OPTION(); + } + + // Read input vector fields + unique_ptr x, y; + + if (x_fname) { + x.reset(new ImageType(x_fname)); + } else { + ImageType xx(xx_fname); + ImageType xy(xy_fname); + unique_ptr xz(xz_fname ? new ImageType(xz_fname) : NULL); + + ImageAttributes attr = xx.Attributes(); + if (attr._t > 1) { + cerr << "Separate input component images must be scalar fields." << endl; + cerr << "Try calling this program with single vector field input and output images instead." << endl; + exit(1); + } + if (xy.Attributes() != attr || (xz.get() && xz->Attributes() != attr)) { + cerr << "Image attributes of input vector component images do not match!" << endl; + exit(1); + } + attr._t = (xz.get() ? 3 : 2); + + x.reset(new ImageType(attr)); + for (int k = 0; k < attr._z; ++k) + for (int j = 0; j < attr._y; ++j) + for (int i = 0; i < attr._x; ++i) { + x->Put(i, j, k, 0, xx.Get(i, j, k)); + x->Put(i, j, k, 1, xy.Get(i, j, k)); + if (xz.get()) x->Put(i, j, k, 2, xz->Get(i, j, k)); + } + } + + if (y_fname) { + y.reset(new ImageType(y_fname)); + } else { + ImageType yx(yx_fname); + ImageType yy(yy_fname); + unique_ptr yz(yz_fname ? new ImageType(yz_fname) : NULL); + + ImageAttributes attr = yx.Attributes(); + if (attr._t > 1) { + cerr << "Separate input component images must be scalar fields." << endl; + cerr << "Try calling this program with single vector field input and output images instead." << endl; + exit(1); + } + if (yy.Attributes() != attr || (yz.get() && yz->Attributes() != attr)) { + cerr << "Image attributes of input vector component images do not match!" << endl; + exit(1); + } + attr._t = (yz.get() ? 3 : 2); + + y.reset(new ImageType(attr)); + for (int k = 0; k < attr._z; ++k) + for (int j = 0; j < attr._y; ++j) + for (int i = 0; i < attr._x; ++i) { + y->Put(i, j, k, 0, yx.Get(i, j, k)); + y->Put(i, j, k, 1, yy.Get(i, j, k)); + if (yz.get()) y->Put(i, j, k, 2, yz->Get(i, j, k)); + } + } + + if (!y->HasSpatialAttributesOf(x.get())) { + cerr << "Spatial attributes of input vector fields do not match!" << endl; + exit(1); + } + if (y->T() != x->T()) { + cerr << "Number of vector components of input vector fields differ!" << endl; + exit(1); + } + + // Compute Lie bracket + ImageType z(x->Attributes()); + liebracket(&z, x.get(), y.get(), usejac); + + // Save output velocities + if (z_fname) { + z.Write(z_fname); + } else { + z.GetFrame(0).Write(zx_fname); + z.GetFrame(1).Write(zy_fname); + if (zz_fname) z.GetFrame(2).Write(zz_fname); + } + + // Clean up + return 0; +} diff --git a/Applications/calculate-logarithmic-map.cc b/Applications/calculate-logarithmic-map.cc new file mode 100644 index 000000000..6426d392b --- /dev/null +++ b/Applications/calculate-logarithmic-map.cc @@ -0,0 +1,215 @@ +/* + * Medical Image Registration ToolKit (MIRTK) + * + * Copyright 2013-2015 Imperial College London + * Copyright 2013-2015 Andreas Schuh + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +#include +#include +#include +#include + +using namespace mirtk; + + +// ============================================================================= +// Global variables with default arguments +// ============================================================================= + +char IntegrationMethod[32] = {"SS"}; +int NumberOfSteps = 64; +int NumberOfBCHSteps = 8; +int NumberOfBCHTerms = 6; +bool SmoothBCHApproximation = false; +bool UseJacobian = false; + +// ============================================================================= +// Help +// ============================================================================= + +// ----------------------------------------------------------------------------- +void PrintHelp(const char *name) +{ + cout << endl; + cout << "Usage: " << name << " [options]" << endl; + cout << " " << name << " [] [] [options]" << endl; + cout << endl; + cout << "Description:" << endl; + cout << " Reads a dense 3D (2D) displacement field, computes the corresponding" << endl; + cout << " stationary 3D (2D) velocity field, and writes the resulting vector field." << endl; + cout << endl; + cout << "Options:" << endl; + cout << " -ss Use scaling and squaring (SS) method for exponentiation. (default)" << endl; + cout << " -euler Use forward Euler method for exponentiation." << endl; + cout << " -steps Number of integration steps (2^n in case of SS). " + << "(default: " << NumberOfSteps << ")" << endl; + cout << " -terms Number of BCH approximation terms (either 2 or 3). (default: " << NumberOfBCHTerms << ")" << endl; + cout << " -iters Number of BCH update steps. (default: " << NumberOfBCHSteps << ")" << endl; + cout << " -smooth Smooth velocities before each update. (default: off)" << endl; + cout << " -jac Use Jacobian for Lie bracket computation. (default: off)" << endl; + PrintCommonOptions(cout); + cout << endl; +} + +// ============================================================================= +// Main +// ============================================================================= + +// ----------------------------------------------------------------------------- +int main(int argc, char *argv[]) +{ + // Print help or version and exit if requested + HANDLE_HELP_OR_VERSION(); + + // Parse positional arguments + const char *disp_fname = NULL; + const char *velo_fname = NULL; + const char *dx_fname = NULL; + const char *dy_fname = NULL; + const char *dz_fname = NULL; + const char *vx_fname = NULL; + const char *vy_fname = NULL; + const char *vz_fname = NULL; + + int posargc = 0; + while (posargc+1 < argc && argv[posargc+1][0] != '-') posargc++; + + switch (posargc) { + case 2: + disp_fname = argv[1]; + velo_fname = argv[2]; + break; + case 4: + dx_fname = argv[1]; + dy_fname = argv[2]; + vx_fname = argv[3]; + vy_fname = argv[4]; + break; + case 6: + dx_fname = argv[1]; + dy_fname = argv[2]; + dz_fname = argv[3]; + vx_fname = argv[4]; + vy_fname = argv[5]; + vz_fname = argv[6]; + break; + default: + PrintHelp(EXECNAME); + exit(1); + } + + // Parse optional arguments + for (ARGUMENTS_AFTER(posargc)) { + if (OPTION("-steps" )) NumberOfSteps = atoi(ARGUMENT); + else if (OPTION("-iters" )) NumberOfBCHSteps = atoi(ARGUMENT); + else if (OPTION("-terms" )) NumberOfBCHTerms = atoi(ARGUMENT); + else if (OPTION("-smooth")) SmoothBCHApproximation = true; + else if (OPTION("-jac")) UseJacobian = true; + else if (OPTION("-int")) { + const char *arg = ARGUMENT; + unsigned int len = strlen(arg); + if (len > 99) len = 99; + for (unsigned int i = 0; i < len; i++) { + if ('a' <= arg[i] && arg[i] <= 'z') IntegrationMethod[i] = 'A' + (arg[i] - 'a'); + else IntegrationMethod[i] = arg[i]; + } + IntegrationMethod[len] = '\0'; + } else if (OPTION("-ss")) { + IntegrationMethod[0] = 'S'; + IntegrationMethod[1] = 'S'; + IntegrationMethod[2] = '\0'; + } else if (OPTION("-euler")) { + IntegrationMethod[0] = 'R'; + IntegrationMethod[1] = 'K'; + IntegrationMethod[2] = 'E'; + IntegrationMethod[3] = '1'; + IntegrationMethod[4] = '\0'; + } else { + HANDLE_COMMON_OR_UNKNOWN_OPTION(); + } + } + + // Read input displacements + unique_ptr d; + + if (disp_fname) { + d.reset(new RealImage(disp_fname)); + } else { + RealImage dx(dx_fname); + RealImage dy(dy_fname); + unique_ptr dz(dz_fname ? new RealImage(dz_fname) : nullptr); + + ImageAttributes attr = dx.Attributes(); + if (attr._t > 1) { + FatalError("Separate input component images must be scalar displacement fields.\n" + "Try calling this program with single vector field input and output images instead."); + } + if (dy.Attributes() != attr || (dz.get() && dz->Attributes() != attr)) { + FatalError("Image attributes of input displacement component images do not match!"); + } + + d.reset(new RealImage(attr, dz.get() ? 3 : 2)); + for (int k = 0; k < attr._z; ++k) + for (int j = 0; j < attr._y; ++j) + for (int i = 0; i < attr._x; ++i) { + d->Put(i, j, k, 0, dx.Get(i, j, k)); + d->Put(i, j, k, 1, dy.Get(i, j, k)); + if (dz.get()) d->Put(i, j, k, 2, dz->Get(i, j, k)); + } + } + + // Instantiate exponential filter + unique_ptr > exp; + if (strcmp(IntegrationMethod, "SS") == 0) { + exp.reset(new VelocityToDisplacementFieldSS()); + } else if (strcmp(IntegrationMethod, "RKE1") == 0 || + strcmp(IntegrationMethod, "FORWARDEULER") == 0 || + strcmp(IntegrationMethod, "EULER") == 0) { + exp.reset(new VelocityToDisplacementFieldEuler()); + } else { + FatalError("Unknown integration method: " << IntegrationMethod); + } + + // Compute stationary velocity field + DisplacementToVelocityFieldBCH dtov; + RealImage v; + + dtov.ExponentialFilter (exp.get()); + dtov.NumberOfSteps (NumberOfSteps); + dtov.NumberOfIterations(NumberOfBCHSteps); + dtov.NumberOfTerms (NumberOfBCHTerms); + dtov.SmoothVelocities (SmoothBCHApproximation); + dtov.UseJacobian (UseJacobian); + + dtov.Input(d.get()); + dtov.Output(&v); + + dtov.Run(); + + // Save output velocities + if (velo_fname) { + v.Write(velo_fname); + } else { + v.GetFrame(0).Write(vx_fname); + v.GetFrame(1).Write(vy_fname); + if (vz_fname) v.GetFrame(2).Write(vz_fname); + } + + return 0; +} diff --git a/Applications/calculate-surface-attributes.cc b/Applications/calculate-surface-attributes.cc new file mode 100644 index 000000000..4976d393d --- /dev/null +++ b/Applications/calculate-surface-attributes.cc @@ -0,0 +1,342 @@ +/* + * Medical Image Registration ToolKit (MIRTK) + * + * Copyright 2013-2015 Imperial College London + * Copyright 2013-2015 Andreas Schuh + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +using namespace mirtk; + + +// ============================================================================= +// Help +// ============================================================================= + +// ----------------------------------------------------------------------------- +void PrintHelp(const char *name) +{ + cout << endl; + cout << "Usage: " << name << " [options]" << endl; + cout << endl; + cout << "Description:" << endl; + cout << "Calculate attributes of input surface such as normals and curvature." << endl; + cout << "If required, as in case of the curvature calculations, the input mesh" << endl; + cout << "is triangulated beforehand if it contains non-triangular faces." << endl; + cout << endl; + cout << "Arguments:" << endl; + cout << " input Input surface mesh." << endl; + cout << " output Output surface mesh." << endl; + cout << endl; + cout << "Optional arguments:" << endl; + cout << " -normals, -point-normals Surface point normals." << endl; + cout << " -cell-normals Surface cell normals." << endl; + cout << " -[no]auto-orient Enable/disable auto-orientation of normals. (default: on)" << endl; + cout << " -[no]splitting Enable/disable splitting of sharp edges. (default: off)" << endl; + cout << " -[no]consistency Enable/disable enforcement of vertex order consistency. (default: on)" << endl; + cout << endl; + cout << "Curvature output options:" << endl; + cout << " -H [] Mean curvature." << endl; + cout << " -K [] Gauss curvature." << endl; + cout << " -C [] Curvedness." << endl; + cout << " -k1 [] Minimum curvature." << endl; + cout << " -k2 [] Maximum curvature." << endl; + cout << " -k1k2 [] [] Principal curvatures." << endl; + cout << " -e1 [] Direction of minimum curvature." << endl; + cout << " -e2 [] Direction of maximum curvature." << endl; + cout << " -normalize Normalize curvature using volume of convex hull." << endl; + cout << " -vtk-curvatures Use vtkCurvatures when possible." << endl; + cout << " -robust-curvatures Do not use vtkCurvatures. Instead, estimate the curvature" << endl; + cout << " tensor field and decompose it to obtain principle curvatures. (default)" << endl; + cout << endl; + cout << " -smooth [] [] []" << endl; + cout << " Smooth calculated scalar curvature measures using a Gaussian smoothing kernel." << endl; + cout << " If sigma2 is specified, an anisotropic kernel with standard deviation" << endl; + cout << " sigma along the direction of minimum curvature, and sigma2 in the" << endl; + cout << " direction of maximum curvature is used. If the value of sigma2 is \"tensor\"" << endl; + cout << " instead of a numeric value, the isotropic Gaussian kernel is oriented" << endl; + cout << " and scaled along each local geometry axis using the curvature tensor." << endl; + PrintCommonOptions(cout); + cout << endl; +} + +// ============================================================================= +// Main +// ============================================================================= + +// ----------------------------------------------------------------------------- +int main(int argc, char *argv[]) +{ + EXPECTS_POSARGS(2); + + const char *input_name = POSARG(1); + const char *output_name = POSARG(2); + + const char *kmin_name = PolyDataCurvature::MINIMUM; + const char *kmax_name = PolyDataCurvature::MAXIMUM; + const char *gauss_name = PolyDataCurvature::GAUSS; + const char *mean_name = PolyDataCurvature::MEAN; + const char *curvedness_name = PolyDataCurvature::CURVEDNESS; + const char *e1_name = PolyDataCurvature::MINIMUM_DIRECTION; + const char *e2_name = PolyDataCurvature::MAXIMUM_DIRECTION; + const char *tensor_name = PolyDataCurvature::TENSOR; + const char *inverse_name = PolyDataCurvature::INVERSE_TENSOR; + + bool point_normals = false; + bool cell_normals = false; + bool auto_orient_normals = true; + bool splitting = false; + bool consistency = true; + bool use_vtkCurvatures = false; + int curvatures = 0; + int tensor_averaging = 3; + bool normalize = false; + int smooth_iterations = 0; + double smooth_sigma = .0; + double smooth_sigma2 = .0; + bool smooth_along_tensor = false; + bool smooth_anisotropic = false; + + for (ALL_OPTIONS) { + if (OPTION("-point-normals") || OPTION("-normals")) point_normals = true; + else if (OPTION("-cell-normals")) cell_normals = true; + else if (OPTION("-auto-orient")) auto_orient_normals = true; + else if (OPTION("-noauto-orient")) auto_orient_normals = false; + else if (OPTION("-consistency")) consistency = true; + else if (OPTION("-noconsistency")) consistency = false; + else if (OPTION("-k1")) { + curvatures |= PolyDataCurvature::Minimum; + if (HAS_ARGUMENT) kmin_name = ARGUMENT; + } + else if (OPTION("-k2")) { + curvatures |= PolyDataCurvature::Maximum; + if (HAS_ARGUMENT) kmax_name = ARGUMENT; + } + else if (OPTION("-k1k2")) { + curvatures |= (PolyDataCurvature::Minimum | PolyDataCurvature::Maximum); + if (HAS_ARGUMENT) { + kmin_name = ARGUMENT; + kmax_name = ARGUMENT; + } + } + else if (OPTION("-H")) { + curvatures |= PolyDataCurvature::Mean; + if (HAS_ARGUMENT) mean_name = ARGUMENT; + } + else if (OPTION("-K")) { + curvatures |= PolyDataCurvature::Gauss; + if (HAS_ARGUMENT) gauss_name = ARGUMENT; + } + else if (OPTION("-C")) { + curvatures |= PolyDataCurvature::Curvedness; + if (HAS_ARGUMENT) curvedness_name = ARGUMENT; + } + else if (OPTION("-n") || OPTION("-normal")) { + curvatures |= PolyDataCurvature::Normal; + } + else if (OPTION("-e1")) { + curvatures |= PolyDataCurvature::MinimumDirection; + if (HAS_ARGUMENT) e1_name = ARGUMENT; + } + else if (OPTION("-e2")) { + curvatures |= PolyDataCurvature::MaximumDirection; + if (HAS_ARGUMENT) e2_name = ARGUMENT; + } + else if (OPTION("-tensor")) { + curvatures |= PolyDataCurvature::Tensor; + if (HAS_ARGUMENT) tensor_name = ARGUMENT; + } + else if (OPTION("-inverse-tensor")) { + curvatures |= PolyDataCurvature::InverseTensor; + if (HAS_ARGUMENT) inverse_name = ARGUMENT; + } + else if (OPTION("-tensor-averaging")) { + PARSE_ARGUMENT(tensor_averaging); + } + else if (OPTION("-normalize")) normalize = true; + else if (OPTION("-vtk-curvatures")) use_vtkCurvatures = true; + else if (OPTION("-robust-curvatures")) use_vtkCurvatures = false; + else if (OPTION("-smooth")) { + smooth_iterations = 1; + if (HAS_ARGUMENT) PARSE_ARGUMENT(smooth_iterations); + if (HAS_ARGUMENT) PARSE_ARGUMENT(smooth_sigma); + if (HAS_ARGUMENT) { + const char *arg = ARGUMENT; + if (strcmp(arg, "tensor") == 0) { + smooth_along_tensor = true; + } else { + smooth_along_tensor = false; + if (!FromString(arg, smooth_sigma2)) { + FatalError("Invalid -smooth argument: " << arg); + } + } + smooth_anisotropic = true; + } + } + else HANDLE_COMMON_OR_UNKNOWN_OPTION(); + } + + if (curvatures == 0 && !point_normals && !cell_normals) { + point_normals = true; + curvatures = PolyDataCurvature::Scalars; + } + + int curvature_type = curvatures; + if (smooth_along_tensor) { + curvature_type |= PolyDataCurvature::Tensor; + } else if (smooth_sigma2 != .0) { + curvature_type |= PolyDataCurvature::MinimumDirection; + curvature_type |= PolyDataCurvature::MaximumDirection; + } + + // Read input surface + vtkSmartPointer input = ReadPolyData(input_name); + + // Triangulate input surface + vtkSmartPointer surface; + if (curvature_type != 0) surface = Triangulate(input); + else surface = input; + surface->BuildLinks(); + + // Calculate normals + if (point_normals || cell_normals) { + if (verbose) cout << "Calculating surface normals...", cout.flush(); + vtkNew filter; + SetVTKInput(filter, surface); + filter->SetSplitting(splitting); + filter->SetComputePointNormals(point_normals); + filter->SetComputeCellNormals(cell_normals); + filter->SetConsistency(consistency); + filter->SetAutoOrientNormals(auto_orient_normals); + filter->FlipNormalsOff(); + filter->Update(); + surface = filter->GetOutput(); + if (verbose) cout << " done" << endl; + } + + // Calculate curvatures + if (curvature_type != 0) { + if (verbose) cout << "Calculating surface curvature measure(s)...", cout.flush(); + + // Compute curvature + PolyDataCurvature curvature; + curvature.Input(surface); + curvature.CurvatureType(curvature_type); + curvature.VtkCurvatures(use_vtkCurvatures); + curvature.TensorAveraging(tensor_averaging); + curvature.Normalize(normalize); + curvature.Run(); + surface = curvature.Output(); + + // Get output arrays + vtkPointData *pd = surface->GetPointData(); + vtkDataArray *kmin = pd->GetArray(PolyDataCurvature::MINIMUM); + vtkDataArray *kmax = pd->GetArray(PolyDataCurvature::MAXIMUM); + vtkDataArray *mean = pd->GetArray(PolyDataCurvature::MEAN); + vtkDataArray *gauss = pd->GetArray(PolyDataCurvature::GAUSS); + vtkDataArray *curvedness = pd->GetArray(PolyDataCurvature::CURVEDNESS); + vtkDataArray *e1 = pd->GetArray(PolyDataCurvature::MINIMUM_DIRECTION); + vtkDataArray *e2 = pd->GetArray(PolyDataCurvature::MAXIMUM_DIRECTION); + vtkDataArray *tensor = pd->GetArray(PolyDataCurvature::TENSOR); + vtkDataArray *inverse = pd->GetArray(PolyDataCurvature::INVERSE_TENSOR); + + // Rename output arrays + if (kmin) kmin ->SetName(kmin_name); + if (kmax) kmax ->SetName(kmax_name); + if (mean) mean ->SetName(mean_name); + if (gauss) gauss ->SetName(gauss_name); + if (curvedness) curvedness->SetName(curvedness_name); + if (e1) e1 ->SetName(e1_name); + if (e2) e2 ->SetName(e2_name); + if (tensor) tensor ->SetName(tensor_name); + if (inverse) inverse ->SetName(inverse_name); + + if (verbose) cout << " done" << endl; + + // Smooth calculated attributes + if (smooth_iterations) { + if (verbose) cout << "Smoothing scalar curvature measures...", cout.flush(); + + PolyDataSmoothing smoother; + smoother.Input(surface); + smoother.SmoothPointsOff(); + if (kmin) smoother.SmoothArray(kmin_name); + if (kmax) smoother.SmoothArray(kmax_name); + if (mean) smoother.SmoothArray(mean_name); + if (gauss) smoother.SmoothArray(gauss_name); + if (curvedness) smoother.SmoothArray(curvedness_name); + if (e1) smoother.SmoothArray(e1_name); + if (e2) smoother.SmoothArray(e2_name); + smoother.NumberOfIterations(smooth_iterations); + smoother.Sigma(-smooth_sigma); // negative: multiple of avg. edge length + if (smooth_anisotropic) { + smoother.Weighting(PolyDataSmoothing::AnisotropicGaussian); + if (smooth_along_tensor) { + smoother.GeometryTensorName(tensor_name); + } else { + smoother.MinimumDirectionName(e2_name); + smoother.MaximumDirectionName(e1_name); + } + smoother.MaximumDirectionSigma(-smooth_sigma2); + } else { + smoother.Weighting(PolyDataSmoothing::Gaussian); + } + smoother.Run(); + vtkPointData *pd = smoother.Output()->GetPointData(); + if (kmin) kmin->DeepCopy(pd->GetArray(kmin_name)); + if (kmax) kmax->DeepCopy(pd->GetArray(kmax_name)); + if (mean) mean->DeepCopy(pd->GetArray(mean_name)); + if (gauss) gauss->DeepCopy(pd->GetArray(gauss_name)); + if (curvedness) curvedness->DeepCopy(pd->GetArray(curvedness_name)); + if (e1) e1->DeepCopy(pd->GetArray(e1_name)); + if (e2) e2->DeepCopy(pd->GetArray(e2_name)); + + if (verbose) cout << " done" << endl; + } + } + + // Remove not requested output arrays which were used for anisotropic smoothing + if ((curvatures & PolyDataCurvature::Tensor) == 0) { + surface->GetPointData()->RemoveArray(tensor_name); + } + if ((curvatures & PolyDataCurvature::MinimumDirection) == 0) { + surface->GetPointData()->RemoveArray(e1_name); + } + if ((curvatures & PolyDataCurvature::MaximumDirection) == 0) { + surface->GetPointData()->RemoveArray(e2_name); + } + + // Write output surface + if (verbose) cout << "Writing output surface to file " << output_name << "...", cout.flush(); + if (!WritePolyData(output_name, surface)) { + if (verbose) cout << " failed" << endl; + FatalError("Failed to write output surface to file " << output_name); + } + if (verbose) cout << " done" << endl; + + return 0; +} diff --git a/Applications/calculate-surface-spectrum.cc b/Applications/calculate-surface-spectrum.cc new file mode 100644 index 000000000..1fb164de5 --- /dev/null +++ b/Applications/calculate-surface-spectrum.cc @@ -0,0 +1,475 @@ +/* + * Medical Image Registration ToolKit (MIRTK) + * + * Copyright 2013-2015 Imperial College London + * Copyright 2013-2015 Andreas Schuh + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +using namespace mirtk; +using namespace mirtk::SpectralDecomposition; + + +// ============================================================================= +// Help +// ============================================================================= + +// ----------------------------------------------------------------------------- +void PrintHelp(const char *name) +{ + cout << endl; + cout << "Usage: " << name << " [options]" << endl; + cout << " " << name << " [options]" << endl; + cout << endl; + cout << "Description:" << endl; + cout << " Performs a spectral analysis of the general graph Laplacian matrix" << endl; + cout << " computed from the given input surface mesh(es)." << endl; + cout << endl; + cout << " If :option:`-target` surface is specified, the eigenvectors are reordered" << endl; + cout << " and their sign flipped to match the eigenvectors of the target surface mesh." << endl; + cout << endl; + cout << " If two input and output surface meshes are specified, a spectral analysis" << endl; + cout << " of the joint graph Laplacian is performed after an initial surface match." << endl; + cout << " The type of initial point correspondences is defined by :option:`-corr`." << endl; + cout << endl; + cout << " The implementation is based on the MATLAB code by Herve Lombaert of [1]_." << endl; + cout << endl; + cout << " .. [1] H. Lombaert, J. Sporring, and K. Siddiqi." << endl; + cout << " Diffeomorphic Spectral Matching of Cortical Surfaces." << endl; + cout << " In the 23rd Image Processing in Medical Imaging (IPMI), 2013." << endl; + cout << endl; + cout << "Arguments:" << endl; + cout << " input File name of input surface (vtkPolyData)." << endl; + cout << " output File name of output surface (vtkPolyData)." << endl; + cout << endl; + cout << "Options:" << endl; + cout << " -target File name of dataset (vtkPolyData) whose spectral components" << endl; + cout << " should be used to adjust the sign and order of the eigenmodes." << endl; + cout << " -corr Type of point correspondences to use for initial correspondences:" << endl; + cout << " \"closest point\" (\"CP\"), \"spectral match\" (\"SM\"). (default: SM)" << endl; + cout << " -corrpar = Name/value pair of correspondence parameter." << endl; + cout << " -corrin Text file naming for each point in input1 the index of the corresponding" << endl; + cout << " point in the second dataset input2, one index per line. (default: :option:`-corr`)" << endl; + cout << " -points Replace output vertex coordinates by first three spectral coordinates." << endl; + PrintCommonOptions(cout); + cout << endl; +} + +// ============================================================================= +// Auxiliaries +// ============================================================================= + +// ----------------------------------------------------------------------------- +/// Write point set with surface eigenmodes +void Write(const char *fname, vtkPointSet *pointset, vtkPolyData *surface, bool compress = true, bool ascii = false, bool as_points = false) +{ + vtkSmartPointer output = pointset; + + vtkDataArray *modes = surface->GetPointData()->GetArray("joint_eigenmodes"); + if (!modes) modes = surface->GetPointData()->GetArray("eigenmodes"); + if (!modes) { + FatalError("Output surface mesh has no eigenmodes point data!"); + } + + if (as_points) { + + vtkSmartPointer points = vtkSmartPointer::New(); + points->SetNumberOfPoints(surface->GetNumberOfPoints()); + + double *p = new double[max(3, int(modes->GetNumberOfComponents()))]; + for (int c = 0; c < 3; ++c) p[c] = .0; + for (vtkIdType i = 0; i < points->GetNumberOfPoints(); ++i) { + modes ->GetTuple(i, p); + points->SetPoint(i, p); + } + delete[] p; + + output = vtkSmartPointer::New(); + output->ShallowCopy(surface); + output->SetPoints(points); + + } else if (pointset != surface) { + + vtkSmartPointer output_modes; + output_modes = vtkSmartPointer::NewInstance(modes); + output_modes->SetName(modes->GetName()); + output_modes->SetNumberOfComponents(modes->GetNumberOfComponents()); + output_modes->SetNumberOfTuples(pointset->GetNumberOfPoints()); + + for (vtkIdType ptId = 0; ptId < pointset->GetNumberOfPoints(); ++ptId) { + for (int j = 0; j < modes->GetNumberOfComponents(); ++j) { + output_modes->SetComponent(ptId, j, .0); + } + } + + vtkIdType origPtId; + vtkDataArray *origPtIds = surface->GetPointData()->GetArray("vtkOriginalPointIds"); + for (vtkIdType ptId = 0; ptId < surface->GetNumberOfPoints(); ++ptId) { + origPtId = static_cast(origPtIds->GetComponent(ptId, 0)); + for (int j = 0; j < modes->GetNumberOfComponents(); ++j) { + output_modes->SetComponent(origPtId, j, modes->GetComponent(ptId, j)); + } + } + + output->GetPointData()->AddArray(output_modes); + } + + WritePointSet(fname, output, compress, ascii); +} + +// ----------------------------------------------------------------------------- +/// Compute eigenmodes of joint intra- and inter-mesh graph Laplacian matrix +int ComputeJointEigenmodes(vtkPolyData *target, const Array &target_sample, + vtkPolyData *source, const Array &source_sample, + const PointCorrespondence *corr, int k, + Matrix &modes, Vector &freq) +{ + // Set intra-mesh affinity weights + const int m = target->GetNumberOfPoints(); + const int n = source->GetNumberOfPoints(); + SparseMatrix::Entries *cols = new SparseMatrix::Entries[m + n]; + AdjacencyMatrix(cols, SparseMatrix::CCS, 0, 0, target); + AdjacencyMatrix(cols, SparseMatrix::CCS, m, m, source); + // Set inter-mesh affinity weights + const FuzzyCorrespondence *fuzzy; + if ((fuzzy = dynamic_cast(corr))) { + ConnectivityMatrix(cols, SparseMatrix::CCS, 0, m, target, &target_sample, + source, &source_sample, &fuzzy->Weight(), false); + ConnectivityMatrix(cols, SparseMatrix::CCS, m, 0, source, &source_sample, + target, &target_sample, &fuzzy->Weight(), true ); + } else { + int i, j; + double w, p1[3], p2[3]; + for (size_t s = 0; s < target_sample.size(); ++s) { + i = target_sample[s]; + j = corr->GetSourceIndex(i); + if (j < 0) { + FatalError("Point correspondence must be fuzzy or map point indices!"); + } + target->GetPoint(i, p1); + source->GetPoint(j, p2); + w = 1.0 / (sqrt(vtkMath::Distance2BetweenPoints(p1, p2)) + EPSILON); + cols[i ].push_back(MakePair(j+m, w)); + cols[j+m].push_back(MakePair(i, w)); + } + for (size_t s = 0; s < source_sample.size(); ++s) { + j = source_sample[s]; + i = corr->GetTargetIndex(j); + if (i < 0) { + FatalError("Point correspondence must be fuzzy or map point indices!"); + } + target->GetPoint(i, p1); + source->GetPoint(j, p2); + w = 1.0 / (sqrt(vtkMath::Distance2BetweenPoints(p1, p2)) + EPSILON); + cols[j+m].push_back(MakePair(i, w)); + cols[i ].push_back(MakePair(j+m, w)); + } + } + // Compute graph Laplacian of joint connectivity graph + SparseMatrix L(SparseMatrix::CCS); + L.Initialize(m + n, m + n, cols); + delete[] cols; + NormalizedLaplacian(L, L); + // Compute eigenmodes of joint graph Laplacian + return ComputeEigenmodes(L, k+1, modes, freq); +} + +// ============================================================================= +// Main +// ============================================================================= + +// ----------------------------------------------------------------------------- +int main(int argc, char *argv[]) +{ + REQUIRES_POSARGS(2); + + const char *input_name [2] = {NULL}; + const char *output_name[2] = {NULL}; + int ndatasets = 0; + + switch (NUM_POSARGS) { + case 2: + input_name [0] = POSARG(1); + output_name[0] = POSARG(2); + ndatasets = 1; + break; + case 4: + input_name [0] = POSARG(1); + input_name [1] = POSARG(2); + output_name[0] = POSARG(3); + output_name[1] = POSARG(4); + ndatasets = 2; + break; + default: + PrintHelp(EXECNAME); + FatalError("Invalid number of positional arguments (" << NUM_POSARGS << ")"); + }; + + // Optional arguments + const char *target_name = NULL; + const char *dofin_name = NULL; + int k = 5; + bool compress = true; + bool ascii = false; + bool as_points = false; + + PointCorrespondence::TypeId ctype = PointCorrespondence::ClosestPoint; + ParameterList cparam; + Array cfeature_name; + Array cfeature_weight; + + for (ALL_OPTIONS) { + if (OPTION("-target")) target_name = ARGUMENT; + else if (OPTION("-dofin")) dofin_name = ARGUMENT; + else if (OPTION("-k") || OPTION("-dim")) k = atoi(ARGUMENT); + else if (OPTION("-corr")) PARSE_ARGUMENT(ctype); + else if (OPTION("-corrin")) { + ctype = PointCorrespondence::FiducialMatch; + Insert(cparam, "Correspondence map", ARGUMENT); + } + else if (OPTION("-corrpar")) { + Insert(cparam, ARGUMENT, ARGUMENT); + } + else if (OPTION("-feature")) { + cfeature_name.push_back(ARGUMENT); + if (HAS_ARGUMENT) { + double w; + PARSE_ARGUMENT(w); + cfeature_weight.push_back(w); + } + else cfeature_weight.push_back(1.0); + } + else if (OPTION("-points")) as_points = true; + else if (OPTION("-compress")) compress = true; + else if (OPTION("-nocompress") ) compress = false; + else if (OPTION("-ascii") || OPTION("-nobinary")) ascii = true; + else if (OPTION("-binary") || OPTION("-noascii") ) ascii = false; + else HANDLE_COMMON_OR_UNKNOWN_OPTION(); + } + + // Default correspondence features and normalization of names to lowercase + if (cfeature_name.empty()) { + cfeature_name .push_back("eigenmodes"); + cfeature_weight.push_back(1.0); + } else { + for (size_t i = 0; i < cfeature_name.size(); ++i) { + transform(cfeature_name[i].begin(), cfeature_name[i].end(), + cfeature_name[i].begin(), ::tolower); + } + } + + // Ensure number of spectral components is positive + if (k <= 0) { + FatalError("Argument for option -k/-dim must be positive!"); + } + if (as_points && k < 3) { + FatalError("Argument for option -k/-dim must be >= 3 when -points should be written!"); + } + + // Read input datasets + vtkSmartPointer dataset[2]; + vtkSmartPointer surface[2]; + for (int i = 0; i < ndatasets; ++i) { + if (verbose > 1) cout << "Reading point set " << (i+1) << " from " << input_name[i] << endl; + dataset[i] = ReadPointSet(input_name[i]); + if (dataset[i] == NULL || dataset[i]->GetNumberOfPoints() == 0) { + FatalError("Failed to read dataset " << (i+1) << " or dataset is empty!"); + } + surface[i] = DataSetSurface(dataset[i], true); + } + + // Transform first surface by specified transformation + if (dofin_name) { + unique_ptr dof(Transformation::New(dofin_name)); + RegisteredSurface transformed; + transformed.InputSurface(surface[0]); + transformed.Transformation(dof.get()); + transformed.Initialize(); + transformed.Update(); + surface[0] = transformed.Surface(); + } + + // Skip individual spectral analysis if joint analysis requested based on + // pre-defined correspondences, non-spectral matches, or known eigenmodes + bool individual_analysis[2] = {false}; + bool spectral_match = false; + if (ndatasets == 1) { + spectral_match = !as_points; // -points option can be applied to single input dataset + } else if (ndatasets > 1) { + switch (ctype) { + case PointCorrespondence::FiducialMatch: + spectral_match = false; + break; + case PointCorrespondence::SpectralMatch: + spectral_match = true; + break; + default: + spectral_match = false; + for (size_t i = 0; i < cfeature_name.size(); ++i) { + if (cfeature_name[i].find("spectral" ) != string::npos || + cfeature_name[i].find("eigenmode") != string::npos) { + spectral_match = true; + break; + } + } + } + } + if (spectral_match) { + for (int i = 0; i < ndatasets; ++i) { + individual_analysis[i] = (!surface[i]->GetPointData()->HasArray("eigenmodes") && + !surface[i]->GetPointData()->HasArray("joint_eigenmodes")); + } + } + + // Individual spectral analysis of input datasets + Matrix modes[2]; + Vector freq [2]; + for (int i = 0; i < ndatasets; ++i) { + if (individual_analysis[i]) { + if (ComputeEigenmodes(surface[i], k, modes[i], freq[i]) < k) { + FatalError("Failed to find " << k << " eigenmodes of dataset " << (i+1)); + } + if (verbose) { + cout << "Frequencies"; + if (ndatasets > 1) cout << " of dataset " << (i+1); + cout << ": "; + for (int c = 0; c < k; ++c) { + if (c > 0) cout << " "; + cout << scientific << freq[i](c); + } + cout << endl; + } + } + } + + // Joint spectral analysis + if (ndatasets == 2) { + + // Adjust sign and order of individually computed eigenmodes + if (spectral_match) { + Vector cost = MatchEigenmodes(surface[0]->GetPoints(), modes[0], freq[0], + surface[1]->GetPoints(), modes[1], freq[1]); + + // TODO: Transform modes[1] using coherent point drift algorithm + + // Add individual eigenmodes to output point data + for (int i = 0; i < ndatasets; ++i) { + SetEigenmodes(surface[i], modes[i], "eigenmodes"); + modes[i].Clear(); + freq [i].Clear(); + } + } + + // Obtain correspondence map for initial links + RegisteredSurface target, source; + target.InputSurface(surface[0]); + source.InputSurface(surface[1]); + target.Initialize(); + source.Initialize(); + target.Update(); + source.Update(); + + unique_ptr cmap(PointCorrespondence::New(ctype)); + cmap->FromTargetToSource(true); + cmap->FromSourceToTarget(true); + cmap->Parameter(cparam); + cmap->Target(&target); + cmap->Source(&source); + for (size_t i = 0; i < cfeature_name.size(); ++i) { + cmap->AddFeature(cfeature_name[i].c_str(), cfeature_weight[i]); + } + cmap->Initialize(); + cmap->Update(); + + // Draw samples for which to add connectivity links + Array target_sample; + Array source_sample; + using PointCorrespondenceUtils::SamplePoints; + SamplePoints(&target, target_sample, target.NumberOfPoints() / 10); + SamplePoints(&source, source_sample, source.NumberOfPoints() / 10); + + // Perform spectral analysis of joint graph Laplacian + Matrix modes; + Vector freq; + if (ComputeJointEigenmodes(target.Surface(), target_sample, + source.Surface(), source_sample, + cmap.get(), k, modes, freq) < k) { + FatalError("Error: Failed to find " << k << " eigenmodes of joint graph Laplacian!"); + } + + // Weight eigenmodes + Vector w(k); + double wsum = .0; + for (int i = 0; i < k; ++i) { + wsum += 1.0 / sqrt(freq(i)); + } + for (int i = 0; i < k; ++i) { + modes.ScaleCol(i, (1.0 / sqrt(freq(i))) / wsum); + } + + // Add eigenmodes to output point data + SetEigenmodes(surface[0], modes, 0, k, "joint_eigenmodes"); + SetEigenmodes(surface[1], modes, target.NumberOfPoints(), k, "joint_eigenmodes"); + + } else { + + // Add individual eigenmodes to output point data + for (int i = 0; i < ndatasets; ++i) { + if (individual_analysis[i]) { + SetEigenmodes(surface[i], modes[i], "eigenmodes"); + modes[i].Clear(); + freq [i].Clear(); + } + } + + } + + // Match spectral components of input datasets with those of target dataset + if (target_name) { + if (verbose > 1) cout << "Reading target from " << target_name << endl; + vtkSmartPointer target = ReadPointSet(target_name); + if (target == NULL || target->GetNumberOfPoints() == 0) { + Warning("Failed to read target dataset or dataset has no points!" + " Skipping matching of spectral components therefore."); + } else { + for (int i = 0; i < ndatasets; ++i) { + // TODO: Match eigenmodes + } + } + } + + // Write datasets together with their spectral components + for (int i = 0; i < ndatasets; ++i) { + Write(output_name[i], dataset[i], surface[i], compress, ascii, as_points); + if (verbose > 1) cout << "Wrote result dataset to " << output_name[i] << endl; + } + + return 0; +} diff --git a/Applications/calculate.cc b/Applications/calculate.cc new file mode 100644 index 000000000..4279c3e6d --- /dev/null +++ b/Applications/calculate.cc @@ -0,0 +1,576 @@ +/* + * Medical Image Registration ToolKit (MIRTK) + * + * Copyright 2013-2015 Imperial College London + * Copyright 2013-2015 Andreas Schuh + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +#include + +#include +#include +#include + +#if MIRTK_Image_WITH_VTK +# include +# include +#endif + +using namespace mirtk; +using namespace mirtk::data; +using namespace mirtk::data::op; +using namespace mirtk::data::statistic; + + +// ============================================================================= +// Help +// ============================================================================= + +// ----------------------------------------------------------------------------- +void PrintHelp(const char *name) +{ + cout << endl; + cout << "Usage: " << name << " [options]" << endl; + cout << endl; + cout << "Description:" << endl; + cout << " This tool can be used for basic calculations from a sequence of data values read" << endl; + cout << " either from an image or a VTK pointset. It can be used, for example, to add two" << endl; + cout << " data sequences and to divide the result by a constant. The current sequence can" << endl; + cout << " be written to an output file again using the -out option. Additionally, statistics" << endl; + cout << " of the current data sequence can be computed such as the mean or variance." << endl; + cout << " The order of the data transformations and calculation of statistics is determined" << endl; + cout << " by the order of the command-line arguments." << endl; + cout << endl; + cout << " By default, data statistics are printed to STDOUT in a human readable format." << endl; + cout << " This output can be appended to a text file using the -append option instead." << endl; + cout << " For a more machine readable output, e.g., as comma separated values (CSV)," << endl; + cout << " specify a delimiting string using the -d option. In this case, a header line" << endl; + cout << " is also printed when the -header option is given with optional user specified" << endl; + cout << " column names for the individual output values." << endl; + cout << endl; + cout << "Input options:" << endl; + cout << " -scalars Name of input point data array. (default: active SCALARS array)" << endl; + cout << endl; + cout << "Data manipulation options:" << endl; + cout << " -mask Exclude values equal a given threshold or with zero input mask value." << endl; + cout << " Note that this does not modify the data values, but only marks them to be" << endl; + cout << " ignored from now on. Use -outside/-pad following this operation to set an outside value." << endl; + cout << " -mask-inside Exclude values which are inside a given closed interval." << endl; + cout << " -mask-outside Exclude values which are outside a given open interval." << endl; + cout << " -mask-below Exclude values less or equal a given threshold." << endl; + cout << " -mask-above Exclude values greator or equal a given threshold." << endl; + cout << " -even Exclude values which are not an even number when cast to an integer." << endl; + cout << " -odd Exclude values which are not an odd number when cast to an integer." << endl; + cout << " -threshold [] Clamp values which are less or equal and optionally greater or equal a threshold." << endl; + cout << " -threshold-below Clamp values less or equal a given threshold." << endl; + cout << " -threshold-above Clamp values greator or equal a given threshold." << endl; + cout << " -rescale Linearly rescale values to the interval [min, max]." << endl; + cout << " -inside Set new value for all currently unmasked data values." << endl; + cout << " -outside Set new value for all currently masked data values." << endl; + cout << " -out Write current data sequence to file in the format of the input file." << endl; + cout << " -add Add constant value or data sequence read from specified file." << endl; + cout << " -sub Subtract constant value or data sequence read from specified file." << endl; + cout << " -mul Multiply by constant value or data sequence read from specified file." << endl; + cout << " -div Divide by constant value or data sequence read from specified file." << endl; + cout << " -abs Take absolute values." << endl; + cout << " -pow Raise values to the power of the given exponent." << endl; + cout << " -sq -square Raise values to the power of 2 (i.e, -pow 2)." << endl; + cout << " -sqrt Calculate square root of each value (i.e, -pow .5)." << endl; + cout << " -exp Calculate exponential of data sequence." << endl; + cout << " -log [] [] Compute logarithm after applying an optional threshold." << endl; + cout << " (default threshold: min double, default base: e)" << endl; + cout << " -lb [] Compute binary logarithm, alias for -log 2." << endl; + cout << " -log2 [] Compute binary logarithm, alias for -log 2." << endl; + cout << " -ln [] Compute natural logarithm, alias for -log []." << endl; + cout << " -lg [] Compute logarithm to base 10, alias for -log 10." << endl; + cout << " -log10 [] Compute logarithm to base 10, alias for -log 10." << endl; + cout << endl; + cout << "Data statistics options:" << endl; + cout << " -append Append output to a file. (default: STDOUT)" << endl; + cout << " -delimiter | -delim | -d | -sep Delimiting character(s). (default: '')" << endl; + cout << " -header [...] Request output of header line if delimiter was" << endl; + cout << " specified as well. If the output is appended to" << endl; + cout << " a text file, the header is only printed if it does not exist." << endl; + cout << " If no or fewer custom column names are given, the default" << endl; + cout << " names for each statistic are printed. (default: none)" << endl; + cout << " -prefix ... One or more prefix strings to print. If no delimiter" << endl; + cout << " is specified, the concatenated strings are printed before" << endl; + cout << " each line of the output. Otherwise, each prefix string" << endl; + cout << " is printed as entry for the first columns in the delimited" << endl; + cout << " output row, separated by the specified delimiter. (default: none)" << endl; + cout << " -precision Number of significant digits. (default: 5)" << endl; + cout << " -mean | -avg Print mean intensity. (default: on)" << endl; + cout << " -variance | -var Print variance of intensity values. (default: off)" << endl; + cout << " -sigma | -std Print standard deviation of intensity values. (default: on)" << endl; + cout << " -minimum | -min Print minimum intensity value. (default: off)" << endl; + cout << " -maximum | -max Print maximum intensity value. (default: off)" << endl; + cout << " -extrema | -minmax Print minimum and maximum intensity value. (default: on)" << endl; + cout << " -range Print range of intensity values (i.e., max - min). (default: off)" << endl; + cout << " -percentile | -pct Print n-th percentile. (default: none)" << endl; + PrintCommonOptions(cout); + cout << endl; + cout << "Examples:" << endl; + cout << endl; + cout << " " << name << " mni305.nii.gz" << endl; + cout << " Mean = 26.9753" << endl; + cout << " Standard deviation = 50.3525" << endl; + cout << " Extrema = [0, 254]" << endl; + cout << " Range = 254" << endl; + cout << endl; + cout << " " << name << " mni305.nii.gz -pct 77" << endl; + cout << " 77th percentile = 25" << endl; + cout << endl; + cout << " " << name << " mni305.nii.gz -padding 25 -range -percentile 25 50 75 -prefix MNI305 '[>25]'" << endl; + cout << " MNI305 [>25] range = 254" << endl; + cout << " MNI305 [>25] 25th percentile = 69" << endl; + cout << " MNI305 [>25] 50th percentile = 113" << endl; + cout << " MNI305 [>25] 75th percentile = 150" << endl; + cout << endl; + cout << " " << name << " mni305.nii.gz -d , -prefix MNI305" << endl; + cout << " MNI305,26.9753,50.3525,0,254,254 [no newline at end of line]" << endl; + cout << endl; + cout << " " << name << " mni305.nii.gz -d , -prefix MNI305 -header" << endl; + cout << " ,Mean,Sigma,Min,Max,Range" << endl; + cout << " MNI305,26.9753,50.3525,0,254,254" << endl; + cout << endl; + cout << " " << name << " mni305.nii.gz -d , -prefix MNI305 -header ID Mean SD" << endl; + cout << " ID,Mean,SD,Min,Max,Range" << endl; + cout << " MNI305,26.9753,50.3525,0,254,254" << endl; + cout << endl; + cout << " " << name << " a.nii.gz + b.nii.gz = c.nii.gz" << endl; + cout << endl; + cout << " " << name << " a.vtk + b.nii.gz - 10 / c.nii = d.vtk" << endl; + cout << " Adds data values at identical sequential memory indices in a and b," << endl; + cout << " subtracts the constant 10, and then divides by the values in image c." << endl; + cout << endl; + cout << " Note: Operations are always executed from left to right," << endl; + cout << " i.e., no mathematical operator precedence is considered!" << endl; + cout << endl; +} + +// ============================================================================= +// Main +// ============================================================================= + +// ----------------------------------------------------------------------------- +// Some special options do not start with a '-' as otherwise required +#undef HAS_ARGUMENT +#define HAS_ARGUMENT \ + _IsArgument(ARGIDX, argc, argv) && \ + strcmp(argv[ARGIDX+1], "+") != 0 && \ + strcmp(argv[ARGIDX+1], "/") != 0 && \ + strcmp(argv[ARGIDX+1], "=") != 0 + +// ----------------------------------------------------------------------------- +int main(int argc, char **argv) +{ + // Initial data values + REQUIRES_POSARGS(1); + + double *data = NULL; + int datatype = MIRTK_VOXEL_DOUBLE; + ImageAttributes attr; + +#if MIRTK_Image_WITH_VTK + const char *scalars_name = NULL; + for (ARGUMENTS_AFTER(1)) { + if (OPTION("-scalars")) scalars_name = ARGUMENT; + } + vtkSmartPointer dataset; + int n = Read(POSARG(1), data, &datatype, &attr, &dataset, scalars_name); +#else // MIRTK_Image_WITH_VTK + int n = Read(POSARG(1), data, &datatype, &attr); +#endif // MIRTK_Image_WITH_VTK + + // Optional arguments + const char *append_name = NULL; + const char *delimiter = NULL; + bool print_header = false; + int digits = 5; + + Array header; + Array prefix; + + Array > ops; + + for (ARGUMENTS_AFTER(1)) { + if (OPTION("-scalars")) scalars_name = ARGUMENT; + else if (OPTION("-append")) append_name = ARGUMENT; + else if (OPTION("-prefix")) { + do { + prefix.push_back(ARGUMENT); + } while (HAS_ARGUMENT); + } else if (OPTION("-header")) { + print_header = true; + while (HAS_ARGUMENT) header.push_back(ARGUMENT); + } else if (OPTION("-mask")) { + double c; + do { + const char *arg = ARGUMENT; + if (FromString(arg, c)) ops.push_back(unique_ptr(new Mask(c))); + else ops.push_back(unique_ptr(new Mask(arg))); + } while (HAS_ARGUMENT); + } else if (OPTION("-mask-outside")) { + const char *arg = ARGUMENT; + double lower, upper; + if (!FromString(arg, lower)) { + cerr << "Invalid -mask-outside lower argument, must be a number!" << endl; + } + arg = ARGUMENT; + if (!FromString(arg, upper)) { + cerr << "Invalid -mask-outside upper argument, must be a number!" << endl; + } + ops.push_back(unique_ptr(new MaskOutsideOpenInterval(lower, upper))); + } else if (OPTION("-mask-inside")) { + const char *arg = ARGUMENT; + double lower, upper; + if (!FromString(arg, lower)) { + cerr << "Invalid -mask-inside lower argument, must be a number!" << endl; + } + arg = ARGUMENT; + if (!FromString(arg, upper)) { + cerr << "Invalid -mask-inside upper argument, must be a number!" << endl; + } + ops.push_back(unique_ptr(new MaskInsideInterval(lower, upper))); + } else if (OPTION("-mask-below")) { + const char *arg = ARGUMENT; + double lower; + if (!FromString(arg, lower)) { + cerr << "Invalid -mask-below argument, must be a number!" << endl; + } + ops.push_back(unique_ptr(new MaskOutsideOpenInterval(lower, numeric_limits::infinity()))); + } else if (OPTION("-mask-above")) { + const char *arg = ARGUMENT; + double upper; + if (!FromString(arg, upper)) { + cerr << "Invalid -mask-above argument, must be a number!" << endl; + } + ops.push_back(unique_ptr(new MaskOutsideOpenInterval(-numeric_limits::infinity(), upper))); + } else if (OPTION("-even")) { + ops.push_back(unique_ptr(new MaskOddValues())); + } else if (OPTION("-odd")) { + ops.push_back(unique_ptr(new MaskEvenValues())); + } else if (OPTION("-threshold")) { + double lower, upper; + if (!FromString(ARGUMENT, lower)) { + cerr << "Invalid -threshold argument, must be a number!" << endl; + exit(1); + } + if (HAS_ARGUMENT) { + if (!FromString(ARGUMENT, upper)) { + cerr << "Invalid -threshold upper argument, must be a number!" << endl; + exit(1); + } + ops.push_back(unique_ptr(new Clamp(lower, upper))); + } else { + ops.push_back(unique_ptr(new LowerThreshold(lower))); + } + } else if (OPTION("-threshold-below") || OPTION("-lower-threshold")) { + double threshold; + if (!FromString(ARGUMENT, threshold)) { + cerr << "Invalid -threshold-below, must be a number!" << endl; + exit(1); + } + ops.push_back(unique_ptr(new LowerThreshold(threshold))); + } else if (OPTION("-threshold-above") || OPTION("-upper-threshold")) { + double threshold; + if (!FromString(ARGUMENT, threshold)) { + cerr << "Invalid -threshold-above, must be a number!" << endl; + exit(1); + } + ops.push_back(unique_ptr(new UpperThreshold(threshold))); + } else if (OPTION("-rescale")) { + double min, max; + if (!FromString(ARGUMENT, min)) { + cerr << "Invalid -rescale minimum, must be a number!" << endl; + exit(1); + } + if (!FromString(ARGUMENT, max)) { + cerr << "Invalid -rescale maximum, must be a number!" << endl; + exit(1); + } + ops.push_back(unique_ptr(new Rescale(min, max))); + } else if (OPTION("-inside")) { + double inside_value; + if (!FromString(ARGUMENT, inside_value)) { + cerr << "Invalid -inside value, must be a number!" << endl; + exit(1); + } + ops.push_back(unique_ptr(new SetInsideValue(inside_value))); + } else if (OPTION("-outside") || OPTION("-pad")) { + double outside_value; + if (!FromString(ARGUMENT, outside_value)) { + cerr << "Invalid -outside value, must be a number!" << endl; + exit(1); + } + ops.push_back(unique_ptr(new SetOutsideValue(outside_value))); + // Data transformations + } else if (OPTION("-add") || OPTION("-plus") || OPTION("+")) { + const char *arg = ARGUMENT; + double c; + if (FromString(arg, c)) { + ops.push_back(unique_ptr(new Add(c))); + } else { + ops.push_back(unique_ptr(new Add(arg))); + } + } else if (OPTION("-sub") || OPTION("-subtract") || OPTION("-minus") || OPTION("-")) { + const char *arg = ARGUMENT; + double c; + if (FromString(arg, c)) { + ops.push_back(unique_ptr(new Sub(c))); + } else { + ops.push_back(unique_ptr(new Sub(arg))); + } + } else if (OPTION("-mul") || OPTION("-multiply-by") || OPTION("-times") || OPTION("*")) { + const char *arg = ARGUMENT; + double c; + if (FromString(arg, c)) { + ops.push_back(unique_ptr(new Mul(c))); + } else { + ops.push_back(unique_ptr(new Mul(arg))); + } + } else if (OPTION("-div") || OPTION("-divide-by") || OPTION("-over") || OPTION("/")) { + const char *arg = ARGUMENT; + double c; + if (FromString(arg, c)) { + if (fequal(c, .0)) { + cerr << "Invalid -div argument, value must not be zero!" << endl; + exit(1); + } + ops.push_back(unique_ptr(new Div(c))); + } else { + ops.push_back(unique_ptr(new Div(arg))); + } + } else if (OPTION("-abs")) { + ops.push_back(unique_ptr(new Abs())); + } else if (OPTION("-pow") || OPTION("-power")) { + const char *arg = ARGUMENT; + double exponent; + if (!FromString(arg, exponent)) { + cerr << "Invalid -power value, must be a number!" << endl; + exit(1); + } + ops.push_back(unique_ptr(new Pow(exponent))); + } else if (OPTION("-sqrt")) { + ops.push_back(unique_ptr(new Pow(.5))); + } else if (OPTION("-square") || OPTION("-sq")) { + ops.push_back(unique_ptr(new Pow(2.0))); + } else if (OPTION("-exp")) { + ops.push_back(unique_ptr(new Exp())); + } else if (OPTION("-log") || OPTION("-log2") || OPTION("-log10") || OPTION("-lb") || OPTION("-ln") || OPTION("-lg")) { + double threshold = numeric_limits::min(); + if (HAS_ARGUMENT) { + if (!FromString(ARGUMENT, threshold) || threshold <= .0) { + cerr << "Invalid -log threshold argument, must be a positive number" << endl; + exit(1); + } + } + Op *op = NULL; + if (strcmp(OPTNAME, "-log") == 0) { + if (HAS_ARGUMENT) { + double base; + if (!FromString(ARGUMENT, base)) { + cerr << "Invalid -log base argument, must be a positive number" << endl; + exit(1); + } + op = new Log(base, threshold); + } else { + op = new Ln(threshold); + } + } else if (strcmp(OPTNAME, "-log2") == 0 || strcmp(OPTNAME, "-lb") == 0) { + op = new Lb(threshold); + } else if (strcmp(OPTNAME, "-log10") == 0 || strcmp(OPTNAME, "-lg") == 0) { + op = new Lg(threshold); + } else if (strcmp(OPTNAME, "-ln") == 0) { + op = new Ln(threshold); + } + ops.push_back(unique_ptr(op)); + } else if (OPTION("-o") || OPTION("-out") || OPTION("-output") || OPTION("=")) { + const char *fname = ARGUMENT; + int dtype = datatype; + if (HAS_ARGUMENT) { + const char *arg = ARGUMENT; + dtype = ToDataType(arg); + if (dtype == MIRTK_VOXEL_UNKNOWN) { + cerr << "Invalid -out data type " << arg << endl; + exit(1); + } + } +#if MIRTK_Image_WITH_VTK + ops.push_back(unique_ptr(new Write(fname, dtype, attr, dataset, scalars_name))); +#else + ops.push_back(unique_ptr(new Write(fname, dtype, attr))); +#endif + // Data statistics + } else if (OPTION("-mean") || OPTION("-average") || OPTION("-avg")) { + ops.push_back(unique_ptr(new Mean())); + } else if (OPTION("-sigma") || OPTION("-stddev") || OPTION("-stdev") || OPTION("-std") || OPTION("-sd")) { + ops.push_back(unique_ptr(new StDev())); + } else if (OPTION("-normal-distribution") || + OPTION("-mean+sigma") || OPTION("-mean+stddev") || OPTION("-mean+stdev") || OPTION("-mean+std") || OPTION("-mean+sd") || + OPTION("-avg+sigma") || OPTION("-avg+stddev") || OPTION("-avg+stdev") || OPTION("-avg+std") || OPTION("-avg+sd")) { + ops.push_back(unique_ptr(new NormalDistribution())); + } else if (OPTION("-variance") || OPTION("-var")) { + ops.push_back(unique_ptr(new Var())); + } else if (OPTION("-minimum") || OPTION("-min")) { + ops.push_back(unique_ptr(new Min())); + } else if (OPTION("-maximum") || OPTION("-max")) { + ops.push_back(unique_ptr(new Max())); + } else if (OPTION("-extrema") || OPTION("-minmax")) { + ops.push_back(unique_ptr(new Extrema())); + } else if (OPTION("-range")) { + ops.push_back(unique_ptr(new Range())); + } else if (OPTION("-percentile") || OPTION("-pct") || OPTION("-p")) { + do { + int p; + if (FromString(ARGUMENT, p) && 0 <= p && p <= 100) { + ops.push_back(unique_ptr(new Percentile(p))); + } else { + cerr << "Invalid -percentile value, must be integer in the range [0, 100]!" << endl; + exit(1); + } + } while (HAS_ARGUMENT); + } else if (OPTION("-lower-percentile-mean") || OPTION("-lpctavg")) { + do { + int p; + if (FromString(ARGUMENT, p) && 0 <= p && p <= 100) { + ops.push_back(unique_ptr(new LowerPercentileMean(p))); + } else { + cerr << "Invalid -lower-percentile-mean value, must be integer in the range [0, 100]!" << endl; + exit(1); + } + } while (HAS_ARGUMENT); + } else if (OPTION("-upper-percentile-mean") || OPTION("-upctavg")) { + do { + int p; + if (FromString(ARGUMENT, p) && 0 <= p && p <= 100) { + ops.push_back(unique_ptr(new UpperPercentileMean(p))); + } else { + cerr << "Invalid -upper-percentile-mean value, must be integer in the range [0, 100]!" << endl; + exit(1); + } + } while (HAS_ARGUMENT); + } else if (OPTION("-delimiter") || OPTION("-delim") || OPTION("-d") || OPTION("-sep")) { + delimiter = ARGUMENT; + } else if (OPTION("-precision") || OPTION("-digits")) { + if (!FromString(ARGUMENT, digits) || digits < 0) { + cerr << "Invalid -precision argument, value must be non-negative integer!" << endl; + exit(1); + } + } else { + HANDLE_COMMON_OR_UNKNOWN_OPTION(); + } + } + + // If delimiter explicitly set to empty string, use none + if (delimiter && delimiter[0] == '\0') delimiter = NULL; + + // Default statistics to compute + if (ops.empty()) { + ops.push_back(unique_ptr(new Mean())); + ops.push_back(unique_ptr(new StDev())); + ops.push_back(unique_ptr(new Extrema())); + ops.push_back(unique_ptr(new Range())); + } + + // Initial data mask + bool *mask = new bool[n]; + for (int i = 0; i < n; ++i) { + if (IsNaN(data[i])) { + mask[i] = false; + } else { + mask[i] = true; + } + } + + // Open output file to append to or use STDOUT if none specified + ofstream ofs; + if (append_name) { + if (print_header) { + ifstream ifs(append_name); + if (ifs.is_open()) { + print_header = false; + ifs.close(); + } + } + ofs.open(append_name, ios_base::app); + if (!ofs.is_open()) { + cerr << "Cannot append to file " << append_name << endl; + delete[] mask; + exit(1); + } + } + ostream &out = (ofs.is_open() ? ofs : cout); + + // Print column names if requested + if (delimiter && print_header) { + size_t c = 0; + for (size_t i = 0; i < prefix.size(); ++i, ++c) { + if (c > 0) out << delimiter; + if (c < header.size()) out << header[c]; + } + for (size_t i = 0; i < ops.size(); ++i) { + Statistic *stat = dynamic_cast(ops[i].get()); + if (stat) { + for (size_t j = 0; j < stat->Names().size(); ++j, ++c) { + if (c > 0) out << delimiter; + if (c < header.size()) out << header[c]; + else out << stat->Names()[j]; + } + } + } + out << endl; + } + + // Process input data, either transform it or compute statistics from it + for (size_t i = 0; i < ops.size(); ++i) ops[i]->Process(n, data, mask); + delete[] mask; + + // Print image statistics + if (delimiter) { + for (size_t i = 0; i < prefix.size(); ++i) { + if (i > 0) out << delimiter; + out << prefix[i]; + } + bool first = prefix.empty(); + for (size_t i = 0; i < ops.size(); ++i) { + Statistic *stat = dynamic_cast(ops[i].get()); + if (stat) { + if (!first) out << delimiter; + else first = false; + stat->PrintValues(out, digits, delimiter); + } + } + // No newline at end of row if printing results to STDOUT which in this + // case is usually assigned to a string in a calling script + if (print_header || ofs.is_open()) out << endl; + } else { + string prefix_string; + for (size_t i = 0; i < prefix.size(); ++i) { + if (i > 0) prefix_string += ' '; + prefix_string += prefix[i]; + } + for (size_t i = 0; i < ops.size(); ++i) { + Statistic *stat = dynamic_cast(ops[i].get()); + if (stat) stat->Print(out, digits, prefix_string.c_str()); + } + } + + ofs.close(); + return 0; +} diff --git a/Applications/close-image.cc b/Applications/close-image.cc new file mode 100644 index 000000000..c3011e2d0 --- /dev/null +++ b/Applications/close-image.cc @@ -0,0 +1,108 @@ +/* + * Medical Image Registration ToolKit (MIRTK) + * + * Copyright 2008-2015 Imperial College London + * Copyright 2008-2013 Daniel Rueckert, Julia Schnabel + * Copyright 2013-2015 Andreas Schuh + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +#include +#include + +using namespace mirtk; + + +// ============================================================================= +// Help +// ============================================================================= + +// ----------------------------------------------------------------------------- +void PrintHelp(const char *name) +{ + cout << "\n"; + cout << "Usage: " << name << " [options]\n"; + cout << "\n"; + cout << "Description:\n"; + cout << " Closes holes by first dilating and then eroding an input image.\n"; + cout << "\n"; + cout << "Arguments:\n"; + cout << " input Input intensity/segmentation image.\n"; + cout << " output Closed output image.\n"; + cout << "\n"; + cout << "Optional arguments:\n"; + cout << " -iterations Number of dilation/erosion iterations. (default: 1)\n"; + cout << " -connectivity Type of voxel connectivity (4, 6, 18, or 26). (default: 26)\n"; + PrintStandardOptions(cout); + cout << "\n"; + cout.flush(); +} + +// ============================================================================= +// Auxiliaries +// ============================================================================= + +// ----------------------------------------------------------------------------- +template +void Close(BaseImage *image, int iterations, ConnectivityType connectivity) +{ + Dilate(image, iterations, connectivity); + Erode (image, iterations, connectivity); +} + +// ============================================================================= +// Main +// ============================================================================= + +// ----------------------------------------------------------------------------- +int main(int argc, char *argv[]) +{ + EXPECTS_POSARGS(2); + + const char *input_name = POSARG(1); + const char *output_name = POSARG(2); + + int iterations = 1; + ConnectivityType connectivity = CONNECTIVITY_26; + + for (ALL_OPTIONS) { + if (OPTION("-iterations") || OPTION("-iter")) { + PARSE_ARGUMENT(iterations); + } + else if (OPTION("-connectivity") || OPTION("-neighbors") || OPTION("-number-of-neighbors")) { + PARSE_ARGUMENT(connectivity); + } + else HANDLE_STANDARD_OR_UNKNOWN_OPTION(); + } + + unique_ptr image(BaseImage::New(input_name)); + + if (verbose) cout << "Closing ... ", cout.flush(); + switch (image->GetDataType()) { + case MIRTK_VOXEL_BINARY: Close(image.get(), iterations, connectivity); break; + case MIRTK_VOXEL_GREY: Close(image.get(), iterations, connectivity); break; + case MIRTK_VOXEL_REAL: Close(image.get(), iterations, connectivity); break; + default: + FatalError("Unsupported voxel type: " << ToString(image->GetDataType())); + } + if (verbose) cout << "done" << endl; + + image->Write(output_name); + + return 0; +} + diff --git a/Applications/concatenate-dofs.cc b/Applications/concatenate-dofs.cc new file mode 100644 index 000000000..28cb7ed28 --- /dev/null +++ b/Applications/concatenate-dofs.cc @@ -0,0 +1,270 @@ +/* + * Medical Image Registration ToolKit (MIRTK) + * + * Copyright 2008-2015 Imperial College London + * Copyright 2008-2013 Daniel Rueckert, Julia Schnabel + * Copyright 2013-2015 Andreas Schuh + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include + + +using namespace mirtk; + + +// ============================================================================= +// Help +// ============================================================================= + +// ----------------------------------------------------------------------------- +// Print help screen +void PrintHelp(const char *name) +{ + cout << endl; + cout << "Usage: " << name << " ... [-target ]" << endl; + cout << endl; + cout << "Description:" << endl; + cout << " Computes the composition T of the given input transformations such that" << endl; + cout << endl; + cout << " .. math::" << endl; + cout << endl; + cout << " T(x) = Tn o ... o T2 o T1(x)" << endl; + cout << endl; + cout << "Optional arguments:" << endl; + cout << " -target Target image on which images will be resampled using" << endl; + cout << " the composed transformation." << endl; + PrintStandardOptions(cout); + cout << endl; +} + +// ============================================================================= +// Auxiliaries +// ============================================================================= + +// ----------------------------------------------------------------------------- +/// Compose current composite transformation with any displacement field +void PushDisplacement(FluidFreeFormTransformation &t1, + Transformation *t2, + ImageAttributes &attr) +{ + if (!attr) { + if (t1.NumberOfLevels() == 0) { + MultiLevelTransformation *mffd; + mffd = dynamic_cast(t2); + if (mffd && mffd->NumberOfLevels() > 0) { + attr = mffd->GetLocalTransformation(-1)->Attributes(); + } else { + // Note: As long as no local transformation was encountered before, + // no global transformation or single-level FFD has to be + // approximated by a displacement field. This error should + // thus never be encountered... + cerr << "Internal error: Specify -target image to circumvent it" << endl; + cerr << " and please consider reporting the issue" << endl; + exit(1); + } + } else { + attr = t1.GetLocalTransformation(-1)->Attributes(); + } + } + GenericImage disp(attr, 3); + t2->Displacement(disp); + t1.PushLocalTransformation(new LinearFreeFormTransformation3D(disp)); +} + +// ----------------------------------------------------------------------------- +/// Compose current composite transformation with rigid/affine transformation +void PushTransformation(FluidFreeFormTransformation &t1, + HomogeneousTransformation *t2, + ImageAttributes &attr, + bool last = false) +{ + if (t2->IsIdentity()) return; + if (t1.NumberOfLevels() == 0) { + HomogeneousTransformation *global = t1.GetGlobalTransformation(); + global->PutMatrix(t2->GetMatrix() * global->GetMatrix()); + } else { + if (last) { + AffineTransformation *post = t1.GetAffineTransformation(); + post->PutMatrix(post->GetMatrix() * t2->GetMatrix()); + } else { + PushDisplacement(t1, t2, attr); + } + } +} + +// ----------------------------------------------------------------------------- +/// Compose current composite transformation with single-level FFD +void PushTransformation(FluidFreeFormTransformation &t1, + FreeFormTransformation *t2, + ImageAttributes &attr, + bool = false) +{ + t1.PushLocalTransformation(t2); +} + +// ----------------------------------------------------------------------------- +/// Compose current composite transformation with multi-level FFD +void PushTransformation(FluidFreeFormTransformation &t1, + MultiLevelFreeFormTransformation *t2, + ImageAttributes &attr, + bool = false) +{ + if (t2->NumberOfLevels() == 0) { + PushTransformation(t1, t2->GetGlobalTransformation(), attr); + } else if (t2->NumberOfLevels() == 1) { + t2->MergeGlobalIntoLocalDisplacement(); + t1.PushLocalTransformation(t2->RemoveLocalTransformation(0)); + } else { + PushDisplacement(t1, t2, attr); + } +} + +// ----------------------------------------------------------------------------- +/// Compose current composite transformation with multi-level SV FFD +void PushTransformation(FluidFreeFormTransformation &t1, + MultiLevelStationaryVelocityTransformation *t2, + ImageAttributes &attr, + bool = false) +{ + if (t2->NumberOfLevels() == 0) { + PushTransformation(t1, t2->GetGlobalTransformation(), attr); + } else if (t2->NumberOfLevels() == 1) { + t2->MergeGlobalIntoLocalDisplacement(); + t1.PushLocalTransformation(t2->RemoveLocalTransformation(0)); + } else { + PushDisplacement(t1, t2, attr); + } +} + +// ----------------------------------------------------------------------------- +/// Compose current composite transformation with another composite transformation +void PushTransformation(FluidFreeFormTransformation &t1, + FluidFreeFormTransformation *t2, + ImageAttributes &attr, + bool last = false) +{ + PushTransformation(t1, t2->GetGlobalTransformation(), attr); + for (int l = 0; l < t2->NumberOfLevels(); ++l) { + PushTransformation(t1, t2->RemoveLocalTransformation(l), attr); + } + if (last) { + AffineTransformation *post = t1.GetAffineTransformation(); + post->PutMatrix(post->GetMatrix() * t2->GetAffineTransformation()->GetMatrix()); + } else { + PushTransformation(t1, t2->GetAffineTransformation(), attr); + } +} + +// ============================================================================= +// Main +// ============================================================================= + +// ----------------------------------------------------------------------------- +int main(int argc, char **argv) +{ + REQUIRES_POSARGS(3); + int N = NUM_POSARGS; + + ImageAttributes attr; + FluidFreeFormTransformation t; + + for (ALL_OPTIONS) { + if (OPTION("-target")) { + GreyImage target(ARGUMENT); + attr = target.Attributes(); + } + else HANDLE_STANDARD_OR_UNKNOWN_OPTION(); + } + + HomogeneousTransformation *aff; + FreeFormTransformation *ffd; + MultiLevelFreeFormTransformation *mffd; + MultiLevelStationaryVelocityTransformation *svmffd; + FluidFreeFormTransformation *fluid; + + // Compose affine transformations at the end + if (verbose) { + cout << "Compose affine transformations at end..." << endl; + } + AffineTransformation *post = t.GetAffineTransformation(); + for (int n = N - 1; n >= 1; --n) { + if (verbose) { + cout << "Reading transformation " << n << ": " << POSARG(n) << endl; + } + unique_ptr dof(Transformation::New(POSARG(n))); + if ((aff = dynamic_cast(dof.get()))) { + post->PutMatrix(post->GetMatrix() * aff->GetMatrix()); + --N; // remove transformation from composition chain + } else { + break; + } + } + if (verbose) { + cout << "Compose affine transformations at end... done" << endl; + } + + // Compose remaining transformations from the start + if (verbose) { + cout << "\nCompose remaining transformations..." << endl; + } + for (int n = 1; n < N; ++n) { + + // Read n-th transformation + if (verbose) { + cout << "Reading transformation " << n << " from " << POSARG(n) << endl; + } + unique_ptr dof(Transformation::New(POSARG(n))); + aff = dynamic_cast(dof.get()); + ffd = dynamic_cast(dof.get()); + mffd = dynamic_cast(dof.get()); + svmffd = dynamic_cast(dof.get()); + fluid = dynamic_cast(dof.get()); + + // Compose current composite transformation with n-th transformation + const bool last = (n == N - 1); + if (aff ) PushTransformation(t, aff, attr, last); + else if (mffd ) PushTransformation(t, mffd, attr, last); + else if (svmffd) PushTransformation(t, svmffd, attr, last); + else if (fluid ) PushTransformation(t, fluid, attr, last); + else if (ffd ) { + PushTransformation(t, ffd, attr, last); + dof.release(); + } + else { + cerr << "Unsupported transformation file " << POSARG(n) << endl; + cerr << " Type name = " << dof->NameOfClass() << endl; + exit(1); + } + + // Transform target image attributes by affine transformation such that + // attributes of consecutively approximated displacment fields overlap + // with the thus far transformed image grid + if (aff && attr) { + attr.PutAffineMatrix(aff->GetMatrix(), true); + } + } + if (verbose) { + cout << "Compose remaining transformations... done" << endl; + } + + // Write composite transformation + const char *output_name = POSARG(NUM_POSARGS); + if (verbose) cout << "Writing composite transformation to " << output_name << endl; + t.Write(output_name); + if (verbose) t.Print(2); +} diff --git a/Applications/concatenate-images.cc b/Applications/concatenate-images.cc new file mode 100644 index 000000000..b85dfccf1 --- /dev/null +++ b/Applications/concatenate-images.cc @@ -0,0 +1,301 @@ +/* + * Medical Image Registration ToolKit (MIRTK) + * + * Copyright 2008-2013 Imperial College London + * Copyright 2008-2013 Daniel Rueckert, Julia Schnabel + * Copyright 2013-2015 Andreas Schuh + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +#include + +using namespace mirtk; + + +// ============================================================================= +// Help +// ============================================================================= + +// ----------------------------------------------------------------------------- +void PrintHelp(const char *name) +{ + cout << "\n"; + cout << "Usage: " << name << " ... [options]\n"; + cout << " " << name << " ... -output [options]\n"; + cout << " " << name << " [options]\n"; + cout << "\n"; + cout << "Description:\n"; + cout << " Concatenate two or more either 2D images to form a 3D volume,\n"; + cout << " or 3D volumes to form a 3D+t temporal sequence. All input images\n"; + cout << " must have the same image attributes, except in either the third (2D)\n"; + cout << " or the third and fourth (3D) image dimension.\n"; + cout << "\n"; + cout << "Optional arguments:\n"; + cout << " -output Output image file. Last positional argument if option missing.\n"; + cout << " -image Add one input image. Option can be given multiple times.\n"; + cout << " -images ... Add one or more input images. Option can be given multiple times.\n"; + cout << " -sort [on|off] Automatically determine correct order of input images\n"; + cout << " based on the distance along the third or fourth image axis.\n"; + cout << " -nosort Disable sorting of input images, same as :option:`-sort` off. (default)\n"; + cout << " -start Position of first slice/frame of output volume/sequence.\n"; + cout << " -spacing Slice/Temporal spacing of output volume/sequence.\n"; + PrintStandardOptions(cout); + cout << "\n"; +} + +// ============================================================================= +// Auxiliaries +// ============================================================================= + +// ----------------------------------------------------------------------------- +bool CheckOrientation(const ImageAttributes &a, const ImageAttributes &b) +{ + if ((a._xaxis[0] * b._xaxis[0] + a._xaxis[1] * b._xaxis[1] + a._xaxis[2] * b._xaxis[2] - 1.0) > .001) return false; + if ((a._yaxis[0] * b._yaxis[0] + a._yaxis[1] * b._yaxis[1] + a._yaxis[2] * b._yaxis[2] - 1.0) > .001) return false; + if ((a._zaxis[0] * b._zaxis[0] + a._zaxis[1] * b._zaxis[1] + a._zaxis[2] * b._zaxis[2] - 1.0) > .001) return false; + return true; +} + +// ============================================================================= +// Main +// ============================================================================= + +// ----------------------------------------------------------------------------- +int main(int argc, char *argv[]) +{ + // Parse arguments + REQUIRES_POSARGS(0); + + Array input_names; + for (int i = 1; i <= NUM_POSARGS; ++i) { + input_names.push_back(POSARG(i)); + } + const char *output_name = nullptr; + + bool sort_input = false; + double spacing = numeric_limits::quiet_NaN(); + double start = numeric_limits::quiet_NaN(); + + for (ALL_OPTIONS) { + if (OPTION("-image")) input_names.push_back(ARGUMENT); + else if (OPTION("-images")) { + do { + input_names.push_back(ARGUMENT); + } while (HAS_ARGUMENT); + } + else if (OPTION("-output")) { + output_name = ARGUMENT; + } + else if (OPTION("-sort")) { + if (HAS_ARGUMENT) PARSE_ARGUMENT(sort_input); + else sort_input = true; + } + else if (OPTION("-nosort")) sort_input = false; + else if (OPTION("-spacing")) PARSE_ARGUMENT(spacing); + else if (OPTION("-start")) PARSE_ARGUMENT(start); + else HANDLE_COMMON_OR_UNKNOWN_OPTION(); + } + + if (!output_name) { + if (input_names.empty()) { + FatalError("No input images and -output file specified!"); + } + output_name = input_names.back(); + input_names.pop_back(); + } + + const int n = static_cast(input_names.size()); + if (n < 2) FatalError("Not enough input images given!"); + + // Get attributes of first input image + unique_ptr ref(BaseImage::New(input_names[0])); + if (ref->T() > 1) { + FatalError("Input images must be either 2D (nz == 1) or 3D (nz == 1 and nt == 1)!"); + } + const bool twod = (ref->Z() == 1); + + // Sort input images + // + // Note: Images are not first read all at once into memory to + // not be limited by the amount of memory we have available. + // We better read in the images twice from disk if necessary. + Array pos(n); + for (int i = 0; i < n; ++i) pos[i] = i; + + if (sort_input) { + if (verbose) cout << "Sorting input images...", cout.flush(); + // Compute distance to axis origin + Point p; + Array origin(n); + if (twod) { + p = ref->GetOrigin(); + ref->WorldToImage(p); + origin[0] = p._z; + } else { + origin[0] = ref->GetTOrigin(); + } + for (int i = 1; i < n; ++i) { + unique_ptr image(BaseImage::New(input_names[i])); + p = image->GetOrigin(); + if (twod) { + Point p = image->GetOrigin(); + ref->WorldToImage(p); + origin[i] = p._z; + } else { + origin[i] = image->GetTOrigin(); + } + } + // Determine order of images sorted by ascending distance value + sort(pos.begin(), pos.end(), SortIndicesOfArray(origin)); + if (pos[0] != 0) ref.reset(BaseImage::New(input_names[pos[0]])); + if (verbose) cout << " done\n"; + } + + // Create output image + if (verbose) { + cout << "Making "; + if (twod) cout << "volume"; + else cout << "sequence"; + cout << " from " << n << " input images..."; + } + + ImageAttributes out_attr = ref->Attributes(); + if (twod) out_attr._z = n; // spacing and origin set later below + else out_attr._t = n; + unique_ptr output(BaseImage::New(ref->GetDataType())); + output->Initialize(out_attr); + + // Concatenate images + double avg_spacing = .0; + double min_zorigin = numeric_limits::max(); + + if (twod) { + + for (int j = 0; j < out_attr._y; ++j) + for (int i = 0; i < out_attr._x; ++i) { + output->PutAsDouble(i, j, 0, 0, ref->GetAsDouble(i, j)); + } + for (int k = 1; k < n; ++k) { + unique_ptr image(BaseImage::New(input_names[pos[k]])); + if (image->Z() != 1 || image->T() != 1) { + if (verbose) cout << " failed" << endl; + FatalError("Input image " << k+1 << " is not 2D!"); + } + if (image->X() != ref->X() || image->Y() != ref->Y()) { + if (verbose) cout << " failed" << endl; + FatalError("Input image " << k+1 << " has differing number of voxels!"); + } + if (!fequal(image->GetXSize(), ref->GetXSize()) || + !fequal(image->GetYSize(), ref->GetYSize())) { + if (verbose) cout << " failed" << endl; + FatalError("Input image " << k+1 << " has differing voxel size!"); + } + if (!fequal(image->GetTSize(), ref->GetTSize())) { + Warning("Input image " << k+1 << " has differing temporal spacing."); + } + if (!CheckOrientation(image->Attributes(), ref->Attributes())) { + if (verbose) cout << " failed" << endl; + FatalError("Input image " << k+1 << " has different orientation!"); + } + for (int j = 0; j < out_attr._y; ++j) + for (int i = 0; i < out_attr._x; ++i) { + output->PutAsDouble(i, j, k, 0, image->GetAsDouble(i, j)); + } + min_zorigin = min(min_zorigin, image->GetOrigin()._z); + avg_spacing += image->GetZSize(); + } + + } else { + + for (int k = 0; k < out_attr._z; ++k) + for (int j = 0; j < out_attr._y; ++j) + for (int i = 0; i < out_attr._x; ++i) { + output->PutAsDouble(i, j, k, 0, ref->GetAsDouble(i, j, k)); + } + for (int l = 1; l < n; ++l) { + unique_ptr image(BaseImage::New(input_names[pos[l]])); + if (image->T() != 1) { + if (verbose) cout << " failed" << endl; + FatalError("Input image " << l+1 << " is not 3D!"); + } + if (image->X() != ref->X() || image->Y() != ref->Y() || image->Z() != ref->Z()) { + if (verbose) cout << " failed" << endl; + FatalError("Input image " << l+1 << " has differing number of voxels!"); + } + if (!fequal(image->GetXSize(), ref->GetXSize()) || + !fequal(image->GetYSize(), ref->GetYSize()) || + !fequal(image->GetZSize(), ref->GetZSize())) { + if (verbose) cout << " failed" << endl; + FatalError("Input image " << l+1 << " has differing voxel size!"); + } + if (!CheckOrientation(image->Attributes(), ref->Attributes())) { + if (verbose) cout << " failed" << endl; + FatalError("Input image " << l+1 << " has different orientation!"); + } + for (int k = 0; k < out_attr._z; ++k) + for (int j = 0; j < out_attr._y; ++j) + for (int i = 0; i < out_attr._x; ++i) { + output->PutAsDouble(i, j, k, l, image->GetAsDouble(i, j, k)); + } + avg_spacing += image->GetZSize(); + } + + } + + // Set output spacing to average of input images by default + if (IsNaN(spacing)) spacing = avg_spacing; + + if (twod) { + if (spacing < .0) { // must be positive + out_attr._zaxis[0] *= -1.0; + out_attr._zaxis[1] *= -1.0; + out_attr._zaxis[2] *= -1.0; + spacing *= -1.0; + } + double dx, dy, dz; + output->GetPixelSize(dx, dy, dz); + output->PutPixelSize(dx, dy, spacing); + } else { + output->PutTSize(spacing); // can be negative + } + + // Set output origin + if (twod) { + if (IsNaN(start)) start = min_zorigin; + Point origin = output->GetOrigin(); + origin._z = start; + // Move to center of image stack, i.e. translate along z axis + double s = .5 * (output->Z() * output->GetZSize()); + origin._x += s * output->Attributes()._zaxis[0]; + origin._y += s * output->Attributes()._zaxis[1]; + origin._z += s * output->Attributes()._zaxis[2]; + output->PutOrigin(origin); + } else { + if (!IsNaN(start)) out_attr._torigin = start; + // else torigin is time of first (sorted) input volume + } + + if (verbose) cout << " done" << endl; + + // Write output image + if (verbose) cout << "Writing concatenated images to " << output_name << "..."; + output->Write(output_name); + if (verbose) cout << " done" << endl; + + return 0; +} + diff --git a/Applications/convert-dof.cc b/Applications/convert-dof.cc new file mode 100644 index 000000000..6594a6fe2 --- /dev/null +++ b/Applications/convert-dof.cc @@ -0,0 +1,1476 @@ +/* + * Medical Image Registration ToolKit (MIRTK) + * + * Copyright 2013-2015 Imperial College London + * Copyright 2013-2015 Andreas Schuh + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +#include +#include +#include + +#ifdef HAVE_MIRTK_NIfTI +# include +# include +#endif // HAVE_MIRTK_NIfTI + + +using namespace mirtk; + + +//////////////////////////////////////////////////////////////////////////////// +// Help +//////////////////////////////////////////////////////////////////////////////// + +// ----------------------------------------------------------------------------- +void PrintHelp(const char *name) +{ + cout << "\n"; + cout << "Usage: " << name << " [options]\n"; + cout << " " << name << " -input-format image [options] (3D+3 input image)\n"; + cout << " " << name << " -output-format image [options] (3D+3 output image)\n"; + cout << " " << name << " -input-format image [options] (3x3D input image)\n"; + cout << " " << name << " -output-format image [options] (3x3D output image)\n"; + cout << "\n"; + cout << "Description:\n"; + cout << " Converts between different transformation file formats:\n"; + cout << "\n"; + cout << " ===================== =================================================================================\n"; + cout << " unknown Unknown, try to guess it from file header/type.\n"; + cout << " disp_world|disp|image Dense displacement field image with world space displacement vectors [mm].\n"; + cout << " disp_voxel Dense displacement field image with target space displacement vectors [voxel].\n"; + cout << " mirtk|irtk2 MIRTK/IRTK 2 transformation file format.\n"; + cout << " mirtk_rigid|rigid Rigid MIRTK/IRTK 2 transformation file format (6 DoFs).\n"; + cout << " mirtk_similarity Similarity MIRTK/IRTK 2 transformation file format (7 DoFs).\n"; + cout << " mirtk_affine|affine Affine MIRTK/IRTK 2 transformation file format (12 DoFs).\n"; + cout << " mirtk_linear_ffd Linear free-form deformation.\n"; + cout << " mirtk_bspline_ffd Cubic B-spline free-form deformation.\n"; + cout << " mni_xfm|xfm Linear FreeSurfer transformation (.xfm file).\n"; + cout << " fsl Guess/choose FSL output file format.\n"; + cout << " flirt FSL FLIRT output file format.\n"; + cout << " fnirt FSL FNIRT output file format.\n"; + cout << " nreg Guess/choose Nifty Reg transformation output file format.\n"; + cout << " aladin Nifty Reg Aladin output file format.\n"; + cout << " f3d Nifty Reg reg_f3d output file format with nifti1.intent_p1 code.\n"; + cout << " f3d_def_field Nifty Reg reg_f3d output image deformation field.\n"; + cout << " f3d_disp_field Nifty Reg reg_f3d output image displacement field.\n"; + cout << " f3d_spline_grid Nifty Reg reg_f3d output control point displacement field.\n"; + cout << " f3d_def_vel_field Nifty Reg reg_f3d output image deformation field as stationary velocity field.\n"; + cout << " f3d_disp_vel_field Nifty Reg reg_f3d output image displacement field as stationary velocity field.\n"; + cout << " f3d_spline_vel_grid Nifty Reg reg_f3d output control point velocity field.\n"; + cout << " ===================== =================================================================================\n"; +#ifndef HAVE_MIRTK_NIfTI + cout << "\n"; + cout << " Note: Cannot convert from/to the following formats because of missing NIfTI module.\n"; + cout << " Rebuild the MIRTK with the NIfTI module enabled to use these formats.\n"; + cout << "\n"; + cout << " Not available: f3d*, fnirt\n"; + cout << "\n"; +#endif // HAVE_MIRTK_NIfTI + cout << "\n"; + cout << "Arguments:\n"; + cout << " input Input transformation file.\n"; + cout << " output Output transformation file.\n"; + cout << "\n"; + cout << "Options:\n"; + cout << " -input-format Format of input file. (default: unknown)\n"; + cout << " -output-format Format of output file. (default: mirtk)\n"; + cout << " -format Short for :option:`-output-format`.\n"; + cout << " -dofin Affine transformation component in MIRTK format which\n"; + cout << " shall be removed from Nifty Reg's FFD such that the\n"; + cout << " resulting transformation is a MIRTK MFFD consisting\n"; + cout << " of this affine component plus the remaining non-linear\n"; + cout << " component of the input FFD. (default: none)\n"; + cout << " -target Target image. Required for from/to FSL format conversion.\n"; + cout << " Also required when converting linear transformations to displacment\n"; + cout << " fields or (M)FFD transformation types. (default: none)\n"; + cout << " -source Source image. Required for from/to FSL format conversion. (default: none)\n"; + cout << " -Ts