Orfeo Toolbox  4.0
otbPolygonListToRCC8GraphFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 #ifndef __otbPolygonListToRCC8GraphFilter_txx
19 #define __otbPolygonListToRCC8GraphFilter_txx
20 
22 #include "itkProgressReporter.h"
23 
24 namespace otb
25 {
29 template <class TPolygonList, class TOutputGraph>
32 {
33  this->SetNumberOfRequiredInputs(1);
34  m_Optimisation = false;
35  m_UseInverted = false;
36 }
40 template <class TPolygonList, class TOutputGraph>
43 {}
44 
45 template <class TPolygonList, class TOutputGraph>
46 void
49 {
50  // Process object is not const-correct so the const_cast is required here
52  const_cast<PolygonListType *>(input));
53 }
54 
55 template <class TPolygonList, class TOutputGraph>
59 {
60  if (this->GetNumberOfInputs() < 1)
61  {
62  return 0;
63  }
64 
65  return static_cast<const TPolygonList *>
67 }
68 
69 template <class TPolygonList, class TOutputGraph>
70 unsigned int
73 {
74  return m_Accumulator[val];
75 }
76 
77 template <class TPolygonList, class TOutputGraph>
78 unsigned int
81 {
82  unsigned int result = 0;
83  for (unsigned int i = 0; i < 8; ++i)
84  {
85  result += m_Accumulator[i];
86 
87  }
88  return result;
89 }
90 template <class TPolygonList, class TOutputGraph>
92 ::KnowledgeStateType
95 {
96  m_Accumulator[0] = 0;
97  m_Accumulator[1] = 0;
98  m_Accumulator[2] = 0;
99  m_Accumulator[3] = 0;
100  m_Accumulator[4] = 0;
101  m_Accumulator[5] = 0;
102  m_Accumulator[6] = 0;
103  m_Accumulator[7] = 0;
104 
105  // otbMsgDebugMacro(<<"RCC8GraphFilter: entering GetKnowledge method.");
106  // This is the RCC8 composition table
107  const int knowledge[8][8]
108  =
109  { {-3, -2, -2, -2, 0, -2, 0, 0}, {-1, -3, -2, -3, -1, -3, 0, 1}, {-1, -1, -3, -3, -1, -3, -1, 2}, { 0, -1, -2, -3, -3, 5, -1, 3}, {-1, -1, -1, -3, -1, -3, 6, 4}, { 0, 0, -2, 5, -2, 5, -3, 5}, {-1, -1, -1, -1, 6, -3, 6, 6}, { 0, 1, 2, 3, 4, 5, 6, 7}
119  };
120 
121  int value = knowledge[r1][r2];
122  // Each negative case correspond to a level of knowledge
123  if (value >= 0)
124  {
125  // otbMsgDebugMacro(<<"RCC8GraphFilter: leaving GetKnowledge method: FULL");
126  return KnowledgeStateType(FULL, static_cast<RCC8ValueType>(value));
127  }
128  else if (value == -1)
129  {
130  // otbMsgDebugMacro(<<"RCC8GraphFilter: leaving GetKnowledge method: LEVEL_1");
131  return KnowledgeStateType(LEVEL_1, OTB_RCC8_DC);
132  }
133  else if (value == -2)
134  {
135  // otbMsgDebugMacro(<<"RCC8GraphFilter: leaving GetKnowledge method.LEVEL_3");
136  return KnowledgeStateType(LEVEL_3, OTB_RCC8_DC);
137  }
138  else
139  {
140  // otbMsgDebugMacro(<<"RCC8GraphFilter: leaving GetKnowledge method.NO_INFO");
141  return KnowledgeStateType(NO_INFO, OTB_RCC8_DC);
142  }
143 }
147 template <class TPolygonList, class TOutputGraph>
148 void
151 {
152  // Call a method that can be overridden by a subclass to perform
153  // some calculations prior to splitting the main computations into
154  // separate threads
155  this->BeforeThreadedGenerateData();
156 
157  // Set up the multithreaded processing
158  ThreadStruct str;
159  str.Filter = this;
160 
161  // // Initializing edges vectors per thread
162  EdgeMapType defaultEdgeMap;
163  m_EdgesPerThread = EdgeMapVectorType(this->GetNumberOfThreads(), defaultEdgeMap);
164 
165  // Setting up multithreader
166  this->GetMultiThreader()->SetNumberOfThreads(this->GetNumberOfThreads());
167  this->GetMultiThreader()->SetSingleMethod(this->ThreaderCallback, &str);
168 
169  // multithread the execution
170  this->GetMultiThreader()->SingleMethodExecute();
171 
172  // Call a method that can be overridden by a subclass to perform
173  // some calculations after all the threads have completed
174  this->AfterThreadedGenerateData();
175 }
176 
177 template <class TPolygonList, class TOutputGraph>
178 void
181 {
182  // Ouptut graph pointer
183  OutputGraphPointerType graph = this->GetOutput();
184  PolygonListConstPointerType inputPtr = this->GetInput();
185 
186 // Initializing output graph
187  graph->SetNumberOfVertices(inputPtr->Size());
188  graph->Build();
189 
190  // Vertex indexes
191  unsigned int nbVertices = 0;
192 
193  typedef typename PolygonListType::ConstIterator PolygonListConstIteratorType;
194 
195  // Loads the polygons list to graph nodes
196  for (PolygonListConstIteratorType it = inputPtr->Begin();
197  it != inputPtr->End(); ++it)
198  {
199  // Create a new vertex
200  VertexPointerType vertex = VertexType::New();
201  // Set its properties
202  vertex->SetPath(it.Get());
203 
204  // look for the appropriate segmentation index
205 
206  unsigned int segIndex = 1;
207 
208  while (segIndex < m_SegmentationRanges.size() &&
209  (nbVertices < m_SegmentationRanges[segIndex - 1] || nbVertices >= m_SegmentationRanges[segIndex]))
210  {
211  ++segIndex;
212  }
213 
214  vertex->SetSegmentationLevel(segIndex - 1);
215 
216  // Put it in the graph
217  graph->SetVertex(nbVertices, vertex);
218  otbMsgDevMacro(<< "Adding vertex: " << nbVertices);
219  ++nbVertices;
220  }
221 }
222 
223 template <class TPolygonList, class TOutputGraph>
224 void
226 ::ThreadedGenerateData(unsigned int startIndex, unsigned int stopIndex, itk::ThreadIdType threadId)
227 {
228  //std::cout<<"Starting thread "<<threadId <<" to work on range ["<<startIndex<<", "<<stopIndex<<"]"<<std::endl;
229 
230  // Ouptut graph pointer
231  OutputGraphPointerType graph = this->GetOutput();
232  PolygonListConstPointerType inputPtr = this->GetInput();
233 
234  // invert value vector
237 
238  unsigned int nbVertices = graph->GetNumberOfVertices();
239  itk::ProgressReporter progress(this, threadId, ((int) stopIndex - (int) startIndex)*nbVertices);
240 
241  otbMsgDevMacro(<< "Adjacency matrix size: " << nbVertices * nbVertices);
242 
243  VertexIteratorType vIt1(graph);
244  VertexIteratorType vIt2(graph);
245 
246  unsigned int count = 0;
247 
248  // For each couple of vertices
249  for (vIt1.GoToBegin(); !vIt1.IsAtEnd(); ++vIt1)
250  {
251  // TODO: this is not correct and should be replaced by a
252  // vIt1(graph, index1, index2), which does not exist for the moment
253  if ((count >= startIndex) && (count < stopIndex))
254  {
255  for (vIt2.GoToBegin(); !vIt2.IsAtEnd(); ++vIt2)
256  {
257  //We do not examine each couple because of the RCC8 symmetry
258  if (vIt1.GetIndex() < vIt2.GetIndex())
259  {
260 
261  // Compute the RCC8 relation
262  typename RCC8CalculatorType::Pointer calc = RCC8CalculatorType::New();
263  calc->SetPolygon1(vIt1.Get()->GetPath());
264  calc->SetPolygon2(vIt2.Get()->GetPath());
265  RCC8ValueType value = OTB_RCC8_DC;
266 
267  // if the optimisations are activated
268  if (m_Optimisation)
269  {
270  // otbMsgDebugMacro(<<"RCC8GraphFilter: Entering optimisation loop");
271  InEdgeIteratorType inIt1(vIt1.GetIndex(), graph);
272  InEdgeIteratorType inIt2(vIt2.GetIndex(), graph);
273  // otbMsgDebugMacro(<<"Optimisation loop: iterators initialised");
274  VertexDescriptorType betweenIndex;
275  KnowledgeStateType know(NO_INFO, OTB_RCC8_DC);
276  inIt1.GoToBegin();
277 
278  // Iterate through the edges going to the first vertex
279  while (!inIt1.IsAtEnd() && (know.first != FULL))
280  {
281  betweenIndex = inIt1.GetSourceIndex();
282  inIt2.GoToBegin();
283  bool edgeFound = false;
284  while (!inIt2.IsAtEnd() && (know.first != FULL))
285  {
286  // try to find an intermediate vertex between the two ones which
287  // we vant to compute the relationship
288  if (inIt2.GetSourceIndex() == betweenIndex)
289  {
290  // if an intermediate vertex is found
291  edgeFound = true;
292  // otbMsgDebugMacro(<<"Optimisation loop: found an intermediary vertex:" <<betweenIndex);
293  // See if it brings some info on the RCCC8 value
294  know = GetKnowledge(invert[inIt1.GetValue()], inIt2.GetValue());
295  calc->SetLevel1APrioriKnowledge(know.first == LEVEL_1);
296  calc->SetLevel3APrioriKnowledge(know.first == LEVEL_3);
297  // otbMsgDebugMacro(<<"Optimisation loop: knowledge: "<<know.first<<","<<know.second);
298  }
299  ++inIt2;
300  }
301  // If no intermediate was found
302  if (!edgeFound)
303  {
304  // otbMsgDebugMacro(<<"Optimisation loop: found an intermediary vertex:" <<betweenIndex);
305  // Try using a DC relationship
306  know = GetKnowledge(invert[inIt1.GetValue()], OTB_RCC8_DC);
307  calc->SetLevel1APrioriKnowledge(know.first == LEVEL_1);
308  calc->SetLevel3APrioriKnowledge(know.first == LEVEL_3);
309  // otbMsgDebugMacro(<<"Optimisation loop: knowledge: "<<know.first<<","<<know.second);
310  }
311  ++inIt1;
312  }
313  // If the search has fully determined the RCC8
314  if (know.first == FULL)
315  {
316  // Get the value
317  value = know.second;
318  }
319  else
320  {
321  // Else trigger the computation
322  // (which will take the optimisation phase info into account)
323  calc->Compute();
324  value = calc->GetValue();
325  }
326  // otbMsgDebugMacro(<<"RCC8GraphFilter: Leaving optimisation loop");
327  }
328  // If the optimisations are not activated
329  else
330  {
331  calc->Compute();
332  value = calc->GetValue();
333  }
334  m_Accumulator[value] += 1;
335  m_Accumulator[invert[value]] += 1;
336  // If the vertices are connected
337  if (value > OTB_RCC8_DC)
338  {
339  // Add the edge to the graph.
340  otbMsgDevMacro(<< "Adding edge: " << vIt1.GetIndex() << " -> " << vIt2.GetIndex() << ": " << value);
341  if (value == OTB_RCC8_NTPPI)
342  {
343  m_EdgesPerThread[threadId][EdgePairType(vIt2.GetIndex(), vIt1.GetIndex())] = OTB_RCC8_NTPP;
344  }
345  else if (value == OTB_RCC8_TPPI)
346  {
347  m_EdgesPerThread[threadId][EdgePairType(vIt2.GetIndex(), vIt1.GetIndex())] = OTB_RCC8_TPP;
348  }
349  else
350  {
351  m_EdgesPerThread[threadId][EdgePairType(vIt1.GetIndex(), vIt2.GetIndex())] = value;
352  }
353  if (m_UseInverted)
354  {
355  m_EdgesPerThread[threadId][EdgePairType(vIt2.GetIndex(), vIt1.GetIndex())] = invert[value];
356  }
357  }
358  }
359  progress.CompletedPixel();
360  progress.CompletedPixel();
361  }
362  }
363  ++count;
364  }
365  otbMsgDebugMacro(<< "End thread " << threadId);
366 }
367 
368 template <class TPolygonList, class TOutputGraph>
369 void
372 {
373  // in order to have the same output graph whatever the number of
374  // thread is, we use a map to sort the edges in lexicographical
375  // order
376 
377  OutputGraphPointerType graph = this->GetOutput();
378  EdgeMapType globalEdgeMap;
379 
380  // merge all edges
381  for (typename EdgeMapVectorType::iterator vIt = m_EdgesPerThread.begin();
382  vIt != m_EdgesPerThread.end(); ++vIt)
383  {
384  for (typename EdgeMapType::iterator mIt = (*vIt).begin();
385  mIt != (*vIt).end(); ++mIt)
386  {
387  globalEdgeMap[mIt->first] = mIt->second;
388  }
389  }
390 
391  // Report edges to the graph
392  for (typename EdgeMapType::iterator mIt = globalEdgeMap.begin();
393  mIt != globalEdgeMap.end(); ++mIt)
394  {
395  graph->AddEdge(mIt->first.first, mIt->first.second, mIt->second);
396  }
397 }
398 
399 // Callback routine used by the threading library. This routine just calls
400 // the ThreadedGenerateData method after setting the correct region for this
401 // thread.
402 
403 template <class TPolygonList, class TOutputGraph>
407 {
408  ThreadStruct *str;
409  int threadId, threadCount;
410  unsigned int total, start, stop;
411 
412  threadId = ((itk::MultiThreader::ThreadInfoStruct *) (arg))->ThreadID;
413  threadCount = ((itk::MultiThreader::ThreadInfoStruct *) (arg))->NumberOfThreads;
414  str = (ThreadStruct *) (((itk::MultiThreader::ThreadInfoStruct *) (arg))->UserData);
415 
416  total = str->Filter->GetOutput()->GetNumberOfVertices();
417 
418  if (threadId < static_cast<int>(total))
419  {
420 
421  // Split the adjacency matrix in strip of equal dimension
422  start =
423  static_cast<unsigned int>(vcl_floor(total *
424  vcl_sqrt(static_cast<double>(threadId) /
425  static_cast<double>(threadCount)) + 0.5));
426  stop =
427  static_cast<unsigned int>(vcl_floor(total *
428  vcl_sqrt(static_cast<double>(threadId +
429  1) / static_cast<double>(threadCount)) + 0.5));
430  if (stop > total) stop = total;
431 
432  // For very small graphs it might occur that start = stop. In this
433  // case the vertex at that index will be processed in the next strip.
434  if (start != stop)
435  {
436  str->Filter->ThreadedGenerateData(start, stop, threadId);
437  }
438  }
439  // else
440  // {
441  // otherwise don't use this thread. Sometimes the threads dont
442  // break up very well and it is just as efficient to leave a
443  // few threads idle.
444  // }
445 
447 }
448 
449 template <class TPolygonList, class TOutputGraph>
450 void
452 ::PrintSelf(std::ostream& os, itk::Indent indent) const
453 {
454  Superclass::PrintSelf(os, indent);
455 }
456 } // end namespace otb
457 #endif

Generated at Sat Mar 8 2014 16:13:25 for Orfeo Toolbox with doxygen 1.8.3.1