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
29 using namespace mirtk;
30
31
32 // =============================================================================
33 // Global variables with default arguments
34 // =============================================================================
35
36 char IntegrationMethod[32] = {"SS"};
37 int NumberOfSteps = 128;
38
39 // =============================================================================
40 // Help
41 // =============================================================================
42
43 // -----------------------------------------------------------------------------
PrintHelp(const char * name)44 void PrintHelp(const char *name)
45 {
46 cout << endl;
47 cout << "Usage: " << name << " <velo> <disp> [options]" << endl;
48 cout << " " << name << " <vx> <vy> [<vz>] <dx> <dy> [<dz>] [options]" << endl;
49 cout << endl;
50 cout << "Description:" << endl;
51 cout << " Reads a dense 3D (2D) stationary velocity field, computes the corresponding" << endl;
52 cout << " dense 3D (2D) displacement field, and writes the resulting vector field." << endl;
53 cout << endl;
54 cout << "Options:" << endl;
55 cout << " -ss Use scaling and squaring (SS) method. (default)" << endl;
56 cout << " -euler Use forward Euler method." << endl;
57 cout << " -steps <int> Number of integration steps (2^n in case of SS). "
58 << "(default: " << NumberOfSteps << ")" << endl;
59 PrintCommonOptions(cout);
60 cout << endl;
61 }
62
63 // =============================================================================
64 // Main
65 // =============================================================================
66
67 // -----------------------------------------------------------------------------
main(int argc,char ** argv)68 int main(int argc, char **argv)
69 {
70 // Print help or version and exit if requested
71 HANDLE_HELP_OR_VERSION();
72
73 InitializeIOLibrary();
74
75 // Parse positional arguments
76 const char *velo_fname = NULL;
77 const char *disp_fname = NULL;
78 const char *vx_fname = NULL;
79 const char *vy_fname = NULL;
80 const char *vz_fname = NULL;
81 const char *dx_fname = NULL;
82 const char *dy_fname = NULL;
83 const char *dz_fname = NULL;
84
85 int posargc = 0;
86 while (posargc+1 < argc && argv[posargc+1][0] != '-') posargc++;
87
88 switch (posargc) {
89 case 2:
90 velo_fname = argv[1];
91 disp_fname = argv[2];
92 break;
93 case 4:
94 vx_fname = argv[1];
95 vy_fname = argv[2];
96 dx_fname = argv[3];
97 dy_fname = argv[4];
98 break;
99 case 6:
100 vx_fname = argv[1];
101 vy_fname = argv[2];
102 vz_fname = argv[3];
103 dx_fname = argv[4];
104 dy_fname = argv[5];
105 dz_fname = argv[6];
106 break;
107 default:
108 PrintHelp(EXECNAME);
109 exit(1);
110 }
111
112 // Parse optional arguments
113 for (ARGUMENTS_AFTER(posargc)) {
114 if (OPTION("-int")) {
115 const char *arg = ARGUMENT;
116 size_t len = strlen(arg);
117 if (len > 99u) len = 99u;
118 for (size_t i = 0u; i < len; i++) {
119 if ('a' <= arg[i] && arg[i] <= 'z') IntegrationMethod[i] = 'A' + (arg[i] - 'a');
120 else IntegrationMethod[i] = arg[i];
121 }
122 IntegrationMethod[len] = '\0';
123 } else if (OPTION("-ss")) {
124 IntegrationMethod[0] = 'S';
125 IntegrationMethod[1] = 'S';
126 IntegrationMethod[2] = '\0';
127 } else if (OPTION("-euler")) {
128 IntegrationMethod[0] = 'R';
129 IntegrationMethod[1] = 'K';
130 IntegrationMethod[2] = 'E';
131 IntegrationMethod[3] = '1';
132 IntegrationMethod[4] = '\0';
133 } else if (OPTION("-steps")) {
134 NumberOfSteps = atoi(ARGUMENT);
135 } else {
136 HANDLE_COMMON_OR_UNKNOWN_OPTION();
137 }
138 }
139
140 // Read input velocities
141 UniquePtr<RealImage> v;
142
143 if (velo_fname) {
144 v.reset(new RealImage(velo_fname));
145 } else {
146 RealImage vx(vx_fname);
147 RealImage vy(vy_fname);
148 UniquePtr<RealImage> vz(vz_fname ? new RealImage(vz_fname) : nullptr);
149
150 ImageAttributes attr = vx.Attributes();
151 if (attr._t > 1) {
152 FatalError("Separate input component images must be scalar velocity fields.\n"
153 "Try calling this program with single vector field input and output images instead.");
154 }
155 if (vy.Attributes() != attr || (vz.get() && vz->Attributes() != attr)) {
156 FatalError("Image attributes of input velocity component images do not match!");
157 }
158 attr._t = (vz.get() ? 3 : 2);
159
160 v.reset(new RealImage(attr));
161 for (int k = 0; k < attr._z; ++k)
162 for (int j = 0; j < attr._y; ++j)
163 for (int i = 0; i < attr._x; ++i) {
164 v->Put(i, j, k, 0, vx.Get(i, j, k));
165 v->Put(i, j, k, 1, vy.Get(i, j, k));
166 if (vz.get()) v->Put(i, j, k, 2, vz->Get(i, j, k));
167 }
168 }
169
170 // Instantiate filter which implements the desired integration method
171 UniquePtr<VelocityToDisplacementField<RealPixel> > vtod;
172
173 if (strcmp(IntegrationMethod, "SS") == 0) {
174 vtod.reset(new VelocityToDisplacementFieldSS<RealPixel>());
175 } else if (strcmp(IntegrationMethod, "RKE1") == 0 || strcmp(IntegrationMethod, "FORWARDEULER") == 0 || strcmp(IntegrationMethod, "EULER") == 0) {
176 vtod.reset(new VelocityToDisplacementFieldEuler<RealPixel>());
177 } else {
178 FatalError("Unknown integration method: " << IntegrationMethod);
179 }
180
181 // Compute displacement field
182 RealImage d;
183
184 vtod->Input(v.get());
185 vtod->Output(&d);
186 vtod->NumberOfSteps(NumberOfSteps);
187 vtod->Run();
188
189 // Save output velocities
190 if (disp_fname) {
191 d.Write(disp_fname);
192 } else {
193 d.GetFrame(0).Write(dx_fname);
194 d.GetFrame(1).Write(dy_fname);
195 if (dz_fname) d.GetFrame(2).Write(dz_fname);
196 }
197
198 return 0;
199 }
200