Basic Image AlgorithmS Library  2.8.0
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
clfRecursiveGauss.cpp
1 /*
2  * clfRecursiveGauss.cpp
3  *
4  * Created on: Oct 28, 2011
5  * Author: fkellner
6  */
7 
8 #include <OpenCLFramework/Algorithm/clfRecursiveGauss.hh>
9 #include <Base/Common/BIASpragma.hh>
10 #include <Base/Image/ImageIO.hh>
11 
12 using namespace BIAS;
13 using namespace BIAS;
14 
15 clfRecursiveGauss::clfRecursiveGauss(float sigma, clfContext *ctx, bool sharedGL, unsigned int device) :
16  clfAlgorithm(ctx,sharedGL,device) {
17  intermediateImage_ = NULL;
18  oneDirFilteredImage_ = NULL;
19  packedGreyImage_ = NULL;
20 
21  // add code to our program
22  programCL_->AddSource("/OpenCLFramework/Algorithm/cl/recursivegauss.cl");
23  programCL_->Build();
24 
25  kernelNames_.clear();
26  kernelNames_.push_back("RecursiveGaussianForwardX");
27  kernelNames_.push_back("RecursiveGaussianForwardY");
28  kernelNames_.push_back("RecursiveGaussianForwardXPacked");
29  kernelNames_.push_back("RecursiveGaussianForwardYPacked");
30  kernelNames_.push_back("RecursiveGaussianBackwardX");
31  kernelNames_.push_back("RecursiveGaussianBackwardY");
32  kernelNames_.push_back("RecursiveGaussianBackwardXPacked");
33  kernelNames_.push_back("RecursiveGaussianBackwardYPacked");
34 
35  programCL_->AddKernels(kernelNames_);
36  SetSigma(sigma);
37 
38  programCL_->AddKernel("RowPackGreyToRGBA");
39  programCL_->AddKernel("ColPackGreyToRGBA");
40  programCL_->AddKernel("RowUnPackRGBAToGrey");
41  programCL_->AddKernel("ColUnPackRGBAToGrey");
42  programCL_->AddKernel("ColPackRGBAToRowPackRGBA");
43  programCL_->AddKernel("RowPackRGBAToColPackRGBA");
44 }
45 
46 clfRecursiveGauss::~clfRecursiveGauss() {
47  Reset();
48 }
49 
50 void clfRecursiveGauss::Reset() {
51  delete intermediateImage_;
52  intermediateImage_ = NULL;
53  delete oneDirFilteredImage_;
54  oneDirFilteredImage_ = NULL;
55  delete packedGreyImage_;
56  packedGreyImage_ = NULL;
57 }
58 
59 
60 int clfRecursiveGauss::Filter(clfImage2D *src, clfImage2D *&dst) {
61 
62  unsigned int w,h;
63  src->GetImageDim(w,h);
64 
65  if (w % 16 == 0 && h % 16 == 0) {
66  FilterX(src, oneDirFilteredImage_);
67  FilterY(oneDirFilteredImage_, dst);
68  } else {
69  BIASERR("recursive gauss needs images divisible by 16");
70  BIASABORT;
71  }
72  return 0;
73 }
74 
75 int clfRecursiveGauss::FilterX(clfImage2D *src, clfImage2D *&dst) {
76  programCL_->KernelSetArgument("RecursiveGaussianForwardX", 0, *src);
77  programCL_->KernelSetArgument("RecursiveGaussianBackwardX", 0, *src);
78 
79  unsigned int w,h;
80  src->GetImageDim(w,h);
81 
82  if (dst == NULL) {
83  dst = context_->CreateImage2D();
84  dst->AllocateFromTemplate(*src);
85  }
86  programCL_->KernelSetArgument("RecursiveGaussianBackwardX", 2, *dst);
87 
88  if (intermediateImage_ == NULL) {
89  intermediateImage_ = context_->CreateImage2D();
90  intermediateImage_->AllocateFromTemplate(*src);
91  }
92  programCL_->KernelSetArgument("RecursiveGaussianForwardX", 1, *intermediateImage_);
93  programCL_->KernelSetArgument("RecursiveGaussianBackwardX", 1, *intermediateImage_);
94 
95  context_->RunOn1DRange(*programCL_, "RecursiveGaussianForwardX", h, 16);
96  context_->RunOn1DRange(*programCL_, "RecursiveGaussianBackwardX", h, 16);
97  return 0;
98 }
99 
100 int clfRecursiveGauss::FilterY(clfImage2D *src, clfImage2D *&dst) {
101 
102  programCL_->KernelSetArgument("RecursiveGaussianForwardY", 0, *src);
103  programCL_->KernelSetArgument("RecursiveGaussianBackwardY", 0, *src);
104 
105  unsigned int w,h;
106  src->GetImageDim(w,h);
107 
108  if (dst == NULL) {
109  dst = context_->CreateImage2D();
110  dst->AllocateFromTemplate(*src);
111  }
112  programCL_->KernelSetArgument("RecursiveGaussianBackwardY", 2, *dst);
113 
114  if (intermediateImage_ == NULL) {
115  intermediateImage_ = context_->CreateImage2D();
116  intermediateImage_->AllocateFromTemplate(*src);
117  }
118  programCL_->KernelSetArgument("RecursiveGaussianForwardY", 1, *intermediateImage_);
119  programCL_->KernelSetArgument("RecursiveGaussianBackwardY", 1, *intermediateImage_);
120 
121  context_->RunOn1DRange(*programCL_, "RecursiveGaussianForwardY", w, 16);
122  context_->RunOn1DRange(*programCL_, "RecursiveGaussianBackwardY", w, 16);
123  return 0;
124 }
125 
126 
127 
128 int clfRecursiveGauss::FilterGrey(clfImage2D *src, clfImage2D *&dst) {
129 
130  { // Allocate temporary images
131  if (packedGreyImage_ == NULL) {
132  packedGreyImage_ = context_->CreateImage2D();
133  packedGreyImage_->AllocateFromTemplate(*src);
134  }
135  if (intermediateImage_ == NULL) {
136  intermediateImage_ = context_->CreateImage2D();
137  intermediateImage_->AllocateFromTemplate(*src);
138  }
139  if (oneDirFilteredImage_ == NULL) {
140  oneDirFilteredImage_ = context_->CreateImage2D();
141  oneDirFilteredImage_->AllocateFromTemplate(*src);
142  }
143  if (dst == NULL) {
144  dst = context_->CreateImage2D();
145  dst->AllocateFromTemplate(*src);
146  } else {
147  if (dst->IsReadOnly()) {
148  THROW_CL_EXCEPTION(cl::Error(-100, "Need to use dst as tempory read/write"));
149  }
150  }
151  }
152  unsigned int w,h;
153  { // make sure everything is divisible by 16*4
154  src->GetImageDim(w,h);
155  if (h % 64 != 0) {
156  h = h+(64-(h%64));
157  }
158  if (w % 64 != 0) {
159  w = w+(64-(w%64));
160  }
161  }
162  { // pack grey image from pixels R000R000R000 to RGBARGBARGBA
163  // rgbargbargbargba rgba
164  // R000R000R000R000 RRRR
165  // R000R000R000R000 -> RRRR
166  // R000R000R000R000 RRRR
167  programCL_->KernelSetArgument("ColPackGreyToRGBA", 0, *src);
168  programCL_->KernelSetArgument("ColPackGreyToRGBA", 1, *packedGreyImage_);
169  context_->RunOn2DRange(*programCL_, "ColPackGreyToRGBA", w/4, h, 16, 16);
170  }
171 
172  { // filter packed image, needing only w/4 read/write operations per row
173  programCL_->KernelSetArgument("RecursiveGaussianForwardXPacked", 0, *packedGreyImage_);
174  programCL_->KernelSetArgument("RecursiveGaussianBackwardXPacked", 0, *packedGreyImage_);
175 
176  programCL_->KernelSetArgument("RecursiveGaussianForwardXPacked", 1, *intermediateImage_);
177  programCL_->KernelSetArgument("RecursiveGaussianBackwardXPacked", 1, *intermediateImage_);
178  programCL_->KernelSetArgument("RecursiveGaussianBackwardXPacked", 2, *dst);
179 
180  context_->RunOn1DRange(*programCL_, "RecursiveGaussianForwardXPacked", h, 16);
181  context_->RunOn1DRange(*programCL_, "RecursiveGaussianBackwardXPacked", h, 16);
182  }
183 
184  { // "repack" image such that the grey image is now stored (RGBA)^T
185  // rgba rgba
186  // RGBA RRRR
187  // RGBA -> GGGG
188  // RGBA BBBB
189  // RGBA AAAA
190  programCL_->KernelSetArgument("ColPackRGBAToRowPackRGBA", 0, *dst);
191  programCL_->KernelSetArgument("ColPackRGBAToRowPackRGBA", 1, *oneDirFilteredImage_);
192  context_->RunOn2DRange(*programCL_, "ColPackRGBAToRowPackRGBA", w/4, h/4, 16, 16);
193  }
194 
195  { // filter packed image, needing only h/4 read/write operations per column
196  programCL_->KernelSetArgument("RecursiveGaussianForwardYPacked", 0, *oneDirFilteredImage_);
197  programCL_->KernelSetArgument("RecursiveGaussianBackwardYPacked", 0, *oneDirFilteredImage_);
198 
199  programCL_->KernelSetArgument("RecursiveGaussianForwardYPacked", 1, *intermediateImage_);
200  programCL_->KernelSetArgument("RecursiveGaussianBackwardYPacked", 1, *intermediateImage_);
201 
202  programCL_->KernelSetArgument("RecursiveGaussianBackwardYPacked", 2, *packedGreyImage_);
203 
204  context_->RunOn1DRange(*programCL_, "RecursiveGaussianForwardYPacked", w, 16);
205  context_->RunOn1DRange(*programCL_, "RecursiveGaussianBackwardYPacked", w, 16);
206  }
207 
208  { // unpack to grey image again
209  programCL_->KernelSetArgument("RowUnPackRGBAToGrey", 0, *packedGreyImage_);
210  programCL_->KernelSetArgument("RowUnPackRGBAToGrey", 1, *dst);
211  context_->RunOn2DRange(*programCL_, "RowUnPackRGBAToGrey", w, h/4, 16, 16);
212  }
213  return 0;
214 }
215 
216 int clfRecursiveGauss::FilterGreyColPacked(clfImage2D *src, clfImage2D *&dst) {
217 // static clfImage2D *temp1 = context_->CreateImage2D();
218 // static clfImage2D *temp2 = context_->CreateImage2D();
219 // temp1->Allocate(src->StorageType(), src->ColorModel(), src->Width()*4, src->Height(), false, false);
220 // temp2->Allocate(src->StorageType(), src->ColorModel(), src->Width()*4, src->Height(), false, false);
221 // programCL_->KernelSetArgument("ColUnPackRGBAToGrey", 0, *src);
222 // programCL_->KernelSetArgument("ColUnPackRGBAToGrey", 1, *temp1);
223 // context_->RunOn2DRange(*programCL_, "ColUnPackRGBAToGrey", src->Width(), src->Height(), 16, 16);
224 // context_->Finish();
225 // FilterGrey(temp1, temp2);
226 // programCL_->KernelSetArgument("ColPackGreyToRGBA", 0, *temp2);
227 // programCL_->KernelSetArgument("ColPackGreyToRGBA", 1, *dst);
228 // context_->RunOn2DRange(*programCL_, "ColPackGreyToRGBA", src->Width(), src->Height(), 16, 16);
229 // return 0;
230 //
231  unsigned int w,h, worg, horg;
232  src->GetImageDim(worg,horg);
233  src->GetImageDim(w,h);
234  w *= 4;
235 
236  if (h % 64 != 0) {
237  h = h+(64-(h%64));
238  }
239  if (w % 64 != 0) {
240  w = w+(64-(w%64));
241  }
242 
243  { // Allocate temporary images
244  if (packedGreyImage_ == NULL) {
245  packedGreyImage_ = context_->CreateImage2D();
246  packedGreyImage_->Allocate(src->StorageType(), src->ColorModel(), w, h, 0, false, false);
247  }
248  if (intermediateImage_ == NULL) {
249  intermediateImage_ = context_->CreateImage2D();
250  intermediateImage_->AllocateFromTemplate(*packedGreyImage_);
251  }
252  if (oneDirFilteredImage_ == NULL) {
253  oneDirFilteredImage_ = context_->CreateImage2D();
254  oneDirFilteredImage_->AllocateFromTemplate(*packedGreyImage_);
255  }
256  if (dst == NULL) {
257  dst = context_->CreateImage2D();
258  dst->AllocateFromTemplate(*src);
259  } else {
260  if (dst->IsReadOnly()) {
261  THROW_CL_EXCEPTION(cl::Error(-100, "Need to use dst as temporary read/write"));
262  }
263  }
264  }
265 
266  { // filter packed image, needing only w/4 read/write operations per row
267  programCL_->KernelSetArgument("RecursiveGaussianForwardXPacked", 0, *src);
268  programCL_->KernelSetArgument("RecursiveGaussianBackwardXPacked", 0, *src);
269 
270  programCL_->KernelSetArgument("RecursiveGaussianForwardXPacked", 1, *intermediateImage_);
271  programCL_->KernelSetArgument("RecursiveGaussianBackwardXPacked", 1, *intermediateImage_);
272  programCL_->KernelSetArgument("RecursiveGaussianBackwardXPacked", 2, *packedGreyImage_);
273 
274  context_->RunOn1DRange(*programCL_, "RecursiveGaussianForwardXPacked", h, 16);
275  context_->RunOn1DRange(*programCL_, "RecursiveGaussianBackwardXPacked", h, 16);
276  context_->Finish();
277  }
278 
279  { // "repack" image such that the grey image is now stored (RGBA)^T
280  // rgba rgba
281  // RGBA RRRR
282  // RGBA -> GGGG
283  // RGBA BBBB
284  // RGBA AAAA
285  programCL_->KernelSetArgument("ColPackRGBAToRowPackRGBA", 0, *packedGreyImage_);
286  programCL_->KernelSetArgument("ColPackRGBAToRowPackRGBA", 1, *oneDirFilteredImage_);
287  context_->RunOn2DRange(*programCL_, "ColPackRGBAToRowPackRGBA", w/4, h/4, 16, 16);
288  context_->Finish();
289  }
290 
291  { // filter packed image, needing only h/4 read/write operations per column
292  programCL_->KernelSetArgument("RecursiveGaussianForwardYPacked", 0, *oneDirFilteredImage_);
293  programCL_->KernelSetArgument("RecursiveGaussianBackwardYPacked", 0, *oneDirFilteredImage_);
294 
295  programCL_->KernelSetArgument("RecursiveGaussianForwardYPacked", 1, *intermediateImage_);
296  programCL_->KernelSetArgument("RecursiveGaussianBackwardYPacked", 1, *intermediateImage_);
297 
298  programCL_->KernelSetArgument("RecursiveGaussianBackwardYPacked", 2, *packedGreyImage_);
299 
300  context_->RunOn1DRange(*programCL_, "RecursiveGaussianForwardYPacked", w, 16);
301  context_->RunOn1DRange(*programCL_, "RecursiveGaussianBackwardYPacked", w, 16);
302  }
303  { // unpack to grey image again
304  programCL_->KernelSetArgument("RowPackRGBAToColPackRGBA", 0, *packedGreyImage_);
305  programCL_->KernelSetArgument("RowPackRGBAToColPackRGBA", 1, *intermediateImage_);
306  context_->RunOn2DRange(*programCL_, "RowPackRGBAToColPackRGBA", w/4, h/4, 16, 16);
307  }
308  // copy to dst
309  intermediateImage_->CopyToImage(*dst, 0,0, 0,0, worg, horg);
310  return 0;
311 }
312 
313 
314 void clfRecursiveGauss::SetSigma(float sigma, int order) {
315 // float &a0=params_[0], float &a1, float &a2, float &a3=params_[3],
316 // float &b1=params_[4], float &b1, float &coefp=params_[6], float &coefn=params_[7]);
317  // pre-compute filter coefficients
318  float alpha = 1.695f / sigma;
319  float ema = exp(-alpha);
320  float ema2 = exp(-2.0f * alpha);
321  params_.resize(8);
322  params_[0] = 0.0f;
323  params_[1] = 0.0f;
324  params_[2] = 0.0f;
325  params_[3] = 0.0f;
326  params_[4] = -2.0f * ema;
327  params_[5] = ema2;
328  params_[6] = 0.0f;
329  params_[7] = 0.0f;
330  switch (order) {
331  case 0: {
332  const float k = (1.0f - ema)*(1.0f - ema)/(1.0f + (2.0f * alpha * ema) - ema2);
333  params_[0] = k;
334  params_[1] = k * (alpha - 1.0f) * ema;
335  params_[2] = k * (alpha + 1.0f) * ema;
336  params_[3] = -k * ema2;
337  }
338  break;
339  case 1: {
340  params_[0] = (1.0f - ema) * (1.0f - ema);
341  params_[1] = 0.0f;
342  params_[2] = -params_[0];
343  params_[3] = 0.0f;
344  }
345  break;
346  case 2: {
347  const float ea = exp(-alpha);
348  const float k = -(ema2 - 1.0f)/(2.0f * alpha * ema);
349  float kn = -2.0f * (-1.0f + (3.0f * ea) - (3.0f * ea * ea) + (ea * ea * ea));
350  kn /= (((3.0f * ea) + 1.0f + (3.0f * ea * ea) + (ea * ea * ea)));
351  params_[0] = kn;
352  params_[1] = -kn * (1.0f + (k * alpha)) * ema;
353  params_[2] = kn * (1.0f - (k * alpha)) * ema;
354  params_[3] = -kn * ema2;
355  }
356  break;
357  default:
358  break;
359  }
360  params_[6] = (params_[0] + params_[1])/(1.0f + params_[4] + params_[5]);
361  params_[7] = (params_[2] + params_[3])/(1.0f + params_[4] + params_[5]);
362 
363  for (unsigned int i=0;i<kernelNames_.size();i++) {
364  for (unsigned int j=0;j<params_.size();j++) {
365  if (i<4)
366  programCL_->KernelSetArgument(kernelNames_[i], j+2, params_[j]);
367  else
368  programCL_->KernelSetArgument(kernelNames_[i], j+3, params_[j]);
369  }
370  }
371 }
372 
BIAS::clfRecursiveGauss::FilterX
int FilterX(clfImage2D *src, clfImage2D *&dst)
filter in X-direction
Definition: clfRecursiveGauss.cpp:75
BIAS::clfAlgorithm::programCL_
clfProgram * programCL_
Definition: clfAlgorithm.hh:35
BIAS::clfProgram::AddSource
std::string AddSource(std::string filename)
adds source code from a file
Definition: clfProgram.cpp:64
BIAS::clfImage2D::GetImageDim
void GetImageDim(unsigned int &width, unsigned int &height) const
Definition: clfImage2D.hh:110
BIAS::clfImage2D::Allocate
void Allocate(BIAS::ImageBase::EStorageType st, BIAS::ImageBase::EColorModel cm, unsigned int width, unsigned int height, unsigned int stride=0, bool readonly=false, bool writeonly=false, const void *hostptr=NULL, bool copy=false)
Allocation of a memory buffer as 2D image, either call directly or use wrapper for BIAS::ImageBase...
Definition: clfImage2D.cpp:62
BIAS::clfImage2D
OpenCL Image2D wrapper.
Definition: clfImage2D.hh:46
BIAS::clfRecursiveGauss::params_
std::vector< float > params_
Definition: clfRecursiveGauss.hh:72
BIAS::clfImage2D::ColorModel
BIAS::ImageBase::EColorModel ColorModel() const
Definition: clfImage2D.hh:76
BIAS::clfProgram::Build
void Build(int deviceNr=0, std::string options="")
builds the sources added by AddSource and AddSourceFromString
Definition: clfProgram.cpp:115
BIAS::clfImage2D::StorageType
BIAS::ImageBase::EStorageType StorageType() const
Definition: clfImage2D.hh:77
BIAS::clfContext
OpenCL Context wrapper.
Definition: clfContext.hh:49
BIAS::clfContext::Finish
void Finish()
force finishing the command queue
Definition: clfContext.cpp:289
BIAS::clfRecursiveGauss::Filter
int Filter(clfImage2D *src, clfImage2D *&dst)
gauss filter rgba image
Definition: clfRecursiveGauss.cpp:60
BIAS::clfRecursiveGauss::clfRecursiveGauss
clfRecursiveGauss(float sigma=2.0f, clfContext *ctx=NULL, bool sharedGL=false, unsigned int device=0)
constructor, give context or let
Definition: clfRecursiveGauss.cpp:15
BIAS::clfRecursiveGauss::FilterGrey
int FilterGrey(clfImage2D *src, clfImage2D *&dst)
optimized gauss filter for grey images
Definition: clfRecursiveGauss.cpp:128
BIAS::clfProgram::KernelSetArgument
void KernelSetArgument(std::string kernelname, unsigned int argnumber, clfBuffer &buffer)
set kernel argument
Definition: clfProgram.cpp:177
BIAS::clfMemory::IsReadOnly
bool IsReadOnly()
Definition: clfMemory.cpp:62
BIAS::clfProgram::AddKernel
void AddKernel(std::string kernelname)
adds a kernel to the program.
Definition: clfProgram.cpp:158
BIAS::clfContext::RunOn2DRange
void RunOn2DRange(clfProgram &program, std::string kernelname, unsigned int globalrangeX, unsigned int globalrangeY, unsigned int localrangeX=0, unsigned int localrangeY=0)
run a kernel on a 2D memory range
Definition: clfContext.cpp:234
BIAS::clfRecursiveGauss::Reset
void Reset()
clear intermediate images resets storage of internal images, needs to be called if image sizes differ...
Definition: clfRecursiveGauss.cpp:50
BIAS::clfRecursiveGauss::FilterGreyColPacked
int FilterGreyColPacked(clfImage2D *src, clfImage2D *&dst)
optimized gauss filter for already packed grey images
Definition: clfRecursiveGauss.cpp:216
BIAS::clfProgram::AddKernels
void AddKernels(std::vector< std::string > kernelnames)
adds a vector of kernels to the program.
Definition: clfProgram.cpp:171
BIAS::clfImage2D::AllocateFromTemplate
void AllocateFromTemplate(const clfImage2D &src, bool readonly=false, bool writeonly=false)
Definition: clfImage2D.cpp:142
BIAS::clfAlgorithm::context_
clfContext * context_
Definition: clfAlgorithm.hh:34
BIAS::clfRecursiveGauss::kernelNames_
std::vector< std::string > kernelNames_
Definition: clfRecursiveGauss.hh:74
BIAS::clfContext::CreateImage2D
clfImage2D * CreateImage2D()
create buffer object
Definition: clfContext.hh:87
BIAS::clfContext::RunOn1DRange
void RunOn1DRange(clfProgram &program, std::string kernelname, unsigned int globalrange, unsigned int localrange=0)
run a kernel on a 1D memory range
Definition: clfContext.cpp:220
BIAS::clfRecursiveGauss::FilterY
int FilterY(clfImage2D *src, clfImage2D *&dst)
filter in Y-direction
Definition: clfRecursiveGauss.cpp:100
BIAS::clfRecursiveGauss::SetSigma
void SetSigma(float sigma, int order=0)
Set Sigma and calculate coefficients for iir filter.
Definition: clfRecursiveGauss.cpp:314
BIAS::clfImage2D::CopyToImage
void CopyToImage(clfImage2D &outputimage, unsigned int srcoriginX=0, unsigned int srcoriginY=0, unsigned int dstoriginX=0, unsigned int dstoriginY=0, unsigned int regionX=0, unsigned int regionY=0, bool blocking=true)
Definition: clfImage2D.cpp:233
Generated on Tue Feb 6 2018 02:01:42 for Basic Image AlgorithmS Library by   doxygen 1.8.5