1 /*
2 * Medical Image Registration ToolKit (MIRTK)
3 *
4 * Copyright 2013-2015 Imperial College London
5 * Copyright 2013-2015 Andreas Schuh
6 *
7 * Licensed under the Apache License, Version 2.0 (the "License");
8 * you may not use this file except in compliance with the License.
9 * You may obtain a copy of the License at
10 *
11 * http://www.apache.org/licenses/LICENSE-2.0
12 *
13 * Unless required by applicable law or agreed to in writing, software
14 * distributed under the License is distributed on an "AS IS" BASIS,
15 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16 * See the License for the specific language governing permissions and
17 * limitations under the License.
18 */
19
20 #include "mirtk/Common.h"
21 #include "mirtk/Options.h"
22
23 #include "mirtk/IOConfig.h"
24
25 #include "mirtk/VelocityToDisplacementField.h"
26 #include "mirtk/VelocityToDisplacementFieldEuler.h"
27 #include "mirtk/VelocityToDisplacementFieldSS.h"
28 #include "mirtk/DisplacementToVelocityFieldBCH.h"
29
30 using namespace mirtk;
31
32
33 // =============================================================================
34 // Global variables with default arguments
35 // =============================================================================
36
37 char IntegrationMethod[32] = {"SS"};
38 int NumberOfSteps = 64;
39 int NumberOfBCHSteps = 8;
40 int NumberOfBCHTerms = 6;
41 bool SmoothBCHApproximation = false;
42 bool UseJacobian = false;
43
44 // =============================================================================
45 // Help
46 // =============================================================================
47
48 // -----------------------------------------------------------------------------
PrintHelp(const char * name)49 void PrintHelp(const char *name)
50 {
51 cout << endl;
52 cout << "Usage: " << name << " <disp> <velo> [options]" << endl;
53 cout << " " << name << " <dx> <dy> [<dz>] <vx> <vy> [<vz>] [options]" << endl;
54 cout << endl;
55 cout << "Description:" << endl;
56 cout << " Reads a dense 3D (2D) displacement field, computes the corresponding" << endl;
57 cout << " stationary 3D (2D) velocity field, and writes the resulting vector field." << endl;
58 cout << endl;
59 cout << "Options:" << endl;
60 cout << " -ss Use scaling and squaring (SS) method for exponentiation. (default)" << endl;
61 cout << " -euler Use forward Euler method for exponentiation." << endl;
62 cout << " -steps <int> Number of integration steps (2^n in case of SS). "
63 << "(default: " << NumberOfSteps << ")" << endl;
64 cout << " -terms <int> Number of BCH approximation terms (either 2 or 3). (default: " << NumberOfBCHTerms << ")" << endl;
65 cout << " -iters <int> Number of BCH update steps. (default: " << NumberOfBCHSteps << ")" << endl;
66 cout << " -smooth Smooth velocities before each update. (default: off)" << endl;
67 cout << " -jac Use Jacobian for Lie bracket computation. (default: off)" << endl;
68 PrintCommonOptions(cout);
69 cout << endl;
70 }
71
72 // =============================================================================
73 // Main
74 // =============================================================================
75
76 // -----------------------------------------------------------------------------
main(int argc,char * argv[])77 int main(int argc, char *argv[])
78 {
79 // Print help or version and exit if requested
80 HANDLE_HELP_OR_VERSION();
81
82 InitializeIOLibrary();
83
84 // Parse positional arguments
85 const char *disp_fname = NULL;
86 const char *velo_fname = NULL;
87 const char *dx_fname = NULL;
88 const char *dy_fname = NULL;
89 const char *dz_fname = NULL;
90 const char *vx_fname = NULL;
91 const char *vy_fname = NULL;
92 const char *vz_fname = NULL;
93
94 int posargc = 0;
95 while (posargc+1 < argc && argv[posargc+1][0] != '-') posargc++;
96
97 switch (posargc) {
98 case 2:
99 disp_fname = argv[1];
100 velo_fname = argv[2];
101 break;
102 case 4:
103 dx_fname = argv[1];
104 dy_fname = argv[2];
105 vx_fname = argv[3];
106 vy_fname = argv[4];
107 break;
108 case 6:
109 dx_fname = argv[1];
110 dy_fname = argv[2];
111 dz_fname = argv[3];
112 vx_fname = argv[4];
113 vy_fname = argv[5];
114 vz_fname = argv[6];
115 break;
116 default:
117 PrintHelp(EXECNAME);
118 exit(1);
119 }
120
121 // Parse optional arguments
122 for (ARGUMENTS_AFTER(posargc)) {
123 if (OPTION("-steps" )) NumberOfSteps = atoi(ARGUMENT);
124 else if (OPTION("-iters" )) NumberOfBCHSteps = atoi(ARGUMENT);
125 else if (OPTION("-terms" )) NumberOfBCHTerms = atoi(ARGUMENT);
126 else if (OPTION("-smooth")) SmoothBCHApproximation = true;
127 else if (OPTION("-jac")) UseJacobian = true;
128 else if (OPTION("-int")) {
129 const char *arg = ARGUMENT;
130 size_t len = strlen(arg);
131 if (len > 99u) len = 99u;
132 for (size_t i = 0u; i < len; ++i) {
133 if ('a' <= arg[i] && arg[i] <= 'z') IntegrationMethod[i] = 'A' + (arg[i] - 'a');
134 else IntegrationMethod[i] = arg[i];
135 }
136 IntegrationMethod[len] = '\0';
137 } else if (OPTION("-ss")) {
138 IntegrationMethod[0] = 'S';
139 IntegrationMethod[1] = 'S';
140 IntegrationMethod[2] = '\0';
141 } else if (OPTION("-euler")) {
142 IntegrationMethod[0] = 'R';
143 IntegrationMethod[1] = 'K';
144 IntegrationMethod[2] = 'E';
145 IntegrationMethod[3] = '1';
146 IntegrationMethod[4] = '\0';
147 } else {
148 HANDLE_COMMON_OR_UNKNOWN_OPTION();
149 }
150 }
151
152 // Read input displacements
153 UniquePtr<RealImage> d;
154
155 if (disp_fname) {
156 d.reset(new RealImage(disp_fname));
157 } else {
158 RealImage dx(dx_fname);
159 RealImage dy(dy_fname);
160 UniquePtr<RealImage> dz(dz_fname ? new RealImage(dz_fname) : nullptr);
161
162 ImageAttributes attr = dx.Attributes();
163 if (attr._t > 1) {
164 FatalError("Separate input component images must be scalar displacement fields.\n"
165 "Try calling this program with single vector field input and output images instead.");
166 }
167 if (dy.Attributes() != attr || (dz.get() && dz->Attributes() != attr)) {
168 FatalError("Image attributes of input displacement component images do not match!");
169 }
170
171 d.reset(new RealImage(attr, dz.get() ? 3 : 2));
172 for (int k = 0; k < attr._z; ++k)
173 for (int j = 0; j < attr._y; ++j)
174 for (int i = 0; i < attr._x; ++i) {
175 d->Put(i, j, k, 0, dx.Get(i, j, k));
176 d->Put(i, j, k, 1, dy.Get(i, j, k));
177 if (dz.get()) d->Put(i, j, k, 2, dz->Get(i, j, k));
178 }
179 }
180
181 // Instantiate exponential filter
182 UniquePtr<VelocityToDisplacementField<RealPixel> > exp;
183 if (strcmp(IntegrationMethod, "SS") == 0) {
184 exp.reset(new VelocityToDisplacementFieldSS<RealPixel>());
185 } else if (strcmp(IntegrationMethod, "RKE1") == 0 ||
186 strcmp(IntegrationMethod, "FORWARDEULER") == 0 ||
187 strcmp(IntegrationMethod, "EULER") == 0) {
188 exp.reset(new VelocityToDisplacementFieldEuler<RealPixel>());
189 } else {
190 FatalError("Unknown integration method: " << IntegrationMethod);
191 }
192
193 // Compute stationary velocity field
194 DisplacementToVelocityFieldBCH<RealPixel> dtov;
195 RealImage v;
196
197 dtov.ExponentialFilter (exp.get());
198 dtov.NumberOfSteps (NumberOfSteps);
199 dtov.NumberOfIterations(NumberOfBCHSteps);
200 dtov.NumberOfTerms (NumberOfBCHTerms);
201 dtov.SmoothVelocities (SmoothBCHApproximation);
202 dtov.UseJacobian (UseJacobian);
203
204 dtov.Input(d.get());
205 dtov.Output(&v);
206
207 dtov.Run();
208
209 // Save output velocities
210 if (velo_fname) {
211 v.Write(velo_fname);
212 } else {
213 v.GetFrame(0).Write(vx_fname);
214 v.GetFrame(1).Write(vy_fname);
215 if (vz_fname) v.GetFrame(2).Write(vz_fname);
216 }
217
218 return 0;
219 }
220