1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157
|
/* -*- mia-c++ -*-
*
* This file is part of MIA - a toolbox for medical image analysis
* Copyright (c) Leipzig, Madrid 1999-2014 Gert Wollny
*
* MIA is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with MIA; if not, see <http://www.gnu.org/licenses/>.
*
*/
#define VSTREAM_DOMAIN "3dmyoserial"
#include <fstream>
#include <libxml++/libxml++.h>
#include <itpp/signal/fastica.h>
#include <boost/filesystem.hpp>
#include <mia/core/filetools.hh>
#include <mia/core/msgstream.hh>
#include <mia/core/cmdlineparser.hh>
#include <mia/core/errormacro.hh>
#include <mia/3d/nonrigidregister.hh>
#include <mia/3d/transformfactory.hh>
#include <mia/3d/imageio.hh>
#include <mia/internal/main.hh>
using namespace std;
using namespace mia;
namespace bfs=boost::filesystem;
const SProgramDescription g_general_help = {
{pdi_group, "Registration of series of 3D images"},
{pdi_short, "Serial registration of 3D images."},
{pdi_description, "This program runs the image registration of a consecutively numbered image series. "
"The registration is run in a serial manner, this is, only images in "
"temporal succession (i.e. consecutive numbers) are registered, and the obtained transformations "
"are applied accumulated to reach full registration. "},
{pdi_example_descr, "Run a serial registration of images inputXXXX.v (X digit) to reference "
"image 20 and store the result in regXXXX.v. Optimize the sum of squared differences "
"and spline transformations with coefficient rate 10."},
{pdi_example_code, "-i input0000.v -o 'reg%04d.v' -f spline:rate=10 -r 20 ssd"}
};
int do_main( int argc, char *argv[] )
{
// IO parameters
string in_filename;
string registered_filebase("reg%04d.v");
P3DTransformationFactory transform_creator;
// registration parameters
PMinimizer minimizer;
size_t mg_levels = 3;
int reference_param = -1;
CCmdOptionList options(g_general_help);
options.set_group("\nFile-IO");
options.add(make_opt( in_filename, "in-file", 'i', "input perfusion data set", CCmdOptionFlags::required_input));
options.add(make_opt( registered_filebase, "out-file", 'o', "file name for registered fiels", CCmdOptionFlags::output));
options.set_group("\nRegistration");
options.add(make_opt( minimizer, "gsl:opt=gd,step=0.1", "optimizer", 'O', "Optimizer used for minimization"));
options.add(make_opt( mg_levels, "mg-levels", 'l', "multi-resolution levels"));
options.add(make_opt( transform_creator, "spline", "transForm", 'f', "transformation type"));
options.add(make_opt( reference_param, "ref", 'r', "reference frame (-1 == use image in the middle)"));
if (options.parse(argc, argv, "cost", &C3DFullCostPluginHandler::instance()) != CCmdOptionList::hr_no)
return EXIT_SUCCESS;
// create cost function chain
auto cost_functions = options.get_remaining();
if (cost_functions.empty())
throw invalid_argument("No cost function given - nothing to register");
C3DFullCostList costs;
for (auto c = cost_functions.begin(); c != cost_functions.end(); ++c) {
auto cost = C3DFullCostPluginHandler::instance().produce(*c);
costs.push(cost);
}
size_t start_filenum = 0;
size_t end_filenum = 0;
size_t format_width = 0;
string src_basename = get_filename_pattern_and_range(in_filename, start_filenum, end_filenum, format_width);
C3DImageSeries input_images;
for (size_t i = start_filenum; i < end_filenum; ++i) {
string src_name = create_filename(src_basename.c_str(), i);
P3DImage image = load_image<P3DImage>(src_name);
if (!image)
throw create_exception<runtime_error>( "image ", src_name, " not found");
cvdebug() << "read '" << src_name << "\n";
input_images.push_back(image);
}
// if reference is not given, use half range
size_t reference = reference_param < 0 ? input_images.size() / 2 : reference_param;
// prepare registration framework
C3DNonrigidRegister nrr(costs, minimizer, transform_creator, mg_levels);
if ( input_images.empty() )
throw invalid_argument("No input images to register");
if (reference > input_images.size() - 1) {
reference = input_images.size() - 1;
cvwarn() << "Reference was out of range, adjusted to " << reference << "\n";
}
// run forward registrations
for (size_t i = 0; i < reference; ++i) {
P3DTransformation transform = nrr.run(input_images[i], input_images[i+1]);
for (size_t j = 0; j <=i ; ++j) {
input_images[j] = (*transform)(*input_images[j]);
}
}
// run backward registration
for (size_t i = input_images.size() - 1; i > reference; --i) {
P3DTransformation transform = nrr.run(input_images[i], input_images[i-1]);
for (size_t j = input_images.size() - 1; j >= i ; --j) {
input_images[j] = (*transform)(*input_images[j]);
}
}
bool success = true;
auto ii = input_images.begin();
for (size_t i = start_filenum; i < end_filenum; ++i, ++ii) {
string out_name = create_filename(registered_filebase.c_str(), i);
cvmsg() << "Save image " << i << " to " << out_name << "\n";
success &= save_image(out_name, *ii);
}
return success ? EXIT_SUCCESS : EXIT_FAILURE;
}
MIA_MAIN(do_main);
|