daffodil.filter.gaussian source code

1 /**
2  * A gaussian filter (aka gaussian blur) is a convolution
3  * (:d:mod:`daffodil.filter.convolve`) using a matrix created from a gaussian
4  * distribution.
5  */
6 module daffodil.filter.gaussian;
7 
8 import std.math;
9 
10 import daffodil.image;
11 import daffodil.filter;
12 
13 /**
14  * Evaluate the gaussian/normal distribution for a given ``x``, ``stDev`` and
15  * ``mean``.
16  */
17 real gaussianDistribution(real x, real stDev = 1, real mean = 0) {
18     return 1/(stDev * sqrt(2 * PI))* E.pow(-pow(x - mean, 2)/(2 * pow(stDev, 2)));
19 }
20 
21 @("gaussian distribution")
22 unittest {
23     static void assertEq(real a, real b) {
24         return assert(approxEqual(a, b));
25     }
26 
27     // Standard Normal Distribution
28     assertEq(gaussianDistribution(0), 0.398942);
29     assertEq(gaussianDistribution(1), 0.241971);
30     assertEq(gaussianDistribution(1), gaussianDistribution(-1));
31     assertEq(gaussianDistribution(2), 0.053991);
32     assertEq(gaussianDistribution(2), gaussianDistribution(-2));
33     assertEq(gaussianDistribution(3), 0.004432);
34     assertEq(gaussianDistribution(3), gaussianDistribution(-3));
35 
36     // Changed standard deviation
37     assertEq(gaussianDistribution(0, 2), 0.199471);
38     assertEq(gaussianDistribution(2, 2), 0.120985);
39 
40     // Changed mean
41     assertEq(gaussianDistribution(0, 2, 2), 0.120985);
42     assertEq(gaussianDistribution(2, 2, 2), 0.199471);
43     assertEq(gaussianDistribution(4, 2, 2), 0.120985);
44 }
45 
46 /**
47  * Create a 1D matrix of a discrete gaussian distribution with a given standard
48  * deviation and the number of standard deviations to stop generating at. The
49  * result is mirrored with guaranteed odd length.
50  *
51  * The result can be used to convolve a image.
52  */
53 real[] gaussianMatrix(real stDev = 1, real maxDev = 3) {
54     auto range = cast(uint)ceil(stDev * maxDev);
55     auto ret = new real[1 + 2*range];
56 
57     ret[range] = gaussianDistribution(0, stDev);
58     foreach (i; 1..range + 1) {
59         ret[range + i] = ret[range - i] = gaussianDistribution(i, stDev);
60     }
61 
62     return ret;
63 }
64 
65 @("gaussian matrix")
66 unittest {
67     const matrix = [0.004432, 0.053991, 0.241971, 0.398942, 0.241971, 0.053991, 0.004432];
68     assert(approxEqual(gaussianMatrix(), matrix));
69 
70     assert(gaussianMatrix(10).length == 61);
71 }
72 
73 /**
74  * Return a copy of ``image`` with a gaussian blur applied across axies
75  * ``axis`` with a given standard deviation and the number of standard
76  * deviations to stop at.
77  */
78 auto gaussianBlurred(string axis = "xy", V)(const Image!V image, real stDev = 1, real maxDev = 3) {
79     auto matrix = gaussianMatrix(stDev, maxDev);
80 
81     return image.convolved!axis(matrix);
82 }
83 
84 @("gaussian blur")
85 unittest {
86     import daffodil;
87 
88     auto image = new Image!ubyte(2, 2, 3, RGB);
89     image[0, 0] = [1f, 1f, 1f];
90     image[0, 1] = [1f, 0f, 0f];
91     image[1, 0] = [0f, 1f, 0f];
92     image[1, 1] = [0f, 0f, 1f];
93 
94     auto blurred = image.gaussianBlurred(2);
95     assert(blurred[0, 0] == [130, 130, 122]);
96     assert(blurred[0, 1] == [130, 117, 122]);
97     assert(blurred[1, 0] == [114, 130, 122]);
98     assert(blurred[1, 1] == [114, 117, 122]);
99 }