SMIL  0.9.1
DMorphImageOperations.hxx
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 _MORPH_IMAGE_OPERATIONS_HXX
31 #define _MORPH_IMAGE_OPERATIONS_HXX
32 
33 #include "DMorphImageOperations.hpp"
34 
35 #ifdef USE_OPEN_MP
36 #include <omp.h>
37 #endif // USE_OPEN_MP
38 
39 namespace smil
40 {
41 
42  template <class T_in, class T_out>
43  RES_T MorphImageFunctionBase<T_in, T_out>::initialize(const imageInType &imIn, imageOutType &imOut, const StrElt &se)
44  {
45  this->imageIn = &imIn;
46  this->imageOut = &imOut;
47 
48  imIn.getSize(imSize);
49  slicesIn = imIn.getSlices();
50  slicesOut = imOut.getSlices();
51  pixelsIn = imIn.getPixels();
52  pixelsOut = imOut.getPixels();
53 
54  sePoints = se.points;
55 
56  oddSe = se.odd;
57 
58  sePointNbr = sePoints.size();
59  relativeOffsets.clear();
60  vector<IntPoint>::iterator pt = sePoints.begin();
61  se_xmin = ImDtTypes<int>::max();
62  se_xmax = ImDtTypes<int>::min();
63  se_ymin = ImDtTypes<int>::max();
64  se_ymax = ImDtTypes<int>::min();
65  se_zmin = ImDtTypes<int>::max();
66  se_zmax = ImDtTypes<int>::min();
67  while(pt!=sePoints.end())
68  {
69  if(pt->x < se_xmin) se_xmin = pt->x;
70  if(pt->x > se_xmax) se_xmax = pt->x;
71  if(pt->y < se_ymin) se_ymin = pt->y;
72  if(pt->y > se_ymax) se_ymax = pt->y;
73  if(pt->z < se_zmin) se_zmin = pt->z;
74  if(pt->z > se_zmax) se_zmax = pt->z;
75 
76  relativeOffsets.push_back(pt->x + pt->y*imSize[0] + pt->z*imSize[0]*imSize[1]);
77  pt++;
78  }
79  return RES_OK;
80  }
81 
82 
83  template <class T_in, class T_out>
84  RES_T MorphImageFunctionBase<T_in, T_out>::finalize(const imageInType & /*imIn*/, imageOutType & /*imOut*/, const StrElt & /*se*/)
85  {
86  return RES_OK;
87  }
88 
89  template <class T_in, class T_out>
90  RES_T MorphImageFunctionBase<T_in, T_out>::_exec(const imageInType &imIn, imageOutType &imOut, const StrElt &se)
91  {
92  ASSERT_ALLOCATED(&imIn)
93  ASSERT_SAME_SIZE(&imIn, &imOut)
94 
95  if ((void*)&imIn==(void*)&imOut)
96  {
97  Image<T_in> tmpIm(imIn, true); // clone
98  return _exec(imIn, imOut, se);
99  }
100 
101  ImageFreezer freeze(imOut);
102 
103  StrElt se2;
104  if (se.size>1)
105  se2 = se.homothety(se.size);
106  else se2 = se;
107 
108  initialize(imIn, imOut, se2);
109 
110  RES_T retVal;
111 
112  retVal = processImage(imIn, imOut, se2);
113 
114  finalize(imIn, imOut, se);
115 
116  return retVal;
117  }
118  template <class T_in, class T_out>
119  RES_T MorphImageFunctionBase<T_in, T_out>::_exec(const imageInType &imIn, const StrElt &se)
120  {
121  ASSERT_ALLOCATED(&imIn)
122 
123  imageOutType &_imOut = (imageOutType&)imIn;
124 
125  StrElt se2;
126  if (se.size>1)
127  se2 = se.homothety(se.size);
128  else se2 = se;
129 
130  this->initialize(imIn, _imOut, se2);
131 
132  RES_T retVal;
133 
134  retVal = processImage(imIn, _imOut, se2);
135 
136  finalize(imIn, _imOut, se);
137 
138  return retVal;
139  }
140 
141  template <class T_in, class T_out>
142  RES_T MorphImageFunctionBase<T_in, T_out>::processImage(const imageInType &/*imIn*/, imageOutType &/*imOut*/, const StrElt &se)
143  {
144  for(size_t curSlice=0;curSlice<imSize[2];curSlice++)
145  {
146  processSlice(*slicesIn, *slicesOut, imSize[1], se);
147  slicesIn++;
148  slicesOut++;
149  }
150  return RES_OK;
151  }
152 // virtual RES_T processImage(imageInType &imIn, imageOutType &imOut, hSE &se)
153 // {
154 // }
155 
156  template <class T_in, class T_out>
157  void MorphImageFunctionBase<T_in, T_out>::processSlice(sliceInType linesIn, sliceOutType linesOut, size_t &lineNbr, const StrElt &se)
158  {
159  for(size_t curLine=0;curLine<lineNbr;curLine++)
160  processLine(linesIn[curLine], linesOut[curLine], imSize[0], se);
161  }
162 
163  // TODO : offset list for 3D odd SE !!
164  template <class T_in, class T_out>
165  void MorphImageFunctionBase<T_in, T_out>::processLine(lineInType pixIn, lineOutType /*pixOut*/, size_t &pixNbr, const StrElt &se)
166  {
167  int x, y, z;
168  IntPoint p;
169  size_t offset = pixIn - pixelsIn;
170  vector<IntPoint> ptList;
171  vector<int> relOffsetList;
172  vector<int> offsetList;
173 
174  size_t curSlice = offset / (imSize[0]*imSize[1]);
175  size_t curLine = (offset / imSize[0])%imSize[1];
176  size_t curPixel = 0;
177 
178  bool oddLine = se.odd && (curLine)%2;
179  // int dx;
180 
181  // Remove points wich are outside the image
182  for (UINT i=0;i<sePointNbr;i++)
183  {
184  p = sePoints[i];
185  y = curLine + p.y;
186  z = curSlice + p.z;
187  if (y>=0 && y<int(imSize[1]) && z>=0 && z<int(imSize[2]))
188  {
189  if (oddLine && ((y+1)%2)!=0)
190  p.x += 1;
191  ptList.push_back(p);
192  if (oddLine && ((y+1)%2)!=0)
193  relOffsetList.push_back(relativeOffsets[i]+1);
194  else
195  relOffsetList.push_back(relativeOffsets[i]);
196  }
197  }
198  UINT ptNbr = ptList.size();
199 
200  // Left border
201  while((int)curPixel < -se_xmin)
202  {
203  offsetList.clear();
204  for (UINT i=0;i<ptNbr;i++)
205  {
206  x = curPixel + ptList[i].x;
207 
208  if (x>=0 && x<(int)imSize[0])
209  offsetList.push_back(relOffsetList[i]);
210  }
211  processPixel(offset, offsetList);
212  curPixel++;
213  offset++;
214  }
215 
216  // Middle
217  offsetList.clear();
218  for (UINT i=0;i<ptNbr;i++)
219  offsetList.push_back(relOffsetList[i]);
220  while(curPixel < pixNbr-se_xmax)
221  {
222  processPixel(offset, offsetList);
223  curPixel++;
224  offset++;
225  }
226 
227  // Right border
228  while(curPixel<pixNbr)
229  {
230  offsetList.clear();
231  for (UINT i=0;i<ptNbr;i++)
232  {
233  x = curPixel + ptList[i].x;
234 
235  if (x>=0 && x<int(imSize[0]))
236  offsetList.push_back(relOffsetList[i]);
237  }
238  #ifdef USE_OPEN_MP
239  #pragma omp critical
240  #endif // USE_OPEN_MP
241  {
242  processPixel(offset, offsetList);
243  }
244  curPixel++;
245  offset++;
246  }
247  }
248 
249  template <class T_in, class T_out>
250  void MorphImageFunctionBase<T_in, T_out>::processPixel(size_t /*pointOffset*/, vector<int> &dOffsets)
251  {
252  // Example: dilation function
253  vector<int>::iterator it = dOffsets.begin();
254  while(it!=dOffsets.end())
255  {
256 // pixelsOut[pointOffset] = max(pixelsOut[pointOffset], pixelsIn[pointOffset + *dOffset]);
257  it++;
258  }
259  }
260 
261 
266 
267  template <class T_in, class lineFunction_T, class T_out, bool Enable>
269  {
270  int st = se.getType();
271 
272  switch(st)
273  {
274  case SE_Horiz:
275  return true;
276  case SE_Vert:
277  return true;
278  case SE_Hex:
279  return true;
280  case SE_Squ:
281  return true;
282  case SE_Cube:
283  return true;
284  case SE_Cross:
285  return true;
286  default:
287  return false;
288  }
289  }
290 
291  template <class T_in, class lineFunction_T, class T_out, bool Enable>
292  RES_T MorphImageFunction<T_in, lineFunction_T, T_out, Enable>::initialize(const imageInType &imIn, imageOutType &/*imOut*/, const StrElt &/*se*/)
293  {
294  this->lineLen = imIn.getWidth();
295 
296  this->borderBuf = ImDtTypes<T_in>::createLine(this->lineLen);
297  this->cpBuf = ImDtTypes<T_in>::createLine(this->lineLen);
298  fillLine<T_in> f;
299  f(this->borderBuf, this->lineLen, this->borderValue);
300 
301  return RES_OK;
302  }
303 
304  template <class T_in, class lineFunction_T, class T_out, bool Enable>
305  RES_T MorphImageFunction<T_in, lineFunction_T, T_out, Enable>::finalize(const imageInType &/*imIn*/, imageOutType &/*imOut*/, const StrElt &/*se*/)
306  {
307  ImDtTypes<T_in>::deleteLine(this->borderBuf);
308  ImDtTypes<T_in>::deleteLine(this->cpBuf);
309 
310  return RES_OK;
311  }
312 
313  template <class T_in, class lineFunction_T, class T_out, bool Enable>
314  RES_T MorphImageFunction<T_in, lineFunction_T, T_out, Enable>::_exec(const imageInType &imIn, imageOutType &imOut, const StrElt &se)
315  {
316  ASSERT_ALLOCATED(&imIn)
317  ASSERT_SAME_SIZE(&imIn, &imOut)
318 
319  if ((void*)&imIn==(void*)&imOut)
320  {
321  Image<T_in> tmpIm(imIn, true); // clone
322  return _exec(imIn, imOut, se);
323  }
324 
325  ImageFreezer freeze(imOut);
326 
327  StrElt se2;
328  if (se.size>1)
329  se2 = se.homothety(se.size);
330  else se2 = se;
331 
332  this->initialize(imIn, imOut, se2);
333 
334  RES_T retVal;
335 
336  retVal = _exec_single(imIn, imOut, se2);
337 
338  this->finalize(imIn, imOut, se);
339 
340  return retVal;
341  }
342 
343  template <class T_in, class lineFunction_T>
345  {
346  ASSERT_ALLOCATED(&imIn, &imOut);
347 
348  int seSize = se.size;
349  int seType = se.getType () ;
350 
351  if (seSize==0)
352  {
353  copy(imIn, imOut);
354  return RES_OK;
355  }
356 
357  this->initialize(imIn, imOut, se);
358 
359  if (seType==SE_Rhombicuboctahedron)
360  {
361  _exec_rhombicuboctahedron (imIn, imOut, se.size);
362  this->finalize(imIn, imOut, se);
363  return RES_OK;
364  }
365 
366 
367  ImageFreezer freezer(imOut);
368 
369  Image<T_in> *inImage;
370  Image<T_in> *outImage, *tmpImage = NULL;
371 
372  inImage = (Image<T_in>*)&imIn;
373 
374  if (this->isInplaceSafe(se) || (&imIn!=&imOut && seSize==1))
375  outImage = (Image<T_in> *)&imOut;
376  else
377  outImage = tmpImage = new Image<T_in>(imOut);
378 
379 
380  // else
381  {
382  _exec_single(*inImage, *outImage, se);
383 
384  if (seSize>1)
385  {
386  if (tmpImage)
387  inImage = tmpImage;
388  else inImage = &imOut;
389 
390  outImage = &imOut;
391 
392  for (int i=1;i<seSize;i++)
393  {
394  _exec_single(*inImage, *outImage, se);
395  if (i<seSize-1)
396  swap(inImage, outImage);
397  }
398  }
399  }
400 
401  if (tmpImage)
402  {
403  if (outImage!=&imOut)
404  copy(*outImage, imOut);
405  delete tmpImage;
406  }
407 
408  this->finalize(imIn, imOut, se);
409 
410  imOut.modified();
411  return RES_OK;
412  }
413 
414 
415 
416  template <class T_in, class lineFunction_T>
418  {
419  int st = se.getType();
420 
421  switch(st)
422  {
423  case SE_Hex:
424  return _exec_single_hexagonal_SE(imIn, imOut);
425  case SE_Squ:
426  return _exec_single_square_SE(imIn, imOut);
427  case SE_Horiz:
428  return _exec_single_horizontal_segment(imIn, 1, imOut);
429  case SE_Cross:
430  return _exec_single_cross(imIn, imOut);
431  case SE_Vert:
432  return _exec_single_vertical_segment(imIn, imOut);
433  case SE_Cube:
434  return _exec_single_cube_SE(imIn, imOut);
435  // case SE_Cross3D:
436  // return _exec_single_cross_3d(imIn, imOut);
437  default:
438  return _exec_single_generic(imIn, imOut, se);
439  }
440 
441  return RES_ERR_NOT_IMPLEMENTED;
442  }
443 
444 
445 
446  template <class T_in, class lineFunction_T>
448  {
449  int sePtsNumber = se.points.size();
450  if (sePtsNumber==0)
451  return RES_OK;
452 
453  int nSlices = imIn.getSliceCount();
454  int nLines = imIn.getHeight();
455 
456  int nthreads = Core::getInstance()->getNumberOfThreads();
457  lineType *_bufs = this->createAlignedBuffers(2*nthreads, this->lineLen);
458  lineType tmpBuf = _bufs[0];
459  lineType tmpBuf2 = _bufs[nthreads];
460 
461  const Image<T_in> *tmpIm;
462 
463  if (&imIn==&imOut)
464  tmpIm = new Image<T_in>(imIn, true); // clone
465  else tmpIm = &imIn;
466 
467  //volType srcSlices = tmpIm->getSlices();
468  volType destSlices = imOut.getSlices();
469 
470  //lineInType *srcLines;
471  lineType *destLines, lineOut;
472 
473  bool oddSe = se.odd;
474  int oddLine = 0;
475 
476  int l, p;
477  #ifdef USE_OPEN_MP
478  int tid;
479  #endif // USE_OPEN_MP
480  int x, y, z;
481  vector<IntPoint> pts = se.points;
482 
483 
484  for (int s=0;s<nSlices;s++)
485  {
486  destLines = destSlices[s];
487 
488  #ifdef USE_OPEN_MP
489  #pragma omp parallel private(tid,tmpBuf,tmpBuf2,x,y,z,lineOut,p) firstprivate(pts,oddLine) num_threads(nthreads)
490  #endif // USE_OPEN_MP
491  {
492  #ifdef USE_OPEN_MP
493  tid = omp_get_thread_num();
494  tmpBuf = _bufs[tid];
495  tmpBuf2 = _bufs[tid+nthreads];
496  #endif // _OPENMP
497 
498 
499  #ifdef USE_OPEN_MP
500  #pragma omp for
501  #endif // USE_OPEN_MP
502  for (l=0;l<nLines;l++)
503  {
504  if (oddSe)
505  oddLine = ((l+1)%2 && (s+1)%2);
506  z = s - pts[0].z;
507  y = l - pts[0].y;
508  x = pts[0].x + (oddLine && y%2);
509 
510  this->_extract_translated_line(tmpIm, x, y, z, tmpBuf);
511 
512  lineOut = destLines[l];
513  for (p=1;p<sePtsNumber;p++)
514  {
515  z = s - pts[p].z;
516  y = l - pts[p].y;
517  x = pts[p].x + (oddLine && y%2);
518 
519  this->_extract_translated_line(tmpIm, x, y, z, tmpBuf2);
520  this->lineFunction._exec(tmpBuf, tmpBuf2, this->lineLen, tmpBuf);
521  }
522 
523  copyLine<T_in>(tmpBuf, this->lineLen, lineOut);
524  }
525  }
526  }
527 
528  if (&imIn==&imOut)
529  delete tmpIm;
530 
531  return RES_OK;
532  }
533 
534 
535  template <class T_in, class lineFunction_T>
537  {
538  int nSlices = imIn.getSliceCount();
539  int nLines = imIn.getHeight();
540 
541  // int nthreads = Core::getInstance()->getNumberOfThreads();
542  typename ImDtTypes<T_in>::matrixType _bufs(5, typename ImDtTypes<T_in>::vectorType(this->lineLen));
543 
544  lineType buf0 = _bufs[0].data();
545  lineType buf1 = _bufs[1].data();
546  lineType buf2 = _bufs[2].data();
547  lineType buf3 = _bufs[3].data();
548  lineType buf4 = _bufs[4].data();
549 
550  lineType tmpBuf;
551 
552  volType srcSlices = imIn.getSlices();
553  volType destSlices = imOut.getSlices();
554 
555  sliceType srcLines;
556  sliceType destLines;
557 
558  lineType curSrcLine;
559  lineType curDestLine;
560 
561  for (int s=0;s<nSlices;s++)
562  {
563  srcLines = srcSlices[s];
564  destLines = destSlices[s];
565 
566  // Process first line
567  this->_exec_shifted_line(srcLines[0], srcLines[0], -1, this->lineLen, buf0, this->cpBuf);
568  this->_exec_shifted_line(buf0, buf0, 1, this->lineLen, buf3, buf4);
569 
570  this->_exec_shifted_line(srcLines[1], srcLines[1], 1, this->lineLen, buf1, this->cpBuf);
571  lineFunction(buf3, buf1, this->lineLen, buf4);
572  lineFunction(this->borderBuf, buf4, this->lineLen, destLines[0]);
573 
574  // int tid;
575  // #pragma omp parallel
576  {
577  int l;
578 
579  // #pragma omp parallel for private(l) shared(tmpIm) ordered
580  for (l=2;l<nLines;l++)
581  {
582 
583  curSrcLine = srcLines[l];
584  curDestLine = destLines[l-1];
585 
586  if(!((l%2==0) ^ (s%2==0)))
587  {
588  this->_exec_shifted_line(curSrcLine, curSrcLine, -1, this->lineLen, buf2, this->cpBuf);
589  this->_exec_shifted_line(buf1, buf1, -1, this->lineLen, buf3, buf4);
590  }
591  else
592  {
593  this->_exec_shifted_line(curSrcLine, curSrcLine, 1, this->lineLen, buf2, this->cpBuf);
594  this->_exec_shifted_line(buf1, buf1, 1, this->lineLen, buf3, buf4);
595  }
596 
597  lineFunction(buf0, buf2, this->lineLen, buf4);
598  lineFunction(buf3, buf4, this->lineLen, curDestLine);
599  tmpBuf = buf0;
600  buf0 = buf1;
601  buf1 = buf2;
602  buf2 = tmpBuf;
603  }
604  }
605 
606  if (!((nLines%2==0) ^ (s%2==0)))
607  this->_exec_shifted_line(buf1, buf1, -1, this->lineLen, buf3, buf4);
608  else
609  this->_exec_shifted_line(buf1, buf1, 1, this->lineLen, buf3, buf4);
610  lineFunction(buf3, buf0, this->lineLen, buf4);
611  lineFunction(this->borderBuf, buf4, this->lineLen, destLines[nLines-1]);
612 
613  }
614 
615  // this->deleteAlignedBuffers();
616 
617  return RES_OK;
618  }
619 
620  template <class T_in, class lineFunction_T>
622  {
623  _exec_single_vertical_segment(imIn, imOut);
624  _exec_single_horizontal_segment(imOut, 1, imOut);
625 
626  return RES_OK;
627  }
628 
629  template <class T_in, class lineFunction_T>
631  {
632  _exec_single_vertical_segment(imIn, imOut);
633  _exec_single_horizontal_segment(imOut, 1, imOut);
634  _exec_single_depth_segment(imOut, 1, imOut);
635 
636  return RES_OK;
637  }
638 
639  template <class T_in, class lineFunction_T>
641  {
642  int lineCount = imIn.getLineCount();
643 
644  int nthreads = Core::getInstance()->getNumberOfThreads();
645  lineType *_bufs = this->createAlignedBuffers(nthreads, this->lineLen);
646  lineType buf = _bufs[0];
647 
648  sliceType srcLines = imIn.getLines();
649  sliceType destLines = imOut.getLines();
650 
651  #ifdef USE_OPEN_MP
652  int tid;
653  #endif // USE_OPEN_MP
654  int l;
655 
656  #ifdef USE_OPEN_MP
657  #pragma omp parallel private(tid, buf) num_threads(nthreads)
658  #endif // USE_OPEN_MP
659  {
660  #ifdef USE_OPEN_MP
661  tid = omp_get_thread_num();
662  buf = _bufs[tid];
663  #pragma omp for
664  #endif
665  for (l=0;l<lineCount;l++)
666  {
667  // Todo: if oddLines...
668  shiftLine<T_in>(srcLines[l], dx, this->lineLen, buf, this->borderValue);
669  this->lineFunction(buf, srcLines[l], this->lineLen, destLines[l]);
670  }
671  }
672  return RES_OK;
673  }
674 
675  template <class T_in, class lineFunction_T>
677  {
678  int imHeight = imIn.getHeight();
679  volType srcSlices = imIn.getSlices();
680  volType destSlices = imOut.getSlices();
681  sliceType srcLines;
682  sliceType destLines;
683 
684  int l;
685 
686  for (size_t s=0;s<imIn.getDepth();s++)
687  {
688  srcLines = srcSlices[s];
689  destLines = destSlices[s];
690 
691  if (dy>0)
692  {
693  for (l=0;l<imHeight-dy;l++)
694  this->lineFunction(srcLines[l], srcLines[l+dy], this->lineLen, destLines[l]);
695  for (l=imHeight-dy;l<imHeight;l++)
696  this->lineFunction(srcLines[l], this->borderBuf, this->lineLen, destLines[l]);
697  }
698  else
699  {
700  for (l=imHeight-1;l>=-dy;l--)
701  this->lineFunction(srcLines[l], srcLines[l+dy], this->lineLen, destLines[l]);
702  for (l=-dy-1;l>=0;l--)
703  this->lineFunction(srcLines[l], this->borderBuf, this->lineLen, destLines[l]);
704  }
705  }
706  return RES_OK;
707  }
708 
709  template <class T_in, class lineFunction_T>
711  {
712  int lineCount = imIn.getLineCount();
713 
714  int nthreads = Core::getInstance()->getNumberOfThreads();
715  lineType *_bufs = this->createAlignedBuffers(2*nthreads, this->lineLen);
716  lineType buf1 = _bufs[0];
717  lineType buf2 = _bufs[nthreads];
718 
719  sliceType srcLines = imIn.getLines();
720  sliceType destLines = imOut.getLines();
721 
722  lineType lineIn;
723 
724  #ifdef USE_OPEN_MP
725  int tid;
726  #endif // USE_OPEN_MP
727  int l, dx = xsize;
728 
729  #ifdef USE_OPEN_MP
730  #pragma omp parallel private(tid,buf1,buf2,lineIn) firstprivate(dx) num_threads(nthreads)
731  #endif // USE_OPEN_MP
732  {
733  #ifdef USE_OPEN_MP
734  tid = omp_get_thread_num();
735  buf1 = _bufs[tid];
736  buf2 = _bufs[tid+nthreads];
737  #pragma omp for
738  #endif
739  for (l=0;l<lineCount;l++)
740  {
741  // Todo: if oddLines...
742  lineIn = srcLines[l];
743  shiftLine<T_in>(lineIn, dx, this->lineLen, buf1, this->borderValue);
744  this->lineFunction(buf1, lineIn, this->lineLen, buf2);
745  shiftLine<T_in>(lineIn, -dx, this->lineLen, buf1, this->borderValue);
746  this->lineFunction(buf1, buf2, this->lineLen, destLines[l]);
747  }
748  }
749 
750  return RES_OK;
751  }
752 
753  // Z-Horizontal segment
754  template <class T_in, class lineFunction_T>
756  {
757  size_t w, h, d;
758  imIn.getSize(&w, &h, &d);
759 
760  int nthreads = Core::getInstance()->getNumberOfThreads();
761  lineType *_bufs = this->createAlignedBuffers(2*nthreads, this->lineLen);
762  lineType buf1 = _bufs[0];
763  lineType buf2 = _bufs[nthreads];
764 
765  volType srcSlices = imIn.getSlices();
766  volType destSlices = imOut.getSlices();
767 
768  #ifdef USE_OPEN_MP
769  int tid;
770  #endif // USE_OPEN_MP
771  size_t y;
772 
773  #ifdef USE_OPEN_MP
774  #pragma omp parallel private(tid,buf1,buf2) num_threads(nthreads)
775  #endif // USE_OPEN_MP
776  {
777  #ifdef USE_OPEN_MP
778  tid = omp_get_thread_num();
779  buf1 = _bufs[tid];
780  buf2 = _bufs[tid+nthreads];
781  #pragma omp for
782  #endif
783  for (y=0;y<h;y++)
784  {
785  this->lineFunction(this->borderBuf, srcSlices[0][y], this->lineLen, buf1);
786 
787  for (size_t z=1;z<d;z++)
788  {
789  // Todo: if oddLines...
790  this->lineFunction(srcSlices[z][y], srcSlices[z-1][y], this->lineLen, buf2);
791  this->lineFunction(buf1, buf2, this->lineLen, destSlices[z-1][y]);
792 
793  swap(buf1, buf2);
794  }
795 
796  this->lineFunction(this->borderBuf, buf1, this->lineLen, destSlices[d-1][y]);
797  }
798  }
799 
800  return RES_OK;
801  }
802 
803  template <class T_in, class lineFunction_T>
805  {
806  UINT imHeight = imIn.getHeight();
807  size_t imWidth = imIn.getWidth();
808  volType srcSlices = imIn.getSlices();
809  volType destSlices = imOut.getSlices();
810  sliceType srcLines;
811  sliceType destLines;
812 
813  int tid = 0, nthreads = MIN(Core::getInstance()->getNumberOfThreads(), imHeight/4);
814  nthreads = MAX(nthreads, 1);
815  int nbufs = 4;
816  lineType *_bufs = this->createAlignedBuffers(nbufs*nthreads, this->lineLen);
817  lineType buf1, buf2, firstLineBuf;
818 
819  size_t firstLine, blockSize;
820 
821  for (size_t s=0;s<imIn.getDepth();s++)
822  {
823  srcLines = srcSlices[s];
824  destLines = destSlices[s];
825 
826  #ifdef USE_OPEN_MP
827  #pragma omp parallel private(tid,blockSize,firstLine,buf1,buf2,firstLineBuf) num_threads(nthreads)
828  #endif
829  {
830  #ifdef USE_OPEN_MP
831  tid = omp_get_thread_num();
832  #endif
833  buf1 = _bufs[tid*nbufs];
834  buf2 = _bufs[tid*nbufs+1];
835  firstLineBuf = _bufs[tid*nbufs+2];
836 
837  blockSize = imHeight/nthreads;
838  firstLine = tid*blockSize;
839  if (tid==nthreads-1)
840  blockSize = imHeight-blockSize*tid;
841 
842 
843  // Process first line
844  copyLine<T_in>(srcLines[firstLine], imWidth, buf1);
845  if (firstLine==0)
846  lineFunction(buf1, this->borderBuf, imWidth, buf2);
847  else
848  lineFunction(buf1, srcLines[firstLine-1], imWidth, buf2);
849  lineFunction(srcLines[firstLine], srcLines[firstLine+1], imWidth, buf1);
850  lineFunction(buf1, buf2, imWidth, firstLineBuf);
851 
852  #ifdef USE_OPEN_MP
853  #pragma omp barrier
854  #endif // USE_OPEN_MP
855 
856  for (size_t i = firstLine+1 ; i<firstLine+blockSize-1 ; i++)
857  {
858  lineFunction(srcLines[i], srcLines[i+1], imWidth, buf2);
859  lineFunction(buf1, buf2, imWidth, destLines[i]);
860 
861  swap(buf1, buf2);
862  }
863 
864  if (firstLine+blockSize==imHeight)
865  lineFunction(srcLines[firstLine+blockSize-1], this->borderBuf, imWidth, buf2);
866  else
867  lineFunction(srcLines[firstLine+blockSize-1], srcLines[firstLine+blockSize], imWidth, buf2);
868  lineFunction(buf1, buf2, imWidth, destLines[firstLine+blockSize-1]);
869 
870  #ifdef USE_OPEN_MP
871  #pragma omp barrier
872  #endif // USE_OPEN_MP
873 
874  // finaly write the first line
875  copyLine<T_in>(firstLineBuf, imWidth, destLines[firstLine]);
876 
877  } // #pragma omp parallel
878  }
879  return RES_OK;
880  }
881 
882  template <class T_in, class lineFunction_T>
884  {
885  UINT imHeight = imIn.getHeight();
886  size_t imWidth = imIn.getWidth();
887  volType srcSlices = imIn.getSlices();
888  volType destSlices = imOut.getSlices();
889  sliceType srcLines;
890  sliceType destLines;
891 
892  int tid = 0, nthreads = MIN(Core::getInstance()->getNumberOfThreads(), imHeight/4);
893  nthreads = MAX(nthreads, 1);
894  int nbufs = 6;
895  lineType *_bufs = this->createAlignedBuffers(nbufs*nthreads, this->lineLen);
896  lineType buf1, buf2, buf3, buf4, tmpBuf, firstLineBuf;
897  lineType swap_buf;
898 
899  size_t firstLine, blockSize;
900 
901  for (size_t s=0;s<imIn.getDepth();s++)
902  {
903  srcLines = srcSlices[s];
904  destLines = destSlices[s];
905 
906  #ifdef USE_OPEN_MP
907  #pragma omp parallel private(tid,blockSize,firstLine,buf1,buf2,buf3,buf4,tmpBuf,firstLineBuf,swap_buf) num_threads(nthreads)
908  #endif
909  {
910  #ifdef USE_OPEN_MP
911  tid = omp_get_thread_num();
912  #endif
913  buf1 = _bufs[tid*nbufs];
914  buf2 = _bufs[tid*nbufs+1];
915  buf3 = _bufs[tid*nbufs+2];
916  buf4 = _bufs[tid*nbufs+3];
917  tmpBuf = _bufs[tid*nbufs+4];
918  firstLineBuf = _bufs[tid*nbufs+5];
919 
920  blockSize = imHeight/nthreads;
921  firstLine = tid*blockSize;
922  if (tid==nthreads-1)
923  blockSize = imHeight-blockSize*tid;
924 
925 
926  // Process first line
927  copyLine<T_in>(srcLines[firstLine], imWidth, buf1);
928  this->_exec_shifted_line_2ways(buf1, 1, imWidth, buf4, tmpBuf);
929 
930  copyLine<T_in>(srcLines[firstLine+1], imWidth, buf2);
931 
932  lineFunction(buf4, buf2, imWidth, tmpBuf);
933  if (firstLine==0)
934  lineFunction(this->borderBuf, tmpBuf, imWidth, firstLineBuf);
935  else
936  lineFunction(srcLines[firstLine-1], tmpBuf, imWidth, firstLineBuf);
937 
938  #ifdef USE_OPEN_MP
939  #pragma omp barrier
940  #endif // USE_OPEN_MP
941 
942  for (size_t i = firstLine+2 ; i<firstLine+blockSize ; i++)
943  {
944  copyLine<T_in>(srcLines[i], imWidth, buf3);
945  this->_exec_shifted_line_2ways(buf2, 1, imWidth, buf4, tmpBuf);
946 
947  lineFunction(buf1, buf3, imWidth, tmpBuf);
948  lineFunction(buf4, tmpBuf, imWidth, destLines[i-1]);
949 
950  swap_buf = buf1;
951  buf1 = buf2;
952  buf2 = buf3;
953  buf3 = swap_buf;
954  }
955 
956  this->_exec_shifted_line_2ways(buf2, 1, imWidth, buf4, tmpBuf);
957  lineFunction(buf1, buf4, imWidth, buf4);
958  if (firstLine+blockSize==imHeight)
959  lineFunction(buf4, this->borderBuf, imWidth, destLines[firstLine+blockSize-1]);
960  else
961  lineFunction(buf4, srcLines[firstLine+blockSize], imWidth, destLines[firstLine+blockSize-1]);
962 
963  #ifdef USE_OPEN_MP
964  #pragma omp barrier
965  #endif // USE_OPEN_MP
966 
967  // finaly write the first line
968  copyLine<T_in>(firstLineBuf, imWidth, destLines[firstLine]);
969 
970  } // #pragma omp parallel
971  }
972  return RES_OK;
973  }
974 
975  template <class T_in, class lineFunction_T>
977  {
978  size_t w,h,d ;
979  imIn.getSize (&w, &h, &d) ;
980 
981  volType srcSlices = imIn.getSlices();
982  volType destSlices = imOut.getSlices();
983  sliceType srcLines;
984  sliceType destLines;
985 
986  int tid=0, nthreads = MIN(Core::getInstance()->getNumberOfThreads (), h/4);
987  nthreads = MAX (nthreads, 1) ;
988  int nbufs = 6;
989  lineType *_bufs = this->createAlignedBuffers ( nbufs * nthreads, this->lineLen ) ;
990  lineType buf1, buf2, buf3, buf4, tmp1, firstLineBuf;
991  lineType swap_buf;
992 
993  size_t firstLine, blockSize;
994 
995  srcLines = srcSlices[0];
996  destLines = destSlices[0];
997  #ifdef USE_OPEN_MP
998  #pragma omp parallel private(tid,blockSize,firstLine,buf1,buf2,buf3,buf4,tmp1,firstLineBuf,swap_buf) num_threads(nthreads)
999  #endif // USE_OPEN_MP
1000  {
1001  #ifdef USE_OPEN_MP
1002  tid = omp_get_thread_num();
1003  #endif
1004  buf1 = _bufs[tid*nbufs];
1005  buf2 = _bufs[tid*nbufs+1];
1006  buf3 = _bufs[tid*nbufs+2];
1007  buf4 = _bufs[tid*nbufs+3];
1008  tmp1 = _bufs[tid*nbufs+4];
1009  firstLineBuf = _bufs[tid*nbufs+5];
1010 
1011  blockSize = h/nthreads;
1012  firstLine = tid*blockSize;
1013  if (tid==nthreads-1)
1014  blockSize = h-blockSize*tid;
1015 
1016  // Process first line.
1017  copyLine<T_in>(srcLines[firstLine], w, buf1);
1018  this->_exec_shifted_line_2ways(buf1, 1, w, tmp1, buf4);
1019  lineFunction(this->borderBuf, srcSlices[0][firstLine], w, buf4);
1020  lineFunction(buf4, srcSlices[1][firstLine], w, buf4);
1021  lineFunction(buf4, tmp1, w, buf4);
1022  copyLine<T_in>(srcLines[firstLine+1], w, buf2);
1023  lineFunction(buf4, buf2, w, tmp1);
1024  if (firstLine==0)
1025  lineFunction(this->borderBuf, tmp1, w, firstLineBuf);
1026  else
1027  lineFunction(srcLines[firstLine-1], tmp1, w, firstLineBuf);
1028  #ifdef USE_OPEN_MP
1029  #pragma omp barrier
1030  #endif // USE_OPEN_MP
1031  for (size_t i=firstLine+2; i<firstLine+blockSize; ++i) {
1032  copyLine<T_in>(srcLines[i], w, buf3);
1033  this->_exec_shifted_line_2ways (buf2, 1, w, buf4, tmp1) ;
1034  lineFunction (buf1, buf3, w, tmp1);
1035  lineFunction (buf4, tmp1, w, tmp1);
1036  lineFunction(this->borderBuf, srcSlices[0][i-1], w, buf4);
1037  lineFunction(buf4, srcSlices[1][i-1], w, buf4);
1038  lineFunction (buf4, tmp1, w, destLines[i-1]);
1039 
1040  swap_buf = buf1;
1041  buf1 = buf2;
1042  buf2 = buf3;
1043  buf3 = swap_buf;
1044  }
1045 
1046  this->_exec_shifted_line_2ways (buf2, 1, w, buf4, tmp1);
1047  lineFunction (buf1, buf4, w, tmp1);
1048  lineFunction(this->borderBuf, srcSlices[0][firstLine+blockSize-1], w, buf4);
1049  lineFunction(buf4, srcSlices[1][firstLine+blockSize-1], w, buf4);
1050  lineFunction(tmp1, buf4, w, buf4);
1051  if (firstLine+blockSize == h)
1052  lineFunction (buf4, this->borderBuf, w, destLines[firstLine+blockSize-1]);
1053  else
1054  lineFunction (buf4, srcLines[firstLine+blockSize], w, destLines[firstLine+blockSize-1]);
1055 
1056  #ifdef USE_OPEN_MP
1057  #pragma omp barrier
1058  #endif // USE_OPEN_MP
1059  // finally write the first line
1060  copyLine<T_in>(firstLineBuf, w, destLines[firstLine]);
1061 
1062  for (size_t s=1; s<d-1; ++s) {
1063  srcLines = srcSlices[s];
1064  destLines = destSlices[s];
1065 
1066  blockSize = h/nthreads;
1067  firstLine = tid*blockSize;
1068  if (tid==nthreads-1)
1069  blockSize = h-blockSize*tid;
1070 
1071  // Process first line.
1072  copyLine<T_in>(srcLines[firstLine], w, buf1);
1073  this->_exec_shifted_line_2ways(buf1, 1, w, tmp1, buf4);
1074  lineFunction(srcSlices[s-1][firstLine], srcSlices[s][firstLine], w, buf4);
1075  lineFunction(buf4, srcSlices[s+1][firstLine], w, buf4);
1076  lineFunction(buf4, tmp1, w, buf4);
1077  copyLine<T_in>(srcLines[firstLine+1], w, buf2);
1078  lineFunction(buf4, buf2, w, tmp1);
1079  if (firstLine==0)
1080  lineFunction(this->borderBuf, tmp1, w, firstLineBuf);
1081  else
1082  lineFunction(srcLines[firstLine-1], tmp1, w, firstLineBuf);
1083  #ifdef USE_OPEN_MP
1084  #pragma omp barrier
1085  #endif // USE_OPEN_MP
1086  for (size_t i=firstLine+2; i<firstLine+blockSize; ++i) {
1087  copyLine<T_in>(srcLines[i], w, buf3);
1088  this->_exec_shifted_line_2ways (buf2, 1, w, buf4, tmp1) ;
1089  lineFunction (buf1, buf3, w, tmp1);
1090  lineFunction (buf4, tmp1, w, tmp1);
1091  lineFunction(srcSlices[s-1][i-1], srcSlices[s][i-1], w, buf4);
1092  lineFunction(buf4, srcSlices[s+1][i-1], w, buf4);
1093  lineFunction (buf4, tmp1, w, destLines[i-1]);
1094 
1095  swap_buf = buf1;
1096  buf1 = buf2;
1097  buf2 = buf3;
1098  buf3 = swap_buf;
1099  }
1100 
1101  this->_exec_shifted_line_2ways (buf2, 1, w, buf4, tmp1);
1102  lineFunction (buf1, buf4, w, tmp1);
1103  lineFunction(srcSlices[s-1][firstLine+blockSize-1], srcSlices[s][firstLine+blockSize-1], w, buf4);
1104  lineFunction(buf4, srcSlices[s+1][firstLine+blockSize-1], w, buf4);
1105  lineFunction(tmp1, buf4, w, buf4);
1106  if (firstLine+blockSize == h)
1107  lineFunction (buf4, this->borderBuf, w, destLines[firstLine+blockSize-1]);
1108  else
1109  lineFunction (buf4, srcLines[firstLine+blockSize], w, destLines[firstLine+blockSize-1]);
1110 
1111  #ifdef USE_OPEN_MP
1112  #pragma omp barrier
1113  #endif // USE_OPEN_MP
1114  // finally write the first line
1115  copyLine<T_in>(firstLineBuf, w, destLines[firstLine]);
1116 
1117  }
1118 
1119  srcLines = srcSlices[d-1];
1120  destLines = destSlices[d-1];
1121 
1122  blockSize = h/nthreads;
1123  firstLine = tid*blockSize;
1124  if (tid==nthreads-1)
1125  blockSize = h-blockSize*tid;
1126 
1127  // Process first line.
1128  copyLine<T_in>(srcLines[firstLine], w, buf1);
1129  this->_exec_shifted_line_2ways(buf1, 1, w, tmp1, buf4);
1130  lineFunction(srcSlices[d-2][firstLine], srcSlices[d-1][firstLine], w, buf4);
1131  lineFunction(buf4, this->borderBuf, w, buf4);
1132  lineFunction(buf4, tmp1, w, buf4);
1133  copyLine<T_in>(srcLines[firstLine+1], w, buf2);
1134  lineFunction(buf4, buf2, w, tmp1);
1135  if (firstLine==0)
1136  lineFunction(this->borderBuf, tmp1, w, firstLineBuf);
1137  else
1138  lineFunction(srcLines[firstLine-1], tmp1, w, firstLineBuf);
1139  #ifdef USE_OPEN_MP
1140  #pragma omp barrier
1141  #endif // USE_OPEN_MP
1142  for (size_t i=firstLine+2; i<firstLine+blockSize; ++i) {
1143  copyLine<T_in>(srcLines[i], w, buf3);
1144  this->_exec_shifted_line_2ways (buf2, 1, w, buf4, tmp1) ;
1145  lineFunction (buf1, buf3, w, tmp1);
1146  lineFunction (buf4, tmp1, w, tmp1);
1147  lineFunction(srcSlices[d-2][i-1], srcSlices[d-1][i-1], w, buf4);
1148  lineFunction(buf4, this->borderBuf, w, buf4);
1149  lineFunction (buf4, tmp1, w, destLines[i-1]);
1150 
1151  swap_buf = buf1;
1152  buf1 = buf2;
1153  buf2 = buf3;
1154  buf3 = swap_buf;
1155  }
1156 
1157  this->_exec_shifted_line_2ways (buf2, 1, w, buf4, tmp1);
1158  lineFunction (buf1, buf4, w, tmp1);
1159  lineFunction(srcSlices[d-2][firstLine+blockSize-1], srcSlices[d-1][firstLine+blockSize-1], w, buf4);
1160  lineFunction(buf4, this->borderBuf, w, buf4);
1161  lineFunction(tmp1, buf4, w, buf4);
1162  if (firstLine+blockSize == h)
1163  lineFunction (buf4, this->borderBuf, w, destLines[firstLine+blockSize-1]);
1164  else
1165  lineFunction (buf4, srcLines[firstLine+blockSize], w, destLines[firstLine+blockSize-1]);
1166 
1167  #ifdef USE_OPEN_MP
1168  #pragma omp barrier
1169  #endif // USE_OPEN_MP
1170  // finally write the first line
1171  copyLine<T_in>(firstLineBuf, w, destLines[firstLine]);
1172 
1173  } // #pragma omp parallel
1174 
1175  return RES_OK;
1176  }
1177 
1178  template <class T_in, class lineFunction_T>
1180  {
1181  double nbSquareDbl = (((double) size)/(1+sqrt(2.)));
1182  double nbSquareFloor = floor(nbSquareDbl);
1183  int nbSquare = (int) (((nbSquareDbl - nbSquareFloor) < 0.5f) ? (nbSquareFloor) : (nbSquareFloor+1));
1184 
1185  ASSERT(_exec_single (imIn, imOut, Cross3DSE ())==RES_OK);
1186 
1187  for (size_t i=1; i<size-nbSquare; ++i)
1188  ASSERT(_exec_single(imOut, imOut, Cross3DSE ())==RES_OK);
1189 
1190  for (int i=0; i<nbSquare; ++i)
1191  ASSERT(_exec_single(imOut, imOut, CubeSE ())==RES_OK);
1192 
1193  return RES_OK;
1194  }
1195 
1196 
1197 } // namespace smil
1198 
1199 # endif // _MORPH_IMAGE_OPERATIONS_HXX
static bool isInplaceSafe(const StrElt &se)
Definition: DMorphImageOperations.hxx:268
Definition: DColorConvert.h:38
sliceType getLines() const
Get an array containing the start offset of each line.
Definition: DImage.hpp:114
RES_T copy(const Image< T1 > &imIn, size_t startX, size_t startY, size_t startZ, size_t sizeX, size_t sizeY, size_t sizeZ, Image< T2 > &imOut, size_t outStartX=0, size_t outStartY=0, size_t outStartZ=0)
Copy image.
Definition: DImageArith.hpp:114
Definition: DBaseImage.h:235
Definition: DMorphImageOperations.hpp:133
vector< IntPoint > points
List of neighbor points.
Definition: DStructuringElement.h:125
StrElt homothety(const UINT s) const
Construct and return an homothetic SE with size s.
Definition: DStructuringElement.cpp:78
size_t getHeight() const
Get image height.
Definition: DBaseImage.h:91
Definition: DBaseImageOperations.hxx:39
Base structuring element.
Definition: DStructuringElement.h:51
virtual void modified()
Trigger modified event (allows to force display update)
Definition: DImage.hxx:223
3D Cubic structuring element (26 neighbors).
Definition: DStructuringElement.h:302
size_t getWidth() const
Get image width.
Definition: DBaseImage.h:87
Definition: DTypes.hpp:78
size_t getLineCount() const
Get the number of lines.
Definition: DBaseImage.h:154
size_t getSliceCount() const
Get the number of slices(for 3D images)
Definition: DBaseImage.h:158
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
3D Cross structuring element (6 neighbors).
Definition: DStructuringElement.h:334