OTB  5.0.0
Orfeo Toolbox
otbVectorDataToImageFilter.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 
19 #ifndef __otbVectorDataToImageFilter_txx
20 #define __otbVectorDataToImageFilter_txx
21 
22 #include <sstream>
24 #include "itkImageRegionIterator.h"
26 #include "otbVectorDataStyle.h"
27 #include "itkRGBAPixel.h"
28 
29 #include "ogr_spatialref.h"
30 
31 #include "otbMapnikAdapter.h"
32 
33 namespace otb
34 {
35 
39 template <class TVectorData, class TImage>
42  m_StyleList(),
43  m_UseAsOverlay(true),
44  m_RenderingStyleType(OSM)
45 {
46  this->SetNumberOfRequiredInputs(1);
47  m_Spacing.Fill(1.0);
48  m_Origin.Fill(0.0);
49  m_Direction.SetIdentity();
50  m_Size.Fill(0);
51  m_StartIndex.Fill(0);
52  m_SensorModelFlip = 1;
53  m_ScaleFactor = 1.0;
54  m_VectorDataProjectionProj4 = "";
55  m_VectorDataProjectionWKT = "";
56 }
58 
59 template <class TVectorData, class TImage>
60 void
63 {
64  // Process object is not const-correct so the const_cast is required here
66  const_cast<VectorDataType *>(input));
67 }
68 
69 template <class TVectorData, class TImage>
70 void
72 ::SetInput(unsigned int idx, const VectorDataType *input)
73 {
74  // Process object is not const-correct so the const_cast is required here
76  const_cast<VectorDataType *>(input));
77 }
78 
79 template <class TVectorData, class TImage>
83 {
84  if (this->GetNumberOfInputs() < 1)
85  {
86  return 0;
87  }
88 
89  return static_cast<const TVectorData *>
91 }
92 
93 template <class TVectorData, class TImage>
96 ::GetInput(unsigned int idx)
97 {
98  return static_cast<const TVectorData *>
99  (this->itk::ProcessObject::GetInput(idx));
100 }
101 
102 //----------------------------------------------------------------------------
103 template <class TVectorData, class TImage>
104 void
106 ::SetSpacing(const SpacingType& spacing)
107 {
108  if (this->m_Spacing != spacing)
109  {
110  this->m_Spacing = spacing;
111  this->Modified();
112  }
113 }
114 
115 //----------------------------------------------------------------------------
116 template <class TVectorData, class TImage>
117 void
119 ::SetSpacing(const double spacing[2])
120 {
121  SpacingType s(spacing);
122  this->SetSpacing(s);
123 }
124 
125 //----------------------------------------------------------------------------
126 template <class TVectorData, class TImage>
127 void
129 ::SetSpacing(const float spacing[2])
130 {
131  itk::Vector<float, 2> sf(spacing);
132  SpacingType s;
133  s.CastFrom(sf);
134  this->SetSpacing(s);
135 }
136 
137 //----------------------------------------------------------------------------
138 template <class TVectorData, class TImage>
139 void
141 ::SetOrigin(const double origin[2])
142 {
143  OriginType p(origin);
144  this->SetOrigin(p);
145 }
146 
147 //----------------------------------------------------------------------------
148 template <class TVectorData, class TImage>
149 void
151 ::SetOrigin(const float origin[2])
152 {
153  itk::Point<float, 2> of(origin);
154  OriginType p;
155  p.CastFrom(of);
156  this->SetOrigin(p);
157 }
158 
162 template <class TVectorData, class TImage>
163 void
166 {
167  // we can't call the superclass method here.
168 
169  // get pointers to the input and output
170  ImagePointer outputPtr = this->GetOutput();
171  if (!outputPtr)
172  {
173  return;
174  }
175 
176  // Set the size of the output region
177  typename TImage::RegionType outputLargestPossibleRegion;
178  outputLargestPossibleRegion.SetSize(m_Size);
179  outputLargestPossibleRegion.SetIndex(m_StartIndex);
180  outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
181 
182  // Set spacing and origin
183  outputPtr->SetSpacing(m_Spacing);
184  outputPtr->SetOrigin(m_Origin);
185  outputPtr->SetDirection(m_Direction);
186 
187  itk::MetaDataDictionary& dict = outputPtr->GetMetaDataDictionary();
188  itk::EncapsulateMetaData<std::string> (dict, MetaDataKey::ProjectionRefKey,
189  static_cast<std::string>(m_VectorDataProjectionWKT));
190 
191  //TODO update or check the projection information
192 
193  return;
194 }
195 
199 template <class TVectorData, class TImage>
200 void
203 {
204 
205  Superclass::BeforeThreadedGenerateData();
206 
207  //Font Handling
208  if(!m_FontFileName.empty())
209  mapnik::freetype_engine::register_font(m_FontFileName);
210 
211  //Handle the style type using helper class
213  styleLoader->SetScaleFactor(m_ScaleFactor);
214 
215  //We assume that all the data are reprojected before using OTB.
216  VectorDataConstPointer input = this->GetInput();
217  //Converting the projection string to the proj.4 format
218  itk::ExposeMetaData<std::string>(
219  input->GetMetaDataDictionary(), MetaDataKey::ProjectionRefKey, m_VectorDataProjectionWKT);
220  otbMsgDebugMacro(<< "WKT -> " << m_VectorDataProjectionWKT);
221 
222  m_SensorModelFlip = 1;
223 
224  if (m_VectorDataProjectionWKT == "")
225  {
226  //We assume that it is an image in sensor model geometry
227  //and tell mapnik that this is utm
228  //(with a resolution of 1m per unit)
229  m_VectorDataProjectionProj4 = "+proj=utm +zone=31 +ellps=WGS84";
230  m_SensorModelFlip = -1;
231  otbMsgDevMacro(<< "The output map will be in sensor geometry");
232  }
233  else
234  {
235  OGRSpatialReferenceH oSRS = OSRNewSpatialReference(m_VectorDataProjectionWKT.c_str());
236  char * pszProj4;
237  OSRExportToProj4(oSRS, &pszProj4);
238  m_VectorDataProjectionProj4 = pszProj4;
239  CPLFree(pszProj4);
240  OSRRelease(oSRS);
241  m_SensorModelFlip = 1;
242  otbMsgDevMacro(<< "The output map will be carto/geo geometry");
243  }
244  otbMsgDebugMacro(<< "Proj.4 -> " << m_VectorDataProjectionProj4);
245 
246  //Internal Tiling Support
247  //Workaround to overcome mapnik maximum tile size limitation
248  ImagePointer output = this->GetOutput();
249  RegionType requestedRegion = output->GetRequestedRegion();
250  otbMsgDevMacro("requestedRegion: " << requestedRegion);
251 
252  m_NbTile = (vcl_pow(std::max(vcl_floor((double)requestedRegion.GetSize()[0] / 16000),
253  vcl_floor((double)requestedRegion.GetSize()[1] / 16000))+1, 2));
254 
255  //std::cout << "nbTile: " << m_NbTile << std::endl;
256  //std::cout << "requestedRegion: " << requestedRegion << std::endl;
257 
258  m_TilingRegions.resize(m_NbTile);
259  m_Maps.resize(m_NbTile);
260  m_VectorDataExtractors.resize(m_NbTile);
261 
262  unsigned int tilingRegionsIdx = 0;
263  unsigned int stdXOffset;
264  unsigned int stdYOffset;
265 
266  stdXOffset = vcl_floor((double)requestedRegion.GetSize()[0]/ (m_NbTile/2))+1;
267  stdYOffset = vcl_floor((double)requestedRegion.GetSize()[1]/ (m_NbTile/2))+1;
268 
269  for(unsigned int i=0; i < vcl_floor((double)(m_NbTile)/2 + 0.5); ++i)
270  {
271  for(unsigned int j=0; j < vcl_floor((double)(m_NbTile)/2 + 0.5); ++j)
272  {
273  //Set Regions
274  SizeType size;
275  IndexType index;
276  if(m_NbTile == 1)
277  {
278  index = requestedRegion.GetIndex();
279  size = requestedRegion.GetSize();
280  }
281  else
282  {
283  index[0] = requestedRegion.GetIndex()[0] + i * stdXOffset;
284  index[1] = requestedRegion.GetIndex()[1] + j * stdYOffset;
285 
286  size[0] = std::min((unsigned int)(requestedRegion.GetSize()[0] - index[0]), stdXOffset);
287  size[1] = std::min((unsigned int)(requestedRegion.GetSize()[1] - index[1]), stdYOffset);
288  }
289  m_TilingRegions[tilingRegionsIdx].SetIndex(index);
290  m_TilingRegions[tilingRegionsIdx].SetSize(size);
291 
292  //std::cout << "tileRegions[" << tilingRegionsIdx << "] : "
293  // << m_TilingRegions[tilingRegionsIdx] << std::endl;
294 
295  //Set Maps
296  m_Maps[tilingRegionsIdx] = mapnik::Map();
297 
299  switch (m_RenderingStyleType)
300  {
301  case OSM:
302  {
303  styleLoader->LoadOSMStyle(m_Maps[tilingRegionsIdx]);
304  if (m_UseAsOverlay)
305  {
306  //Set the default backgroup to transparent
307  m_Maps[tilingRegionsIdx].set_background(mapnik::color(255, 255, 255, 0));
308  }
309  else
310  {
311  m_Maps[tilingRegionsIdx].set_background(mapnik::color("#b5d0d0"));
312  }
313  break;
314  }
315  case Binary:
316  {
317  styleLoader->LoadBinaryRasterizationStyle(m_Maps[tilingRegionsIdx]);
318  //Set the backgroup to white
319  m_Maps[tilingRegionsIdx].set_background(mapnik::color("#ffffff"));
320  break;
321  }
322  default:
323  {
324  itkExceptionMacro(<< "Style Type Not Supported!");
325  break;
326  }
327  }
328 
330  m_Maps[tilingRegionsIdx].set_srs(m_VectorDataProjectionProj4);
331 
332  //Set VectorData extracts
333  m_VectorDataExtractors[tilingRegionsIdx].resize(this->GetNumberOfInputs());
334  for (unsigned int idx = 0; idx < this->GetNumberOfInputs(); ++idx)
335  {
336  if (this->GetInput(idx))
337  {
338  RemoteSensingRegionType rsRegion;
339  SizePhyType sizePhy;
340  sizePhy[0] = size[0] * m_Spacing[0];
341  sizePhy[1] = size[1] * m_Spacing[1];
342  rsRegion.SetSize(sizePhy);
343  OriginType origin;
344  origin[0] = m_Origin[0] + (static_cast<double>(index[0]) - 0.5) * m_Spacing[0];
345  origin[1] = m_Origin[1] + (static_cast<double>(index[1]) - 0.5) * m_Spacing[1];
346  rsRegion.SetOrigin(origin);
347  rsRegion.SetRegionProjection(m_VectorDataProjectionWKT);
348 
349  //std::cout << "m_SensorModelFlip: " << m_SensorModelFlip << std::endl;
350  //std::cout << "rsTileRegions[" << tilingRegionsIdx << "] : "
351  // << rsRegion << std::endl;
352 
353  m_VectorDataExtractors[tilingRegionsIdx][idx] = VectorDataExtractROIType::New();
354  m_VectorDataExtractors[tilingRegionsIdx][idx]->SetRegion(rsRegion);
355  m_VectorDataExtractors[tilingRegionsIdx][idx]->SetInput(this->GetInput(idx));
356  }
357  }
358 
359  tilingRegionsIdx ++;
360  }
361  }
362 }
363 
367 template <class TVectorData, class TImage>
368 void
371 {
372 
373  this->AllocateOutputs();
374 
375  this->BeforeThreadedGenerateData();
376 
377  if (m_StyleList.size() == 0)
378  {
379  switch (m_RenderingStyleType)
380  {
381  case OSM:
382  {
383  //Add default styles
384  itkExceptionMacro(<< "No style is provided for the vector data");
385  break;
386  }
387  case Binary:
388  {
389  //Add default styles
390  this->AddStyle("binary-rasterization");
391  break;
392  }
393  default:
394  {
395  itkExceptionMacro(<< "No style is provided for the vector data");
396  break;
397  }
398  }
399  }
400 
401  ImagePointer output = this->GetOutput();
402 
403  for (unsigned int tileIdx = 0; tileIdx < m_NbTile; ++tileIdx)
404  {
405  //Delete the previous layers from the map
406  int numberLayer = mapnik_otb::get_num_layer(m_Maps[tileIdx]);
407  for (int i = numberLayer - 1; i >= 0; i--) //yes, int.
408  {
409  m_Maps[tileIdx].removeLayer(i);
410  }
411  m_Maps[tileIdx].resize(m_TilingRegions[tileIdx].GetSize()[0], m_TilingRegions[tileIdx].GetSize()[1]);
412 
413  RemoteSensingRegionType rsRegion;
414 
415  for (unsigned int vdIdx = 0; vdIdx < this->GetNumberOfInputs(); ++vdIdx)
416  {
417  if (this->GetInput(vdIdx))
418  {
419  datasource_ptr mDatasource = datasource_ptr(new mapnik::memory_datasource);
420 
421  rsRegion = m_VectorDataExtractors[tileIdx][vdIdx]->GetRegion();
422 
423  m_VectorDataExtractors[tileIdx][vdIdx]->Update();
424  VectorDataConstPointer input = m_VectorDataExtractors[tileIdx][vdIdx]->GetOutput();
425  InternalTreeNodeType * inputRoot = const_cast<InternalTreeNodeType *>(input->GetDataTree()->GetRoot());
426 
427  ProcessNode(inputRoot, mDatasource);
428  otbMsgDevMacro("Datasource size: " << mDatasource->size());
429 
430  std::stringstream layerName;
431  layerName << "layer-" << tileIdx;
432  mapnik::layer lyr(layerName.str());
433  lyr.set_srs(m_VectorDataProjectionProj4);
434  lyr.set_datasource(mDatasource);
435 
436  for (unsigned int i = 0; i < m_StyleList.size(); ++i)
437  {
438  lyr.add_style(m_StyleList[i]);
439  }
440 
441  m_Maps[tileIdx].addLayer(lyr);
442  }
443  }
444  assert((m_SensorModelFlip == 1) || (m_SensorModelFlip == -1));
445 
446  mapnik_otb::box2d envelope(
447  rsRegion.GetOrigin(0),
448  m_SensorModelFlip*(rsRegion.GetOrigin(1) + rsRegion.GetSize(1)),
449  rsRegion.GetOrigin(0) + rsRegion.GetSize(0),
450  m_SensorModelFlip*(rsRegion.GetOrigin(1))
451  );
452 
453  mapnik_otb::zoom_to_box(&m_Maps[tileIdx], envelope);
454  otbMsgDebugMacro(<< "Envelope: " << envelope);
455 
456  otbMsgDebugMacro(<< "Map scale: " << m_Maps[tileIdx].scale_denominator());
457  mapnik::image_32 buf(mapnik_otb::get_width(m_Maps[tileIdx]),
458  mapnik_otb::get_height(m_Maps[tileIdx]));
459  mapnik::agg_renderer<mapnik::image_32> ren(m_Maps[tileIdx], buf);
460  ren.apply();
461 
462  const unsigned char * src = buf.raw_data();
463 
464  itk::ImageRegionIterator<ImageType> it(output, m_TilingRegions[tileIdx]);
465 
466  for (it.GoToBegin(); !it.IsAtEnd(); ++it)
467  {
469  pix[0] = *src;
470  pix[1] = *(src+1);
471  pix[2] = *(src+2);
472  pix[3] = *(src+3);
473  src += 4;
474 
475  it.Set(m_RGBAConverter->Convert(pix));
476  }
477  }
478 }
479 
480 template <class TVectorData, class TImage>
481 void
484 {
485  typedef typename VectorDataType::DataNodeType DataNodeType;
486  typedef typename DataNodeType::Pointer DataNodePointerType;
487 
488  // Get the children list from the input node
489  ChildrenListType children = source->GetChildrenList();
490 
491  // For each child
492  for (typename ChildrenListType::iterator it = children.begin(); it != children.end(); ++it)
493  {
494  // Copy input DataNode info
495  DataNodePointerType dataNode = (*it)->Get();
496 
497  switch (dataNode->GetNodeType())
498  {
499  case otb::ROOT:
500  {
501  ProcessNode((*it), mDatasource);
502  break;
503  }
504  case otb::DOCUMENT:
505  {
506  ProcessNode((*it), mDatasource);
507  break;
508  }
509  case otb::FOLDER:
510  {
511  ProcessNode((*it), mDatasource);
512  break;
513  }
514  case FEATURE_POINT:
515  {
516  mapnik_otb::geom* point = mapnik_otb::create_geom(mapnik::Point);
517 
518  point->move_to(dataNode->GetPoint()[0], m_SensorModelFlip * dataNode->GetPoint()[1]);
519 // std::cout << dataNode->GetPoint()[0] << ", " << dataNode->GetPoint()[1] << std::endl;
520 
521  typedef boost::shared_ptr<mapnik::raster> raster_ptr;
522  typedef mapnik::feature<mapnik_otb::geom, raster_ptr> Feature;
523  typedef boost::shared_ptr<Feature> feature_ptr;
524 
525  feature_ptr mfeature = feature_ptr(new Feature(1));
526  mfeature->add_geometry(point);
527 
528  mapnik::transcoder tr("ISO-8859-15");
529 
530  if (dataNode->HasField("place_name"))
531  boost::put(*mfeature, "name",
532  tr.transcode((dataNode->GetFieldAsString("place_name")).c_str()));
533 
534  boost::put(*mfeature, "place", tr.transcode("city"));
535  boost::put(*mfeature, "capital", tr.transcode("yes")); //FIXME more a question of style
536 
537  boost::put(*mfeature, "geometry", tr.transcode("point"));
538 
539  mDatasource->push(mfeature);
540 
541  break;
542  }
543  case otb::FEATURE_LINE:
544  {
545  mapnik_otb::geom* line = mapnik_otb::create_geom(mapnik::LineString);
546 
547  typedef typename DataNodeType::LineType::VertexListConstIteratorType VertexIterator;
548  VertexIterator itVertex = dataNode->GetLine()->GetVertexList()->Begin();
549  while (itVertex != dataNode->GetLine()->GetVertexList()->End())
550  {
551 // std::cout << itVertex.Value()[0] << ", " << itVertex.Value()[1] << std::endl;
552  line->line_to(itVertex.Value()[0], m_SensorModelFlip * itVertex.Value()[1]);
553  ++itVertex;
554  }
555 
556 // std::cout << "Num points: " << line->num_points() << std::endl;
557 
558  typedef boost::shared_ptr<mapnik::raster> raster_ptr;
559  typedef mapnik::feature<mapnik_otb::geom, raster_ptr> Feature;
560  typedef boost::shared_ptr<Feature> feature_ptr;
561 
562  feature_ptr mfeature = feature_ptr(new Feature(1));
563  mfeature->add_geometry(line);
564 
565  mapnik::transcoder tr("ISO-8859-15");
566 
567  if (dataNode->HasField("name"))
568  boost::put(*mfeature, "name",
569  tr.transcode((dataNode->GetFieldAsString("name")).c_str()));
570  if (dataNode->HasField("NAME"))
571  boost::put(*mfeature, "name",
572  tr.transcode((dataNode->GetFieldAsString("NAME")).c_str()));
573 
574 // std::cout << mfeature->props().size() << std::endl;
575 // std::cout << " -> " << (*mfeature)["name"] << std::endl;
576 
577 // std::cout << "Name: " << dataNode->GetFieldAsString("NAME") << std::endl;
578 // std::cout << "Type: " << dataNode->GetFieldAsString("TYPE") << std::endl;
579 // std::cout << "OSM ID: " << dataNode->GetFieldAsString("osm_id") << std::endl;
580 
581  if (dataNode->HasField("type"))
582  boost::put(*mfeature, "highway",
583  tr.transcode((dataNode->GetFieldAsString("type")).c_str()));
584  if (dataNode->HasField("TYPE"))
585  boost::put(*mfeature, "highway",
586  tr.transcode((dataNode->GetFieldAsString("TYPE")).c_str()));
587 
588  boost::put(*mfeature, "geometry", tr.transcode("line"));
589 
590  mDatasource->push(mfeature);
591 
592  break;
593  }
594  case FEATURE_POLYGON:
595  {
596  mapnik_otb::geom* polygon = mapnik_otb::create_geom(mapnik::Polygon);
597 
598  typedef typename DataNodeType::PolygonType::VertexListConstIteratorType VertexIterator;
599  VertexIterator itVertex = dataNode->GetPolygonExteriorRing()->GetVertexList()->Begin();
600  while (itVertex != dataNode->GetPolygonExteriorRing()->GetVertexList()->End())
601  {
602  polygon->line_to(itVertex.Value()[0], m_SensorModelFlip * itVertex.Value()[1]);
603  ++itVertex;
604  }
605 
606  typedef boost::shared_ptr<mapnik::raster> raster_ptr;
607  typedef mapnik::feature<mapnik_otb::geom, raster_ptr> Feature;
608  typedef boost::shared_ptr<Feature> feature_ptr;
609 
610  feature_ptr mfeature = feature_ptr(new Feature(1));
611  mfeature->add_geometry(polygon);
612 
613  mapnik::transcoder tr("ISO-8859-15");
614 
615  boost::put(*mfeature, "geometry", tr.transcode("polygon"));
616 
617  mDatasource->push(mfeature);
618 
619  break;
620  }
621  case FEATURE_MULTIPOINT:
622  {
623  itkExceptionMacro(
624  << "This type (FEATURE_MULTIPOINT) is not handle (yet) by VectorDataToImageFilter(), please request for it");
625  break;
626  }
627  case FEATURE_MULTILINE:
628  {
629  itkExceptionMacro(
630  << "This type (FEATURE_MULTILINE) is not handle (yet) by VectorDataToImageFilter(), please request for it");
631  break;
632  }
634  {
635  itkExceptionMacro(
636  << "This type (FEATURE_MULTIPOLYGON) is not handle (yet) by VectorDataToImageFilter(), please request for it");
637  break;
638  }
639  case FEATURE_COLLECTION:
640  {
641  itkExceptionMacro(
642  << "This type (FEATURE_COLLECTION) is not handle (yet) by VectorDataToImageFilter(), please request for it");
643  break;
644  }
645  }
646  }
647 }
648 
652 template <class TVectorData, class TImage>
653 void
655 ::PrintSelf(std::ostream& os, itk::Indent indent) const
656 {
657  Superclass::PrintSelf(os, indent);
658 }
659 }
661 
662 #endif
InternalTreeNodeType::ChildrenListType ChildrenListType
virtual void PrintSelf(std::ostream &os, itk::Indent indent) const
virtual void SetOrigin(OriginType _arg)
void Set(const PixelType &value) const
virtual void SetSpacing(const SpacingType &spacing)
static void zoom_to_box(mapnik::Map *map, const mapnik::box2d< double > &envelope)
VectorDataType::DataTreeType::TreeNodeType InternalTreeNodeType
VectorDataType::ConstPointer VectorDataConstPointer
void ProcessNode(InternalTreeNodeType *source, datasource_ptr mDatasource)
virtual void SetInput(const VectorDataType *input)
#define otbMsgDebugMacro(x)
Definition: otbMacro.h:55
const IndexType & GetOrigin() const
void SetOrigin(const IndexType &index)
static vcl_size_t get_num_layer(const mapnik::Map &map)
static Pointer New()
static unsigned get_height(const mapnik::Map &map)
const VectorDataType * GetInput(void)
mapnik::geometry< vertex2d > geom
DataObject * GetInput(const DataObjectIdentifierType &key)
boost::shared_ptr< mapnik::memory_datasource > datasource_ptr
static unsigned get_width(const mapnik::Map &map)
void SetRegionProjection(const std::string &projection)
mapnik::box2d< double > box2d
void SetSize(const SizeType &size)
virtual void SetNthInput(DataObjectPointerArraySizeType num, DataObject *input)
bool IsAtEnd(void) const
char const * ProjectionRefKey
static geom * create_geom(mapnik::eGeomType geom_type)
const SizeType & GetSize() const
const SizeValueType * GetSize() const
#define otbMsgDevMacro(x)
Definition: otbMacro.h:95