SMIL  0.9.1
DImageConvolution.hpp
1 /*
2  * Copyright (c) 2011-2016, Matthieu FAESSEL and ARMINES
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions are met:
7  *
8  * * Redistributions of source code must retain the above copyright
9  * notice, this list of conditions and the following disclaimer.
10  * * Redistributions in binary form must reproduce the above copyright
11  * notice, this list of conditions and the following disclaimer in the
12  * documentation and/or other materials provided with the distribution.
13  * * Neither the name of Matthieu FAESSEL, or ARMINES nor the
14  * names of its contributors may be used to endorse or promote products
15  * derived from this software without specific prior written permission.
16  *
17  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS'' AND ANY
18  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
19  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
20  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS AND CONTRIBUTORS BE LIABLE FOR ANY
21  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
22  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
23  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
24  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
26  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27  */
28 
29 
30 #ifndef _D_IMAGE_CONVOLUTION_HPP
31 #define _D_IMAGE_CONVOLUTION_HPP
32 
33 #include "DLineArith.hpp"
34 #include "Core/include/private/DBufferPool.hpp"
35 
36 
37 namespace smil
38 {
39 
60  // Inplace safe
61  template <class T>
62  RES_T horizConvolve(const Image<T> &imIn, const vector<double> &kernel, Image<T> &imOut)
63  {
64  CHECK_ALLOCATED(&imIn, &imOut);
65  CHECK_SAME_SIZE(&imIn, &imOut);
66 
67  ImageFreezer freeze(imOut);
68 
69  typename ImDtTypes<T>::sliceType linesIn = imIn.getLines();
70  typename ImDtTypes<T>::sliceType linesOut = imOut.getLines();
71  typename ImDtTypes<T>::lineType lOut;
72 
73  int imW = imIn.getWidth();
74  int kernelRadius = (kernel.size()-1)/2;
75 // int kLen = 2*kernelRadius+1;
76 
77  double *partialKernWeights = new double[kernelRadius];
78  double pkwSum = 0;
79  for (int i=0;i<kernelRadius;i++)
80  pkwSum += kernel[i];
81  for (int i=0;i<kernelRadius;i++)
82  {
83  pkwSum += kernel[i+kernelRadius];
84  partialKernWeights[i] = pkwSum;
85  }
86 
87  typedef double bufType; // If float, loops are vectorized
88  BufferPool<bufType> bufferPool(imW);
89 
90  #ifdef USE_OPEN_MP
91  int nthreads = Core::getInstance()->getNumberOfThreads();
92  #pragma omp parallel private(lOut) num_threads(nthreads)
93  #endif // USE_OPEN_MP
94  {
95  typename ImDtTypes<bufType>::lineType lIn = bufferPool.getBuffer();
96  double sum;
97 
98  #ifdef USE_OPEN_MP
99  #pragma omp for
100  #endif // USE_OPEN_MP
101  for (int y=0;y<(int)imIn.getLineCount();y++)
102  {
103  copyLine<T,bufType>(linesIn[y], imW, lIn);
104  lOut = linesOut[y];
105 
106  // left pixels
107  for (int x=0;x<kernelRadius;x++)
108  {
109  sum = 0;
110  for (int i=-x;i<=kernelRadius;i++)
111  sum += kernel[i+kernelRadius]*lIn[x+i];
112  lOut[x] = T(sum / partialKernWeights[x]);
113 
114  }
115 
116  // center pixels
117  for (int x=kernelRadius;x<imW-kernelRadius;x++)
118  {
119  sum = 0;
120  for (int i=-kernelRadius;i<=kernelRadius;i++)
121  sum += kernel[i+kernelRadius]*lIn[x+i];
122  lOut[x] = T(sum);
123  }
124 
125  // right pixels
126  for (int x=imW-kernelRadius;x<imW;x++)
127  {
128  sum = 0;
129  for (int i=-kernelRadius;i<imW-x;i++)
130  sum += kernel[i+kernelRadius]*lIn[x+i];
131  lOut[x] = T(sum / partialKernWeights[imW-1-x]);
132  }
133  }
134 
135  }
136  delete[] partialKernWeights;
137  return RES_OK;
138  }
139 
145  template <class T>
146  RES_T vertConvolve(const Image<T> &imIn, const vector<double> &kernel, Image<T> &imOut)
147  {
148  CHECK_ALLOCATED(&imIn, &imOut);
149  CHECK_SAME_SIZE(&imIn, &imOut);
150 
151  ImageFreezer freeze(imOut);
152 
153  typename ImDtTypes<T>::volType slicesIn = imIn.getSlices();
154  typename ImDtTypes<T>::volType slicesOut = imOut.getSlices();
155  typename ImDtTypes<T>::sliceType sIn, sOut;
156 
157  int imW = imIn.getWidth();
158  int imH = imIn.getHeight();
159  int imD = imIn.getDepth();
160  int kernelRadius = (kernel.size()-1)/2;
161 
162  double *partialKernWeights = new double[kernelRadius];
163  double pkwSum = 0;
164  for (int i=0;i<kernelRadius;i++)
165  pkwSum += kernel[i];
166  for (int i=0;i<kernelRadius;i++)
167  {
168  pkwSum += kernel[i+kernelRadius];
169  partialKernWeights[i] = pkwSum;
170  }
171 
172  typedef double bufType; // If double, loops are vectorized
173  BufferPool<bufType> bufferPool(imW);
174 
175  for (int z=0;z<imD;z++)
176  {
177  #ifdef USE_OPEN_MP
178  int nthreads = Core::getInstance()->getNumberOfThreads();
179  #pragma omp parallel private(sIn, sOut) num_threads(nthreads)
180  #endif // USE_OPEN_MP
181  {
182  sIn = slicesIn[z];
183  sOut = slicesOut[z];
184  double sum;
185 
186  #ifdef USE_OPEN_MP
187  #pragma omp for
188  #endif // USE_OPEN_MP
189  for (int x=0;x<imW;x++)
190  {
191  // Top pixels
192  for (int y=0;y<kernelRadius;y++)
193  {
194  sum = 0;
195  for (int i=-y;i<kernelRadius+1;i++)
196  sum += kernel[i+kernelRadius]*sIn[y+i][x];
197  sOut[y][x] = T(sum / partialKernWeights[y]);
198  }
199 
200  // Center pixels
201  for (int y=kernelRadius;y<imH-kernelRadius;y++)
202  {
203  sum = 0;
204  for (int i=-kernelRadius;i<=kernelRadius;i++)
205  sum += kernel[i+kernelRadius]*sIn[y+i][x];
206  sOut[y][x] = T(sum);
207  }
208 
209  // Bottom pixels
210  for (int y=imH-kernelRadius;y<imH;y++)
211  {
212  sum = 0;
213  for (int i=-kernelRadius;i<imH-y;i++)
214  sum += kernel[i+kernelRadius]*sIn[y+i][x];
215  sOut[y][x] = T(sum / partialKernWeights[imH-1-y]);
216  }
217  }
218  }
219  }
220  delete[] partialKernWeights;
221  return RES_OK;
222  }
223 
229  template <class T>
230  RES_T convolve(const Image<T> &imIn, const vector<double> &kernel, Image<T> &imOut)
231  {
232  if (&imIn==&imOut)
233  {
234  Image<T> tmpIm(imIn, true); // clone
235  return convolve(tmpIm, kernel, imOut);
236  }
237 
238  CHECK_ALLOCATED(&imIn, &imOut);
239 
240  ImageFreezer freeze(imOut);
241 
242  vertConvolve(imIn, kernel, imOut);
243  horizConvolve(imOut, kernel, imOut); // inplace safe
244 
245  return RES_OK;
246  }
247 
248 
254  template <class T>
255  RES_T gaussianFilter(const Image<T> &imIn, size_t radius, Image<T> &imOut)
256  {
257  CHECK_ALLOCATED(&imIn, &imOut);
258 
259  ImageFreezer freeze(imOut);
260 
261  int kernelSize = radius*2+1;
262  vector<double> kernel(kernelSize);
263 
264  double sigma = double(radius)/2.;
265  double sum = 0.;
266 
267  // Determine kernel coefficients
268  for (int i=0;i<kernelSize;i++)
269  {
270  kernel[i] = double(exp( -pow((double(i)-radius)/sigma, 2)/2. ));
271  sum += kernel[i];
272  }
273  // Normalize
274  for (int i=0;i<kernelSize;i++)
275  kernel[i] /= sum;
276 
277 
278  return convolve(imIn, kernel, imOut);
279 
280  }
281 
282 
285 } // namespace smil
286 
287 
288 #endif // _D_IMAGE_DRAW_HPP
289 
Definition: DColorConvert.h:38
sliceType getLines() const
Get an array containing the start offset of each line.
Definition: DImage.hpp:114
Definition: DBaseImage.h:235
RES_T convolve(const Image< T > &imIn, const vector< double > &kernel, Image< T > &imOut)
Convolution in both x and y directions using the same 1D kernel.
Definition: DImageConvolution.hpp:230
RES_T gaussianFilter(const Image< T > &imIn, size_t radius, Image< T > &imOut)
2D Gaussian filter
Definition: DImageConvolution.hpp:255
size_t getHeight() const
Get image height.
Definition: DBaseImage.h:91
RES_T horizConvolve(const Image< T > &imIn, const vector< double > &kernel, Image< T > &imOut)
Horizontal convolution.
Definition: DImageConvolution.hpp:62
size_t getWidth() const
Get image width.
Definition: DBaseImage.h:87
Main Image class.
Definition: DQVtkViewer.hpp:44
size_t getLineCount() const
Get the number of lines.
Definition: DBaseImage.h:154
Definition: DBufferPool.hpp:44
RES_T vertConvolve(const Image< T > &imIn, const vector< double > &kernel, Image< T > &imOut)
Vertical convolution.
Definition: DImageConvolution.hpp:146
volType getSlices() const
Get an array containing the start offset of each slice.
Definition: DImage.hpp:118
size_t getDepth() const
Get image depth (Z)
Definition: DBaseImage.h:95