Orfeo Toolbox  3.16
itkFlatStructuringElement.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Insight Segmentation & Registration Toolkit
4  Module: $RCSfile: itkFlatStructuringElement.txx,v $
5  Language: C++
6  Date: $Date: 2010-01-22 09:19:36 $
7  Version: $Revision: 1.20 $
8 
9  Copyright (c) Insight Software Consortium. All rights reserved.
10  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11 
12  This software is distributed WITHOUT ANY WARRANTY; without even
13  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14  PURPOSE. See the above copyright notices for more information.
15 
16 =========================================================================*/
17 
18 #ifndef __itkFlatStructuringElement_txx
19 #define __itkFlatStructuringElement_txx
20 #include "vnl/vnl_math.h"
21 #include "itkFlatStructuringElement.h"
22 #include <math.h>
23 #include <vector>
24 
25 #ifndef M_PI
26 #define M_PI vnl_math::pi
27 #endif
28 
29 #include "itkImage.h"
30 #include "itkImageRegionIterator.h"
31 #include "itkFloodFilledSpatialFunctionConditionalIterator.h"
32 #include "itkEllipsoidInteriorExteriorSpatialFunction.h"
33 
34 #include "itkVanHerkGilWermanDilateImageFilter.h"
35 
36 
37 namespace itk
38 {
39 template<unsigned int VDimension>
40 FlatStructuringElement<VDimension> FlatStructuringElement<VDimension>
41 ::Poly(RadiusType radius, unsigned lines)
42 {
43  Self res = Self();
44  res = res.PolySub(Dispatch<VDimension>(), radius, lines);
45  res.SetRadius( radius );
46 #if 0
47  float theta, phi, step;
48  theta = phi = 0;
49  step = M_PI/(lines - 1);
50 
51  while (theta < M_PI)
52  {
53 
54  std::cout << "theta= " << theta << " phi = " << phi << std::endl;
55  LType O = res.mkOffset(phi, theta);
56  std::cout << O << std::endl;
57  if (res.checkParallel(O, res.m_Lines))
58  {
59  std::cout << "Already have this line" << std::endl;
60  }
61  else
62  {
63  res.m_Lines.push_back(O);
64  }
65  theta += step;
66  phi += step;
67  }
68  std::cout << "---------------" << std::endl;
69 #endif
70  res.ComputeBufferFromLines();
71  return(res);
72 
73 }
74 
75 template<unsigned int VDimension>
76 FlatStructuringElement<VDimension> FlatStructuringElement<VDimension>
77 ::PolySub(const Dispatch<2> &, RadiusType radius, unsigned lines) const
78 {
79  // radial decomposition method from "Radial Decomposition of Discs
80  // and Spheres" - CVGIP: Graphical Models and Image Processing
81  //std::cout << "2 dimensions" << std::endl;
82  Self res = Self();
83  res.m_Decomposable = true;
84 
85  unsigned int rr = 0;
86  for (unsigned i=0;i<VDimension;i++)
87  {
88  if (radius[i] > rr) rr = radius[i];
89  }
90  if (lines == 0)
91  {
92  // select some default line values
93  if (rr <= 3) lines=2;
94  else if (rr <= 8) lines=4;
95  else lines=6;
96  }
97  // start with a circle - figure out the length of the structuring
98  // element we need -- This method results in a polygon with 2*lines
99  // sides, each side with length k, where k is the structuring
100  // element length. Therefore the value of k we need to produce the
101  // radius we want is: (M_PI * rr * 2)/(2*lines)
102  float k1 = (M_PI * (float)radius[0])/((float)lines);
103  float k2 = (M_PI * (float)radius[1])/((float)lines);
104  //std::cout << "k= " << k << std::endl;
105  float theta, step;
106  step = M_PI/lines;
107  theta = 0;
108  // just to ensure that we get the last one
109  while (theta <= M_PI/2.0 + 0.0001)
110  {
111  LType O;
112  O[0] = k1 * vcl_cos(theta);
113  O[1] = k2 * vcl_sin(theta);
114  if (!res.checkParallel(O, res.m_Lines))
115  {
116  //std::cout << O << std::endl;
117  res.m_Lines.push_back(O);
118  }
119  O[0] = k1 * vcl_cos(-theta);
120  O[1] = k2 * vcl_sin(-theta);
121  if (!res.checkParallel(O, res.m_Lines))
122  {
123  //std::cout << O << std::endl;
124  res.m_Lines.push_back(O);
125  }
126  theta += step;
127  //std::cout << "theta1 = " << theta << " " << M_PI/2.0 << std::endl;
128  }
129 
130  return(res);
131 }
132 
133 // O[0] = k1 * vcl_cos(phi) * vcl_cos(theta);
134 // O[1] = k2 * vcl_cos(phi) * vcl_sin(theta);
135 // O[2] = k3 * vcl_sin(theta);
136 
137 template<unsigned int VDimension>
138 FlatStructuringElement<VDimension> FlatStructuringElement<VDimension>
139 ::PolySub(const Dispatch<3> &, RadiusType radius, unsigned lines) const
140 {
141  Self res = Self();
142  res.m_Decomposable = true;
143  // std::cout << "3 dimensions" << std::endl;
144  unsigned int rr = 0;
145  int iterations = 1;
146  int faces = lines * 2;
147  for (unsigned i=0;i<VDimension;i++)
148  {
149  if (radius[i] > rr) rr = radius[i];
150  }
151  switch (faces)
152  {
153  case 12:
154  {
155  // dodecahedron
156  float phi = (1.0 + vcl_sqrt(5.0)) / 2.0;
157  float b = 1.0/phi;
158  float c = 2.0 - phi;
159  unsigned facets = 12;
160  typedef std::vector<FacetType> FacetArrayType;
161  FacetArrayType FacetArray;
162  FacetArray.resize( facets );
163  // set up vectors normal to the faces - only put in 3 points for
164  // each face:
165  // face 1
166  LType PP(0.0);
167  FacetType Fc;
168  b /= 2.0;
169  c /= 2.0;
170 
171  PP[0]=c;PP[1]=0;PP[2]=0.5;
172  Fc.P1 = PP;
173  PP[0]=-c;PP[1]=0;PP[2]=0.5;
174  Fc.P2 = PP;
175  PP[0]=-b;PP[1]=b;PP[2]=b;
176  Fc.P3 = PP;
177  FacetArray[0] = Fc;
178 
179  PP[0]=-c;PP[1]=0;PP[2]=0.5;
180  Fc.P1 = PP;
181  PP[0]=c;PP[1]=0;PP[2]=0.5;
182  Fc.P2 = PP;
183  PP[0]=b;PP[1]=-b;PP[2]=b;
184  Fc.P3 = PP;
185  FacetArray[1] = Fc;
186 
187  PP[0]=c;PP[1]=0;PP[2]=-0.5;
188  Fc.P1 = PP;
189  PP[0]=-c;PP[1]=0;PP[2]=-0.5;
190  Fc.P2 = PP;
191  PP[0]=-b;PP[1]=-b;PP[2]=-b;
192  Fc.P3 = PP;
193  FacetArray[2] = Fc;
194 
195  PP[0]=-c;PP[1]=0;PP[2]=-0.5;
196  Fc.P1 = PP;
197  PP[0]=c;PP[1]=0;PP[2]=-0.5;
198  Fc.P2 = PP;
199  PP[0]=b;PP[1]=b;PP[2]=-b;
200  Fc.P3 = PP;
201  FacetArray[3] = Fc;
202 
203  PP[0]=0;PP[1]=0.5;PP[2]=-c;
204  Fc.P1 = PP;
205  PP[0]=0;PP[1]=0.5;PP[2]=c;
206  Fc.P2 = PP;
207  PP[0]=b;PP[1]=b;PP[2]=b;
208  Fc.P3 = PP;
209  FacetArray[4] = Fc;
210 
211  PP[0]=0;PP[1]=0.5;PP[2]=c;
212  Fc.P1 = PP;
213  PP[0]=0;PP[1]=0.5;PP[2]=-c;
214  Fc.P2 = PP;
215  PP[0]=-b;PP[1]=b;PP[2]=-b;
216  Fc.P3 = PP;
217  FacetArray[5] = Fc;
218 
219  PP[0]=0;PP[1]=-0.5;PP[2]=-c;
220  Fc.P1 = PP;
221  PP[0]=0;PP[1]=-0.5;PP[2]=c;
222  Fc.P2 = PP;
223  PP[0]=-b;PP[1]=-b;PP[2]=b;
224  Fc.P3 = PP;
225  FacetArray[6] = Fc;
226 
227  PP[0]=0;PP[1]=-0.5;PP[2]=c;
228  Fc.P1 = PP;
229  PP[0]=0;PP[1]=-0.5;PP[2]=-c;
230  Fc.P2 = PP;
231  PP[0]=b;PP[1]=-b;PP[2]=-b;
232  Fc.P3 = PP;
233  FacetArray[7] = Fc;
234 
235  PP[0]=0.5;PP[1]=c;PP[2]=0;
236  Fc.P1 = PP;
237  PP[0]=0.5;PP[1]=-c;PP[2]=0;
238  Fc.P2 = PP;
239  PP[0]=b;PP[1]=-b;PP[2]=b;
240  Fc.P3 = PP;
241  FacetArray[8] = Fc;
242 
243  PP[0]=0.5;PP[1]=-c;PP[2]=0;
244  Fc.P1 = PP;
245  PP[0]=0.5;PP[1]=c;PP[2]=0;
246  Fc.P2 = PP;
247  PP[0]=b;PP[1]=b;PP[2]=-b;
248  Fc.P3 = PP;
249  FacetArray[9] = Fc;
250 
251  PP[0]=-0.5;PP[1]=c;PP[2]=0;
252  Fc.P1 = PP;
253  PP[0]=-0.5;PP[1]=-c;PP[2]=0;
254  Fc.P2 = PP;
255  PP[0]=-b;PP[1]=-b;PP[2]=-b;
256  Fc.P3 = PP;
257  FacetArray[10] = Fc;
258 
259  PP[0]=-0.5;PP[1]=-c;PP[2]=0;
260  Fc.P1 = PP;
261  PP[0]=-0.5;PP[1]=c;PP[2]=0;
262  Fc.P2 = PP;
263  PP[0]=-b;PP[1]=b;PP[2]=b;
264  Fc.P3 = PP;
265  FacetArray[11] = Fc;
266  for (unsigned j = 0; j < facets; j++)
267  {
268  // Find a line perpendicular to each face
269  LType L, A, B;
270  A = FacetArray[j].P2 - FacetArray[j].P1;
271  B = FacetArray[j].P3 - FacetArray[j].P1;
272  L[0] = A[1]*B[2] - B[1]*A[2];
273  L[1] = B[0]*A[2] - A[0]*B[2];
274  L[2] = A[0]*B[1] - B[0]*A[1];
275 
276  L.Normalize();
277  // Scale to required length
278  L *= rr;
279  if (!res.checkParallel(L, res.m_Lines))
280  {
281  res.m_Lines.push_back(L);
282  }
283  }
284  return(res);
285 
286  }
287  break;
288  case 14:
289  {
290  // cube with the corners cut off
291  LType A;
292  // The axes
293  A[0]=1;A[1]=0;A[2]=0;
294  A *= rr;
295  res.m_Lines.push_back(A);
296  A[0]=0;A[1]=1;A[2]=0;
297  A *= rr;
298  res.m_Lines.push_back(A);
299  A[0]=0;A[1]=0;A[2]=1;
300  A *= rr;
301  res.m_Lines.push_back(A);
302  // Diagonals
303  A[0]=1;A[1]=1;A[2]=1;
304  A.Normalize();
305  A *= rr;
306  res.m_Lines.push_back(A);
307 
308  A[0]=-1;A[1]=1;A[2]=1;
309  A.Normalize();
310  A *= rr;
311  res.m_Lines.push_back(A);
312 
313  A[0]=1;A[1]=-1;A[2]=1;
314  A.Normalize();
315  A *= rr;
316  res.m_Lines.push_back(A);
317 
318  A[0]=-1;A[1]=-1;A[2]=1;
319  A.Normalize();
320  A *= rr;
321  res.m_Lines.push_back(A);
322  return(res);
323  }
324  break;
325  case 20:
326  {
327  // Icosahedron
328  float phi = (1.0 + vcl_sqrt(5.0)) / 2.0;
329  float a = 0.5;
330  float b = 1.0/(2.0*phi);
331  unsigned facets = 20;
332  typedef std::vector<FacetType> FacetArrayType;
333  FacetArrayType FacetArray;
334  FacetArray.resize( facets );
335  // set up vectors normal to the faces - only put in 3 points for
336  // each face:
337  // face 1
338  LType PP(0.0);
339  FacetType Fc;
340 
341  PP[0]=0;PP[1]=b;PP[2]=-a;
342  Fc.P1 = PP;
343  PP[0]=b;PP[1]=a;PP[2]=0;
344  Fc.P2 = PP;
345  PP[0]=-b;PP[1]=a;PP[2]=0;
346  Fc.P3 = PP;
347  FacetArray[0] = Fc;
348 
349  PP[0]=0;PP[1]=b;PP[2]=a;
350  Fc.P1 = PP;
351  PP[0]=-b;PP[1]=a;PP[2]=0;
352  Fc.P2 = PP;
353  PP[0]=b;PP[1]=a;PP[2]=0;
354  Fc.P3 = PP;
355  FacetArray[1] = Fc;
356 
357  PP[0]=0;PP[1]=b;PP[2]=a;
358  Fc.P1 = PP;
359  PP[0]=0;PP[1]=-b;PP[2]=a;
360  Fc.P2 = PP;
361  PP[0]=-a;PP[1]=0;PP[2]=b;
362  Fc.P3 = PP;
363  FacetArray[2] = Fc;
364 
365  PP[0]=0;PP[1]=b;PP[2]=a;
366  Fc.P1 = PP;
367  PP[0]=a;PP[1]=0;PP[2]=b;
368  Fc.P2 = PP;
369  PP[0]=0;PP[1]=-b;PP[2]=a;
370  Fc.P3 = PP;
371  FacetArray[3] = Fc;
372 
373  PP[0]=0;PP[1]=b;PP[2]=-a;
374  Fc.P1 = PP;
375  PP[0]=0;PP[1]=-b;PP[2]=-a;
376  Fc.P2 = PP;
377  PP[0]=a;PP[1]=0;PP[2]=-b;
378  Fc.P3 = PP;
379  FacetArray[4] = Fc;
380 
381  PP[0]=0;PP[1]=b;PP[2]=-a;
382  Fc.P1 = PP;
383  PP[0]=-a;PP[1]=0;PP[2]=-b;
384  Fc.P2 = PP;
385  PP[0]=0;PP[1]=-b;PP[2]=-a;
386  Fc.P3 = PP;
387  FacetArray[5] = Fc;
388 
389  PP[0]=0;PP[1]=-b;PP[2]=a;
390  Fc.P1 = PP;
391  PP[0]=b;PP[1]=-a;PP[2]=0;
392  Fc.P2 = PP;
393  PP[0]=-b;PP[1]=-a;PP[2]=0;
394  Fc.P3 = PP;
395  FacetArray[6] = Fc;
396 
397  PP[0]=0;PP[1]=-b;PP[2]=-a;
398  Fc.P1 = PP;
399  PP[0]=-b;PP[1]=-a;PP[2]=0;
400  Fc.P2 = PP;
401  PP[0]=b;PP[1]=-a;PP[2]=0;
402  Fc.P3 = PP;
403  FacetArray[7] = Fc;
404 
405  PP[0]=-b;PP[1]=a;PP[2]=0;
406  Fc.P1 = PP;
407  PP[0]=-a;PP[1]=0;PP[2]=b;
408  Fc.P2 = PP;
409  PP[0]=-a;PP[1]=0;PP[2]=-b;
410  Fc.P3 = PP;
411  FacetArray[8] = Fc;
412 
413  PP[0]=-b;PP[1]=-a;PP[2]=0;
414  Fc.P1 = PP;
415  PP[0]=-a;PP[1]=0;PP[2]=-b;
416  Fc.P2 = PP;
417  PP[0]=-a;PP[1]=0;PP[2]=b;
418  Fc.P3 = PP;
419  FacetArray[9] = Fc;
420 
421  PP[0]=b;PP[1]=a;PP[2]=0;
422  Fc.P1 = PP;
423  PP[0]=a;PP[1]=0;PP[2]=-b;
424  Fc.P2 = PP;
425  PP[0]=a;PP[1]=0;PP[2]=b;
426  Fc.P3 = PP;
427  FacetArray[10] = Fc;
428 
429  PP[0]=b;PP[1]=-a;PP[2]=0;
430  Fc.P1 = PP;
431  PP[0]=a;PP[1]=0;PP[2]=b;
432  Fc.P2 = PP;
433  PP[0]=a;PP[1]=0;PP[2]=-b;
434  Fc.P3 = PP;
435  FacetArray[11] = Fc;
436 
437  PP[0]=0;PP[1]=b;PP[2]=a;
438  Fc.P1 = PP;
439  PP[0]=-a;PP[1]=0;PP[2]=b;
440  Fc.P2 = PP;
441  PP[0]=-b;PP[1]=a;PP[2]=0;
442  Fc.P3 = PP;
443  FacetArray[12] = Fc;
444 
445  PP[0]=0;PP[1]=b;PP[2]=a;
446  Fc.P1 = PP;
447  PP[0]=b;PP[1]=a;PP[2]=0;
448  Fc.P2 = PP;
449  PP[0]=a;PP[1]=0;PP[2]=b;
450  Fc.P3 = PP;
451  FacetArray[13] = Fc;
452 
453  PP[0]=0;PP[1]=b;PP[2]=-a;
454  Fc.P1 = PP;
455  PP[0]=-b;PP[1]=a;PP[2]=0;
456  Fc.P2 = PP;
457  PP[0]=-a;PP[1]=0;PP[2]=-b;
458  Fc.P3 = PP;
459  FacetArray[14] = Fc;
460 
461  PP[0]=0;PP[1]=b;PP[2]=-a;
462  Fc.P1 = PP;
463  PP[0]=a;PP[1]=0;PP[2]=-b;
464  Fc.P2 = PP;
465  PP[0]=b;PP[1]=a;PP[2]=0;
466  Fc.P3 = PP;
467  FacetArray[15] = Fc;
468 
469  PP[0]=0;PP[1]=-b;PP[2]=-a;
470  Fc.P1 = PP;
471  PP[0]=-a;PP[1]=0;PP[2]=-b;
472  Fc.P2 = PP;
473  PP[0]=-b;PP[1]=-a;PP[2]=0;
474  Fc.P3 = PP;
475  FacetArray[16] = Fc;
476 
477  PP[0]=0;PP[1]=-b;PP[2]=-a;
478  Fc.P1 = PP;
479  PP[0]=b;PP[1]=-a;PP[2]=0;
480  Fc.P2 = PP;
481  PP[0]=a;PP[1]=0;PP[2]=-b;
482  Fc.P3 = PP;
483  FacetArray[17] = Fc;
484 
485  PP[0]=0;PP[1]=-b;PP[2]=a;
486  Fc.P1 = PP;
487  PP[0]=-b;PP[1]=-a;PP[2]=0;
488  Fc.P2 = PP;
489  PP[0]=-a;PP[1]=0;PP[2]=b;
490  Fc.P3 = PP;
491  FacetArray[18] = Fc;
492 
493  PP[0]=0;PP[1]=-b;PP[2]=a;
494  Fc.P1 = PP;
495  PP[0]=a;PP[1]=0;PP[2]=b;
496  Fc.P2 = PP;
497  PP[0]=b;PP[1]=-a;PP[2]=0;
498  Fc.P3 = PP;
499  FacetArray[19] = Fc;
500 
501  for (unsigned j = 0; j < facets; j++)
502  {
503  // Find a line perpendicular to each face
504  LType L, A, B;
505  A = FacetArray[j].P2 - FacetArray[j].P1;
506  B = FacetArray[j].P3 - FacetArray[j].P1;
507  L[0] = A[1]*B[2] - B[1]*A[2];
508  L[1] = B[0]*A[2] - A[0]*B[2];
509  L[2] = A[0]*B[1] - B[0]*A[1];
510 
511  L.Normalize();
512  // Scale to required length
513  L *= rr;
514  if (!res.checkParallel(L, res.m_Lines))
515  {
516  res.m_Lines.push_back(L);
517  }
518  }
519  return(res);
520  }
521  break;
522  case 32:
523  {
524  iterations = 1;
525  // 2 iterations leads to 128 faces, which is too many
526  // subdivision of octahedron
527  // create triangular facet approximation to a sphere - begin with
528  // unit sphere
529  // total number of facets is 8 * (4^iterations)
530  unsigned int facets = 8 * (int)vcl_pow((double)4, iterations);
531  float sqrt2 = vcl_sqrt(2.0);
532  // std::cout << facets << " facets" << std::endl;
533  typedef std::vector<FacetType> FacetArrayType;
534  FacetArrayType FacetArray;
535  FacetArray.resize( facets );
536 
537  // original corners of octahedron
538  LType P0(0.0), P1(0.0), P2(0.0), P3(0.0), P4(0.0), P5(0.0);
539  P0[0]=0; P0[1]=0; P0[2]=1;
540  P1[0]=0; P1[1]=0; P1[2]=-1;
541  P2[0]=-1.0/sqrt2;P2[1]=-1/sqrt2;P2[2]=0;
542  P3[0]=1/sqrt2; P3[1]=-1/sqrt2;P3[2]=0;
543  P4[0]=1/sqrt2; P4[1]=1/sqrt2; P4[2]=0;
544  P5[0]=-1/sqrt2; P5[1]=1/sqrt2; P5[2]=0;
545 
546  FacetType F0, F1, F2, F3, F4, F5, F6, F7;
547  F0.P1 = P0; F0.P2 = P3; F0.P3 = P4;
548  F1.P1 = P0; F1.P2 = P4; F1.P3 = P5;
549  F2.P1 = P0; F2.P2 = P5; F2.P3 = P2;
550  F3.P1 = P0; F3.P2 = P2; F3.P3 = P3;
551  F4.P1 = P1; F4.P2 = P4; F4.P3 = P3;
552  F5.P1 = P1; F5.P2 = P5; F5.P3 = P4;
553  F6.P1 = P1; F6.P2 = P2; F6.P3 = P5;
554  F7.P1 = P1; F7.P2 = P3; F7.P3 = P2;
555 
556  FacetArray[0] = F0;
557  FacetArray[1] = F1;
558  FacetArray[2] = F2;
559  FacetArray[3] = F3;
560  FacetArray[4] = F4;
561  FacetArray[5] = F5;
562  FacetArray[6] = F6;
563  FacetArray[7] = F7;
564  int pos = 8;
565  // now subdivide the octahedron
566  for (unsigned it = 0; it<(unsigned)iterations; it++)
567  {
568  // Bisect edges and move to sphere
569  unsigned ntold = pos;
570  for (unsigned i = 0; i < ntold; i++)
571  {
572  LType Pa, Pb, Pc;
573  for (unsigned d = 0; d< VDimension;d++)
574  {
575  Pa[d] = (FacetArray[i].P1[d] + FacetArray[i].P2[d])/2;
576  Pb[d] = (FacetArray[i].P2[d] + FacetArray[i].P3[d])/2;
577  Pc[d] = (FacetArray[i].P3[d] + FacetArray[i].P1[d])/2;
578  }
579  Pa.Normalize();
580  Pb.Normalize();
581  Pc.Normalize();
582  FacetArray[pos].P1 = FacetArray[i].P1;
583  FacetArray[pos].P2 = Pa;
584  FacetArray[pos].P3 = Pc;
585  pos++;
586  FacetArray[pos].P1 = Pa;
587  FacetArray[pos].P2 = FacetArray[i].P2;
588  FacetArray[pos].P3 = Pb;
589  pos++;
590  FacetArray[pos].P1 = Pb;
591  FacetArray[pos].P2 = FacetArray[i].P3;
592  FacetArray[pos].P3 = Pc;
593  pos++;
594  FacetArray[i].P1 = Pa;
595  FacetArray[i].P2 = Pb;
596  FacetArray[i].P3 = Pc;
597  }
598  }
599 
600 
601  for (unsigned j = 0; j < facets; j++)
602  {
603  // Find a line perpendicular to each face
604  LType L, A, B;
605  A = FacetArray[j].P2 - FacetArray[j].P1;
606  B = FacetArray[j].P3 - FacetArray[j].P1;
607  L[0] = A[1]*B[2] - B[1]*A[2];
608  L[1] = B[0]*A[2] - A[0]*B[2];
609  L[2] = A[0]*B[1] - B[0]*A[1];
610 
611  L.Normalize();
612  // Scale to required length
613  L *= rr;
614  if (!res.checkParallel(L, res.m_Lines))
615  {
616  res.m_Lines.push_back(L);
617  }
618  }
619  return(res);
620  }
621  break;
622  default:
623  std::cout << "Unsupported number of lines" << std::endl;
624  return(res);
625  }
626 }
627 
628 template<unsigned int VDimension>
629 FlatStructuringElement<VDimension> FlatStructuringElement<VDimension>
630 ::PolySub(const DispatchBase &, RadiusType , unsigned) const
631 {
632  Self res = Self();
633  res.m_Decomposable = true;
634  std::cout << "Don't know how to deal with this many dimensions" << std::endl;
635  return(res);
636 }
637 
638 
639 template<unsigned int VDimension>
640 FlatStructuringElement<VDimension> FlatStructuringElement<VDimension>
641 ::Box(RadiusType radius)
642 {
643  // this should work for any number of dimensions
644  Self res = Self();
645  res.m_Decomposable = true;
646  res.SetRadius( radius );
647  for (unsigned i = 0;i<VDimension;i++)
648  {
649  if (radius[i] != 0)
650  {
651  LType L;
652  L.Fill(0);
653  L[i] = radius[i]*2+1;
654  res.m_Lines.push_back(L);
655  }
656  }
657  // this doesn't work if one of the dimensions is zero. Suspect an
658  //"inconsistency" in the way
659 // res.ComputeBufferFromLines();
660  Iterator kernel_it;
661  for( kernel_it=res.Begin(); kernel_it != res.End(); ++kernel_it )
662  {
663  *kernel_it= true;
664  }
665 
666  return(res);
667 }
668 
669 
670 template<unsigned int VDimension>
671 FlatStructuringElement<VDimension> FlatStructuringElement<VDimension>
672 ::Cross(RadiusType radius)
673 {
674  // this should work for any number of dimensions
675  Self res = Self();
676  res.m_Decomposable = false;
677  res.SetRadius( radius );
678  Iterator kernel_it;
679  for( kernel_it=res.Begin(); kernel_it != res.End(); ++kernel_it )
680  {
681  *kernel_it= false;
682  }
683  for( int d = 0; d<(int)VDimension; d++ )
684  {
685  OffsetType o;
686  o.Fill( 0 );
687  for( int i=-(int)radius[d]; i<=(int)radius[d]; i++ )
688  {
689  o[d] = i;
690  res[o] = true;
691  }
692  }
693 
694  return(res);
695 }
696 
697 
698 template<unsigned int VDimension>
699 FlatStructuringElement<VDimension> FlatStructuringElement<VDimension>
700 ::Ball(RadiusType radius)
701 {
702  Self res = Self();
703  res.SetRadius( radius );
704  res.m_Decomposable = false;
705 
706  unsigned int i;
707 
708  // Image typedef
709  typedef Image<bool, VDimension> ImageType;
710 
711  // Create an image to hold the ellipsoid
712  //
713  typename ImageType::Pointer sourceImage = ImageType::New();
714  typename ImageType::RegionType region;
715  RadiusType size = radius;
716  for( i=0; i<(int)VDimension; i++ )
717  {
718  size[i] = 2*size[i] + 1;
719  }
720  region.SetSize( size );
721 
722  sourceImage->SetRegions( region );
723  sourceImage->Allocate();
724  // sourceImage->Print( std::cout );
725 
726  // Set the background to be zero
727  //
728  ImageRegionIterator<ImageType> it(sourceImage, region);
729 
730  for(it.GoToBegin(); !it.IsAtEnd(); ++it)
731  {
732  it.Set(false);
733  }
734 
735 
736  // Create the ellipsoid
737  //
738 
739  // Ellipsoid spatial function typedef
740  typedef EllipsoidInteriorExteriorSpatialFunction<VDimension> EllipsoidType;
741 
742  // Create an ellipsoid spatial function for the source image
743  typename EllipsoidType::Pointer spatialFunction = EllipsoidType::New();
744 
745  // Define and set the axes lengths for the ellipsoid
746  typename EllipsoidType::InputType axes;
747  for (i=0; i < VDimension; i++)
748  {
749  axes[i] = res.GetSize(i);
750  }
751  spatialFunction->SetAxes( axes );
752 
753  // Define and set the center of the ellipsoid in physical space
754  typename EllipsoidType::InputType center;
755  for (i=0; i < VDimension; i++)
756  {
757  // put the center of ellipse in the middle of the center pixel
758  center[i] = res.GetRadius(i) + 0.5;
759  }
760  spatialFunction->SetCenter( center );
761 
762  // Define the orientations of the ellipsoid axes, for now, we'll use
763  // the identify matrix
764  typename EllipsoidType::OrientationType orientations;
765  orientations.fill( 0.0 );
766  orientations.fill_diagonal( 1.0 );
767  spatialFunction->SetOrientations(orientations);
768 
769  typename ImageType::IndexType seed;
770  for (i=0; i < VDimension; i++)
771  {
772  seed[i] = res.GetRadius(i);
773  }
774  FloodFilledSpatialFunctionConditionalIterator<ImageType, EllipsoidType>
775  sfi = FloodFilledSpatialFunctionConditionalIterator<ImageType,
776  EllipsoidType>(sourceImage, spatialFunction, seed);
777  sfi.SetCenterInclusionStrategy();
778 
779  // Iterate through the entire image and set interior pixels to 1
780  for(; !sfi.IsAtEnd(); ++sfi)
781  {
782  sfi.Set(true);
783  }
784 
785 
786  // Copy the ellipsoid into the kernel
787  //
788  Iterator kernel_it;
789  for (it.GoToBegin(), kernel_it=res.Begin();!it.IsAtEnd();++it,++kernel_it)
790  {
791  *kernel_it = it.Get();
792  }
793 
794  // Clean up
795  // ...temporary image should be cleaned up by SmartPointers automatically
796 
797  return res;
798 }
799 
800 
801 template<unsigned int NDimension>
802 FlatStructuringElement<NDimension>
803 FlatStructuringElement<NDimension>
804 ::Annulus(RadiusType radius,
805  unsigned int thickness,
806  bool includeCenter)
807 {
808  Self result = Self();
809  result.SetRadius( radius );
810  result.m_Decomposable = false;
811 
812  // Image typedef
813  typedef Image<bool,NDimension> ImageType;
814 
815  // Create an image to hold the ellipsoid
816  //
817  typename ImageType::Pointer kernelImage = ImageType::New();
818  typename ImageType::RegionType region;
819  RadiusType size = radius;
820  for(unsigned int i=0; i<NDimension; i++ )
821  {
822  size[i] = 2*size[i] + 1;
823  }
824  region.SetSize( size );
825 
826  kernelImage->SetRegions( region );
827  kernelImage->Allocate();
828 
829  // Set the background to be zero
830  //
831  ImageRegionIterator<ImageType> it(kernelImage, region);
832 
833  for(it.GoToBegin(); !it.IsAtEnd(); ++it)
834  {
835  it.Set(false);
836  }
837 
838  // Create two ellipsoids
839  //
840 
841  // Ellipsoid spatial function typedef
842  typedef EllipsoidInteriorExteriorSpatialFunction<NDimension> EllipsoidType;
843 
844  // Create an ellipsoid spatial function for the source image
845  typename EllipsoidType::Pointer ellipsoidOuter = EllipsoidType::New();
846  typename EllipsoidType::Pointer ellipsoidInner = EllipsoidType::New();
847 
848  // Define and set the axes lengths for the ellipsoid
849  typename EllipsoidType::InputType axesOuter;
850  typename EllipsoidType::InputType axesInner;
851  for (unsigned int i=0; i < NDimension; i++)
852  {
853  axesOuter[i] = result.GetSize(i);
854  axesInner[i] = std::max( 2*(long)radius[i] + 1 - 2*(long)thickness, (long)1 );
855  }
856  ellipsoidOuter->SetAxes( axesOuter );
857  ellipsoidInner->SetAxes( axesInner );
858 
859  // Define and set the center of the ellipsoid in physical space
860  typename EllipsoidType::InputType center;
861  for (unsigned int i=0; i < NDimension; i++)
862  {
863  // put the center of ellipse in the middle of the center pixel
864  center[i] = result.GetRadius(i) + 0.5;
865  }
866  ellipsoidOuter->SetCenter( center );
867  ellipsoidInner->SetCenter( center );
868 
869  // Define the orientations of the ellipsoid axes, for now, we'll use
870  // the identity matrix
871  typename EllipsoidType::OrientationType orientations;
872  orientations.fill( 0.0 );
873  orientations.fill_diagonal( 1.0 );
874  ellipsoidOuter->SetOrientations( orientations );
875  ellipsoidInner->SetOrientations( orientations );
876 
877  // Create the starting seed
878  typename ImageType::IndexType seed;
879  for (unsigned int i=0; i < NDimension; i++)
880  {
881  seed[i] = result.GetRadius(i);
882  }
883 
884  // Define the iterators for each ellipsoid
885  typedef FloodFilledSpatialFunctionConditionalIterator<ImageType, EllipsoidType>
886  FloodIteratorType;
887  FloodIteratorType itEllipsoidOuter =
888  FloodIteratorType( kernelImage, ellipsoidOuter, seed );
889  FloodIteratorType itEllipsoidInner =
890  FloodIteratorType( kernelImage, ellipsoidInner, seed );
891  itEllipsoidOuter.SetCenterInclusionStrategy();
892  itEllipsoidInner.SetCenterInclusionStrategy();
893 
894  // Iterate through the image and set all outer pixels to 'ON'
895  for(; !itEllipsoidOuter.IsAtEnd(); ++itEllipsoidOuter)
896  {
897  itEllipsoidOuter.Set(true);
898  }
899  // Iterate through the image and set all inner pixels to 'OFF'
900  for(; !itEllipsoidInner.IsAtEnd(); ++itEllipsoidInner)
901  {
902  itEllipsoidInner.Set(false);
903  }
904 
905  // Set center pixel if included
906  kernelImage->SetPixel( seed, includeCenter );
907 
908  // Copy the annulus into the kernel
909  //
910  Iterator kernel_it;
911  for (it.GoToBegin(), kernel_it=result.Begin(); !it.IsAtEnd(); ++it,++kernel_it)
912  {
913  *kernel_it = it.Get();
914  }
915 
916  // Clean up
917  // ...temporary image should be cleaned up by SmartPointers automatically
918 
919  return result;
920 }
921 template<unsigned int VDimension>
922 bool
923 FlatStructuringElement<VDimension>::
924 checkParallel(LType NewVec, DecompType Lines)
925 {
926  LType NN = NewVec;
927  NN.Normalize();
928  for (unsigned i = 0; i < Lines.size(); i++)
929  {
930  LType LL = Lines[i];
931  LL.Normalize();
932  float L = NN*LL;
933  if ((1.0 - vcl_fabs(L)) < 0.000001) return(true);
934  }
935  return(false);
936 }
937 
938 template<unsigned int VDimension>
939 void FlatStructuringElement< VDimension >
940 ::PrintSelf(std::ostream &os, Indent indent) const
941 {
942  //Superclass::PrintSelf(os, indent);
943  if (m_Decomposable)
944  {
945  os << indent << "SE decomposition:" << std::endl;
946  for (unsigned i = 0;i < m_Lines.size(); i++)
947  {
948  os << indent << m_Lines[i] << std::endl;
949  }
950  }
951 }
952 
953 
954 template<unsigned int VDimension>
955 void
956 FlatStructuringElement<VDimension>::
957 ComputeBufferFromLines()
958 {
959  if( m_Decomposable )
960  {
961 /* itkExceptionMacro("Element must be decomposable.");*/
962  }
963 
964  // create an image with a single pixel in the center which will be dilated
965  // by the structuring lines (with AnchorDilateImageFilter) so the content
966  // of the buffer will reflect the shape of the structuring element
967 
968  // Image typedef
969  typedef Image<bool, VDimension> ImageType;
970 
971  // Create an image to hold the ellipsoid
972  //
973  typename ImageType::Pointer sourceImage = ImageType::New();
974  typename ImageType::RegionType region;
975  RadiusType size = this->GetRadius();
976  for( int i=0; i<(int)VDimension; i++ )
977  {
978  size[i] = 2*size[i] + 1;
979  }
980  region.SetSize( size );
981  sourceImage->SetRegions( region );
982  sourceImage->Allocate();
983 
984  // sourceImage->Print(std::cout);
985 
986  // Set the background to be zero
987  //
988  ImageRegionIterator<ImageType> it( sourceImage, region );
989 
990  for( it.GoToBegin(); !it.IsAtEnd(); ++it )
991  {
992  it.Set( false );
993  }
994 
995  // set the center pixel to 1
996  typename ImageType::IndexType center;
997  for( int i=0; i<(int)VDimension; i++)
998  {
999  center[i] = this->GetRadius()[i];
1000  }
1001  sourceImage->SetPixel( center, true );
1002 
1003  // initialize the kernel with everything to false, to avoid warnings in
1004  // valgrind in the SetKernel() method
1005  for( Iterator kernel_it=this->Begin(); kernel_it != this->End(); ++kernel_it )
1006  {
1007  *kernel_it = false;
1008  }
1009 
1010  // dilate the pixel
1011  typedef VanHerkGilWermanDilateImageFilter<ImageType, Self> DilateType;
1012  typename DilateType::Pointer dilate = DilateType::New();
1013  dilate->SetInput( sourceImage );
1014  dilate->SetKernel( *this );
1015  dilate->Update();
1016 
1017  // copy back the image to the kernel
1018  ImageRegionIterator<ImageType> oit( dilate->GetOutput(), region );
1019  Iterator kernel_it;
1020  for( oit.GoToBegin(), kernel_it=this->Begin(); !oit.IsAtEnd(); ++oit,++kernel_it )
1021  {
1022  *kernel_it = oit.Get();
1023  }
1024 
1025 }
1026 
1027 
1028 }
1029 
1030 #endif
Generated at Sat Feb 2 2013 23:38:18 for Orfeo Toolbox with doxygen 1.8.1.1