Skip to content
GitLab
Select Git revision
  • master default
  • enhancements
Find file BlameHistoryPermalink
otbZonalStatistics.cxx 12.22 KiB
1 
/*=========================================================================
2
3 
  Copyright (c) Remi Cresson (IRSTEA). All rights reserved.
4
5
6 
     This software is distributed WITHOUT ANY WARRANTY; without even
7 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
8 
     PURPOSE.  See the above copyright notices for more information.
9
10 
=========================================================================*/
11 
#include "itkFixedArray.h"
12 
#include "itkObjectFactory.h"
13
14 
// Elevation handler
15 
#include "otbWrapperElevationParametersHandler.h"
16 
#include "otbWrapperApplicationFactory.h"
17
18 
// Application engine
19 
#include "otbStandardFilterWatcher.h"
20 
#include "itkFixedArray.h"
21
22 
// Filter
23 
#include "otbVectorDataToLabelImageFilter.h"
24 
#include "otbVectorDataIntoImageProjectionFilter.h"
25
26 
// Vnl vector
27 
#include "vnl/vnl_vector.h"
28
29 
// itk iterator
30 
#include "itkImageRegionConstIterator.h"
31
32 
using namespace std;
33
34 
namespace otb
35 
{
36
37 
namespace Wrapper
38 
{
39
40 
class ZonalStatistics : public Application
41 
{
42 
public:
43 
  /** Standard class typedefs. */
44 
  typedef ZonalStatistics                        Self;
45 
  typedef Application                   Superclass;
46 
  typedef itk::SmartPointer<Self>       Pointer;
47 
  typedef itk::SmartPointer<const Self> ConstPointer;
48
49
50 
  /* vector data filters */
51 
  typedef size_t                                                      LabelValueType;
52 
  typedef otb::Image<LabelValueType, 2>                               LabelImageType;
53 
  typedef otb::VectorData<double, 2>                                  VectorDataType;
54 
  typedef otb::VectorDataIntoImageProjectionFilter<VectorDataType,
55 
      FloatVectorImageType>                                           VectorDataReprojFilterType;
56 
  typedef otb::VectorDataToLabelImageFilter<VectorDataType,
57 
      LabelImageType>                                                 RasterizeFilterType;
58 
  typedef VectorDataType::DataTreeType                                DataTreeType;
59 
  typedef itk::PreOrderTreeIterator<DataTreeType>                     TreeIteratorType;
60 
  typedef VectorDataType::DataNodeType                                DataNodeType;
61 
  typedef DataNodeType::PolygonListPointerType                        PolygonListPointerType;
62
63 
  /** Statistics */
64 
  typedef vnl_matrix<double>                                          RealMatrix;
65 
  typedef vnl_vector<LabelValueType>                                  SizeVector;
66
67 
  /** Image iterator */
68 
  typedef itk::ImageRegionConstIterator<LabelImageType>               LabelIteratorType;
69 
  typedef itk::ImageRegionConstIterator<FloatVectorImageType>         ImageIteratorType;
70
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
/** Standard macro */ itkNewMacro(Self); itkTypeMacro(ZonalStatistics, Application); void DoInit() { SetName("ZonalStatistics"); SetDescription("Computes zonal statistics"); // Documentation SetDocName("ZonalStatistics"); SetDocLongDescription("This application zonal statistics"); SetDocLimitations("None"); SetDocAuthors("Remi Cresson"); AddDocTag(Tags::Manip); // Inputs AddParameter(ParameterType_InputImage, "in", "Input Image"); AddParameter(ParameterType_InputVectorData, "vec", "Input vector data"); // No data value AddParameter(ParameterType_Float, "nodata", "No-data value"); MandatoryOff("nodata"); // Output AddParameter(ParameterType_OutputVectorData, "out", "Output vector data"); AddRAMParameter(); // Doc example parameter settings SetDocExampleParameterValue("in", "input.tif"); SetDocExampleParameterValue("vec", "terrain.shp"); SetDocExampleParameterValue("out", "terain_with_stats.shp"); } void DoUpdateParameters() { // Nothing to do here : all parameters are independent } // Create a string of the kind "prefix_i" std::string CreateFieldName(std::string prefix, unsigned int i) { std::stringstream ss; ss << prefix << "_" << i; return ss.str(); } void DoExecute() { // Get inputs FloatVectorImageType::Pointer img = GetParameterImage("in"); VectorDataType* shp = GetParameterVectorData("vec"); const bool has_nodata = HasValue("nodata"); FloatVectorImageType::InternalPixelType m_InputNoData = 0; if (has_nodata) { m_InputNoData = GetParameterFloat("nodata"); otbAppLogINFO("Using user no-data value: " << m_InputNoData); } // Reproject vector data otbAppLogINFO("Reproject vector data"); m_VectorDataReprojectionFilter = VectorDataReprojFilterType::New(); m_VectorDataReprojectionFilter->SetInputVectorData(shp);
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
m_VectorDataReprojectionFilter->SetInputImage(img); m_VectorDataReprojectionFilter->Update(); // Add a FID field otbAppLogINFO("Adding internal FID"); LabelValueType internalFID = 1; const string internalFIDField = "fid____"; vdWithFid = VectorDataType::New(); DataNodeType::Pointer root1 = vdWithFid->GetDataTree()->GetRoot()->Get(); DataNodeType::Pointer document1 = DataNodeType::New(); document1->SetNodeType(otb::DOCUMENT); vdWithFid->GetDataTree()->Add(document1, root1); DataNodeType::Pointer folder1 = DataNodeType::New(); folder1->SetNodeType(otb::FOLDER); vdWithFid->GetDataTree()->Add(folder1, document1); vdWithFid->SetProjectionRef(m_VectorDataReprojectionFilter->GetOutput()->GetProjectionRef()); TreeIteratorType itVector1(m_VectorDataReprojectionFilter->GetOutput()->GetDataTree()); itVector1.GoToBegin(); while (!itVector1.IsAtEnd()) { if (!itVector1.Get()->IsRoot() && !itVector1.Get()->IsDocument() && !itVector1.Get()->IsFolder()) { DataNodeType::Pointer currentGeometry = itVector1.Get(); currentGeometry->SetFieldAsInt(internalFIDField, internalFID ); vdWithFid->GetDataTree()->Add(currentGeometry, folder1); internalFID++; } ++itVector1; } // next feature // Rasterize vector data otbAppLogINFO("Rasterize vector data"); m_RasterizeFilter = RasterizeFilterType::New(); m_RasterizeFilter->AddVectorData(vdWithFid); m_RasterizeFilter->SetOutputOrigin(img->GetOrigin()); m_RasterizeFilter->SetOutputSpacing(img->GetSignedSpacing()); m_RasterizeFilter->SetOutputSize(img->GetLargestPossibleRegion().GetSize()); m_RasterizeFilter->SetOutputProjectionRef(img->GetProjectionRef()); m_RasterizeFilter->SetBurnAttribute(internalFIDField); m_RasterizeFilter->Update(); // Computing stats otbAppLogINFO("Computing stats"); // Explicit streaming over the input target image, based on the RAM parameter typedef otb::RAMDrivenStrippedStreamingManager<FloatVectorImageType> StreamingManagerType; StreamingManagerType::Pointer m_StreamingManager = StreamingManagerType::New(); m_StreamingManager->SetAvailableRAMInMB(GetParameterInt("ram")); m_StreamingManager->PrepareStreaming(img, img->GetLargestPossibleRegion() ); // Init. stats containers const LabelValueType N = internalFID; const unsigned int nBands = img->GetNumberOfComponentsPerPixel(); RealMatrix minimums(nBands, N, itk::NumericTraits<double>::max()); RealMatrix maximums(nBands, N, itk::NumericTraits<double>::min()); RealMatrix sum (nBands, N, 0); RealMatrix sqSum (nBands, N, 0); RealMatrix means (nBands, N, 0); RealMatrix stdevs (nBands, N, 0); SizeVector counts(N,0); if (m_RasterizeFilter->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels() != img->GetLargestPossibleRegion().GetNumberOfPixels()) { otbAppLogFATAL("Rasterized image and input image don't have the same number of pixels"); } int m_NumberOfDivisions = m_StreamingManager->GetNumberOfSplits(); for (int m_CurrentDivision = 0; m_CurrentDivision < m_NumberOfDivisions; m_CurrentDivision++) { std::cout << "region: " << m_CurrentDivision << std::endl; FloatVectorImageType::RegionType streamRegion = m_StreamingManager->GetSplit(m_CurrentDivision);
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
img->SetRequestedRegion(streamRegion); img->PropagateRequestedRegion(); img->UpdateOutputData(); LabelIteratorType it(m_RasterizeFilter->GetOutput(), streamRegion); ImageIteratorType it_in(img, streamRegion); for (it.GoToBegin(), it_in.GoToBegin(); !it.IsAtEnd(); ++it, ++it_in) { const LabelValueType fid = it.Get(); const FloatVectorImageType::PixelType pix = it_in.Get(); bool valid = true; if (has_nodata) { for (unsigned int band = 0 ; band < nBands ; band++) { const FloatVectorImageType::InternalPixelType val = pix[band]; if (val == m_InputNoData) { valid = false; break; } } } if (valid) { for (unsigned int band = 0 ; band < nBands ; band++) { const FloatVectorImageType::InternalPixelType val = pix[band]; if (val < minimums[band][fid]) minimums[band][fid] = val; if (val > maximums[band][fid]) maximums[band][fid] = val; sum[band][fid]+=val; sqSum[band][fid]+=val*val; } // next band counts[fid]++; } } // next pixel } // Summarize stats otbAppLogINFO("Summarizing stats"); for (LabelValueType fid = 0 ; fid < N ; fid++) { const LabelValueType count = counts[fid]; for (unsigned int band = 0 ; band < nBands ; band++) { if (count>0) { means[band][fid] = sum[band][fid] / static_cast<FloatVectorImageType::InternalPixelType>(count); // Unbiased estimate FloatVectorImageType::InternalPixelType variance = 0; if (count > 1) { variance = (sqSum[band][fid] - (sum[band][fid]*sum[band][fid] / static_cast<FloatVectorImageType::InternalPixelType>(count) ) ) / (static_cast<FloatVectorImageType::InternalPixelType>(count) - 1); if (variance > 0) { stdevs[band][fid] = vcl_sqrt(variance); } } } } } // Add a statistics fields otbAppLogINFO("Writing stats in output vector data"); internalFID = 1;
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330
vdWithStats = VectorDataType::New(); DataNodeType::Pointer root = vdWithStats->GetDataTree()->GetRoot()->Get(); DataNodeType::Pointer document = DataNodeType::New(); document->SetNodeType(otb::DOCUMENT); vdWithStats->GetDataTree()->Add(document, root); DataNodeType::Pointer folder = DataNodeType::New(); folder->SetNodeType(otb::FOLDER); vdWithStats->GetDataTree()->Add(folder, document); vdWithStats->SetProjectionRef(m_VectorDataReprojectionFilter->GetOutput()->GetProjectionRef()); TreeIteratorType itVector(m_VectorDataReprojectionFilter->GetOutput()->GetDataTree()); itVector.GoToBegin(); while (!itVector.IsAtEnd()) { if (!itVector.Get()->IsRoot() && !itVector.Get()->IsDocument() && !itVector.Get()->IsFolder()) { DataNodeType::Pointer currentGeometry = itVector.Get(); currentGeometry->SetFieldAsInt("fid", internalFID ); for (unsigned int band = 0 ; band < nBands ; band++) { currentGeometry->SetFieldAsDouble(CreateFieldName("mean", band), means [band][internalFID] ); currentGeometry->SetFieldAsDouble(CreateFieldName("stdev", band), stdevs [band][internalFID] ); currentGeometry->SetFieldAsDouble(CreateFieldName("min", band), minimums[band][internalFID] ); currentGeometry->SetFieldAsDouble(CreateFieldName("max", band), maximums[band][internalFID] ); // currentGeometry->SetFieldAsDouble(CreateFieldName("sum", band), sum [band][internalFID] ); // currentGeometry->SetFieldAsDouble(CreateFieldName("sqSum", band), sqSum [band][internalFID] ); } currentGeometry->SetFieldAsDouble("count", static_cast<double>(counts[internalFID]) ); vdWithStats->GetDataTree()->Add(currentGeometry, folder); internalFID++; } ++itVector; } // next feature SetParameterOutputVectorData("out", vdWithStats); } VectorDataType::Pointer vdWithFid; VectorDataType::Pointer vdWithStats; VectorDataReprojFilterType::Pointer m_VectorDataReprojectionFilter; RasterizeFilterType::Pointer m_RasterizeFilter; }; } } OTB_APPLICATION_EXPORT( otb::Wrapper::ZonalStatistics )