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 "gtest/gtest.h"
21
22 #include "mirtk/GenericImage.h"
23 #include "mirtk/ConvolutionFunction.h"
24 #include "mirtk/ScalarGaussian.h"
25 #include "mirtk/ScalarFunctionToImage.h"
26
27 using namespace mirtk;
28 using namespace mirtk::ConvolutionFunction;
29
30 static const char *test_file = NULL;
31 static const char *result_file = NULL;
32 static double sigma = 1.0;
33 static double bgvalue = 0.0;
34
35 // ===========================================================================
36 // Kernels
37 // ===========================================================================
38
Gaussian(double sigma,double dx)39 GenericImage<RealPixel> Gaussian(double sigma, double dx)
40 {
41 sigma = sigma / dx;
42 ScalarGaussian gaussian(sigma, 1, 1, 0, 0, 0);
43 GenericImage<RealPixel> kernel(2 * int(round(4.0 * sigma)) + 1, 1, 1);
44 ScalarFunctionToImage<RealPixel> generator;
45 generator.Input (&gaussian);
46 generator.Output(&kernel);
47 generator.Run();
48 RealPixel sum = .0;
49 for (int i = 0; i < kernel.X(); ++i) sum += kernel(i, 0, 0, 0);
50 for (int i = 0; i < kernel.X(); ++i) kernel(i, 0, 0, 0) /= sum;
51 return kernel;
52 }
53
54 // ===========================================================================
55 // Tests
56 // ===========================================================================
57
58 // ---------------------------------------------------------------------------
TEST(ConvolutionFunction,ConvolveMirroredForeground)59 TEST(ConvolutionFunction, ConvolveMirroredForeground)
60 {
61 GenericImage<RealPixel> image (test_file);
62 GenericImage<RealPixel> output(image.GetImageAttributes());
63 GenericImage<RealPixel> xkernel = Gaussian(sigma, image.GetXSize());
64 GenericImage<RealPixel> ykernel = Gaussian(sigma, image.GetYSize());
65 GenericImage<RealPixel> zkernel = Gaussian(sigma, image.GetZSize());
66 image.PutBackgroundValueAsDouble(bgvalue, true);
67 ConvolveMirroredForegroundInX<RealPixel> xconv(&image, &xkernel);
68 ParallelForEachVoxel(image.GetImageAttributes(), image, output, xconv);
69 ConvolveMirroredForegroundInY<RealPixel> yconv(&image, &ykernel);
70 ParallelForEachVoxel(image.GetImageAttributes(), output, image, yconv);
71 ConvolveMirroredForegroundInZ<RealPixel> zconv(&image, &zkernel);
72 ParallelForEachVoxel(image.GetImageAttributes(), image, output, zconv);
73 output.Write(result_file);
74 }
75
76 // ===========================================================================
77 // Main
78 // ===========================================================================
79
80 // ---------------------------------------------------------------------------
main(int argc,char * argv[])81 int main(int argc, char *argv[])
82 {
83 ::testing::InitGoogleTest(&argc, argv);
84 if (argc < 3 || argc > 5) {
85 cerr << "usage: " << argv[0] << " <input> <output> [sigma] [background]" << endl;
86 exit(1);
87 }
88 test_file = argv[1];
89 result_file = argv[2];
90 if (argc > 3) sigma = atof(argv[3]);
91 if (argc > 4) bgvalue = atof(argv[4]);
92 return RUN_ALL_TESTS();
93 }
94