otbComputePolylineFeatureFromImage.cxx 9.89 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
/*
 * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
 *
 * This file is part of Orfeo Toolbox
 *
 *     https://www.orfeo-toolbox.org/
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
20 21 22 23

#include <iostream>

#include "otbVectorImage.h"
24 25 26

#include "otbWrapperApplicationFactory.h"

27 28 29 30 31 32 33 34 35 36 37 38

#include "otbVectorDataFileReader.h"
#include "otbVectorDataFileWriter.h"
#include "itkPreOrderTreeIterator.h"

#include "otbVectorDataIntoImageProjectionFilter.h"
#include "otbVectorDataTransformFilter.h"
#include "itkAffineTransform.h"

#include "otbParserConditionDataNodeFeatureFunction.h"
#include "otbNDVIDataNodeFeatureFunction.h"

39 40 41
// Elevation handler
#include "otbWrapperElevationParametersHandler.h"

42 43
namespace otb
{
44
namespace Wrapper
45 46
{

47
class ComputePolylineFeatureFromImage: public Application
48
{
49 50 51 52 53 54
public:
  /** Standard class typedefs. */
  typedef ComputePolylineFeatureFromImage Self;
  typedef Application Superclass;
  typedef itk::SmartPointer<Self> Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;
55
  // Images
56 57
  typedef FloatVectorImageType::PixelType PixelType;
  typedef FloatVectorImageType ImageType;
58 59

  // VectorData
60 61 62 63 64 65 66
  typedef VectorData<> VectorDataType;
  typedef VectorDataType::DataNodeType DataNodeType;
  typedef DataNodeType::ValuePrecisionType PrecisionType;
  typedef DataNodeType::PrecisionType CoordRepType;
  typedef itk::PreOrderTreeIterator<VectorDataType::DataTreeType> TreeIteratorType;
  typedef VectorDataIntoImageProjectionFilter<VectorDataType, ImageType> VectorDataIntoImageProjType;
  typedef VectorDataProjectionFilter<VectorDataType, VectorDataType> VectorDataProjectionFilterType;
67 68

  typedef ParserConditionDataNodeFeatureFunction<ImageType, CoordRepType, PrecisionType>
69 70 71 72
      ParserConditionFeatureFunctionType;

  /** Standard macro */
  itkNewMacro(Self)
OTB Bot's avatar
OTB Bot committed
73
;
74 75

  itkTypeMacro(ComputePolylineFeatureFromImage, otb::Application)
OTB Bot's avatar
OTB Bot committed
76
;
77 78

private:
79
  void DoInit() override
80 81
  {
    SetName("ComputePolylineFeatureFromImage");
82
    SetDescription("This application computes the chosen descriptors for each studied polyline contained in the input VectorData.");
83

84
    SetDocName("Compute Polyline Feature From Image");
85
    SetDocLongDescription("The first step in the classifier fusion based validation is to compute the chosen descriptors for each studied polyline.");
86
    SetDocLimitations("Since it does not rely on streaming process, take care of the size of input image before launching application.");
87 88 89 90 91
    SetDocAuthors("OTB-Team");
    SetDocSeeAlso(" ");
    AddDocTag(Tags::FeatureExtraction);

    AddParameter(ParameterType_InputImage, "in", "Input Image");
92
    SetParameterDescription("in", "An image to compute the descriptors on.");
93 94

    AddParameter(ParameterType_InputVectorData, "vd", "Vector Data");
95
    SetParameterDescription("vd", "Vector data containing the polylines where the features will be computed.");
96

97 98
    // Elevation
    ElevationParametersHandler::AddElevationParameters(this, "elev");
99 100

    AddParameter(ParameterType_String, "expr", "Feature expression");
101
    SetParameterDescription("expr", "The feature formula (b1 < 0.3) where b1 is the standard name of input image first band");
102 103

    AddParameter(ParameterType_String, "field", "Feature name");
104
    SetParameterDescription("field", "The field name corresponding to the feature codename (NONDVI, ROADSA...)");
105 106

    AddParameter(ParameterType_OutputVectorData, "out", "Output Vector Data");
107
    SetParameterDescription("out", "The output vector data containing polylines with a new field");
108

109 110 111 112 113 114
    // Doc example parameter settings
    SetDocExampleParameterValue("in", "NDVI.TIF");
    SetDocExampleParameterValue("vd", "roads_ground_truth.shp");
    SetDocExampleParameterValue("expr", "\"(b1 > 0.4)\"");
    SetDocExampleParameterValue("field", "NONDVI");
    SetDocExampleParameterValue("out", "PolylineFeatureFromImage_LI_NONDVI_gt.shp");
115

116
    SetOfficialDocLink();
117
  }
118

119
  void DoUpdateParameters() override
120 121 122 123
  {
    // Nothing to do here : all parameters are independent
  }

124
  void DoExecute() override
125 126 127 128
  {
    // Vector Data into Image projection
    FloatVectorImageType::Pointer inImage = GetParameterImage("in");

129 130 131 132 133 134 135
    inImage->UpdateOutputInformation();
    try
      {
      inImage->Update();
      }
    catch (...)
      {
OTB Bot's avatar
OTB Bot committed
136
      itkGenericExceptionMacro(<< "Error during image update, maybe your image is to big.");
137 138
      }

139 140 141 142 143 144
    otbAppLogDEBUG( << "Starting PolylineFeature extraction process" )

    // Vector Data into Image projection
    //// Read the Vectordata

    VectorDataType* inVectorData = GetParameterVectorData("vd");
145
    inVectorData->Update();
146 147 148 149 150
    //// Projection
    VectorDataIntoImageProjType::Pointer vprojIm = VectorDataIntoImageProjType::New();
    vprojIm->SetInputVectorData(inVectorData);
    vprojIm->SetInputImage(inImage);

151 152
    // Setup the DEM Handler
    otb::Wrapper::ElevationParametersHandler::SetupDEMHandlerFromElevationParameters(this,"elev");
153

154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190
    vprojIm->SetUseOutputSpacingAndOriginFromImage(true); // we want index as input;
    vprojIm->Update();

    // Add description
    ParserConditionFeatureFunctionType::Pointer vdescriptor = ParserConditionFeatureFunctionType::New();
    vdescriptor->SetExpression(GetParameterString("expr"));
    vdescriptor->SetInputImage(inImage);

    m_OutVectorData = VectorDataType::New();
    // Retrieving root node
    DataNodeType::Pointer root = m_OutVectorData->GetDataTree()->GetRoot()->Get();
    // Create the document node
    DataNodeType::Pointer document = DataNodeType::New();
    document->SetNodeType(otb::DOCUMENT);
    // Adding the layer to the data tree
    m_OutVectorData->GetDataTree()->Add(document, root);
    // Create the folder node
    DataNodeType::Pointer folder = DataNodeType::New();
    folder->SetNodeType(otb::FOLDER);
    // Adding the layer to the data tree
    m_OutVectorData->GetDataTree()->Add(folder, document);

    TreeIteratorType itVector(vprojIm->GetOutput()->GetDataTree());
    itVector.GoToBegin();

    while (!itVector.IsAtEnd())
      {
      vdescriptor->SetInputImage(inImage);
      if (!itVector.Get()->IsRoot() && !itVector.Get()->IsDocument() && !itVector.Get()->IsFolder())
        {
        DataNodeType::Pointer currentGeometry = itVector.Get();
        currentGeometry->SetFieldAsDouble(GetParameterString("field"),
                                          (double) (vdescriptor->Evaluate(*(currentGeometry.GetPointer()))[0]));
        m_OutVectorData->GetDataTree()->Add(currentGeometry, folder);
        }
      ++itVector;
      }
191

192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
    /*
     * Reprojection of the output VectorData
     *
     * The output of VectorDataToRoadDescription is in image index coordinates
     *
     * 3 cases :
     * - input image has no geo-information : pass through
     * - input image is in cartographic projection : apply image spacing and origin, and set the ProjectionRef
     * - input image is in sensor model geometry : reproject in WGS84
     *
     */

    std::string projRef = inImage->GetProjectionRef();
    ImageKeywordlist kwl;
    itk::ExposeMetaData<ImageKeywordlist>(inImage->GetMetaDataDictionary(), MetaDataKey::OSSIMKeywordlistKey, kwl);

    if (!projRef.empty())
      {
      // image is in cartographic projection
      // apply spacing and origin + set projection WKT
      // The VectorData in output of the chain is in image index coordinate,
      // and the projection information is lost
      // Apply an affine transform to apply image origin and spacing,
      // and arbitrarily set the ProjectionRef to the input image ProjectionRef

      typedef itk::AffineTransform<VectorDataType::PrecisionType, 2> TransformType;
      typedef otb::VectorDataTransformFilter<VectorDataType, VectorDataType> VDTransformType;

      TransformType::ParametersType params;
      params.SetSize(6);
222
      params[0] = inImage->GetSignedSpacing()[0];
223 224
      params[1] = 0;
      params[2] = 0;
225
      params[3] = inImage->GetSignedSpacing()[1];
226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252
      params[4] = inImage->GetOrigin()[0];
      params[5] = inImage->GetOrigin()[1];

      TransformType::Pointer transform = TransformType::New();
      transform->SetParameters(params);

      VDTransformType::Pointer vdTransform = VDTransformType::New();
      vdTransform->SetTransform(transform);
      vdTransform->SetInput(m_OutVectorData);
      vdTransform->Update();

      m_ProjectedVectorData = vdTransform->GetOutput();

      m_ProjectedVectorData->SetProjectionRef(inImage->GetProjectionRef());
      }
    else
      if (kwl.GetSize() > 0)
        {
        // image is in sensor model geometry

        // Reproject VectorData in image projection
        VectorDataProjectionFilterType::Pointer vproj = VectorDataProjectionFilterType::New();
        vproj->SetInput(m_OutVectorData);

        vproj->SetOutputKeywordList(inImage->GetImageKeywordlist());
        vproj->SetOutputProjectionRef(inImage->GetProjectionRef());
        vproj->SetOutputOrigin(inImage->GetOrigin());
253
        vproj->SetOutputSpacing(inImage->GetSignedSpacing());
254

255 256
        // Setup the DEM Handler
        otb::Wrapper::ElevationParametersHandler::SetupDEMHandlerFromElevationParameters(this,"elev");
257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275

        vproj->Update();
        m_ProjectedVectorData = vproj->GetOutput();
        }
      else
        {
        // no georeferencing information
        m_ProjectedVectorData = m_OutVectorData;
        }

    // Set the output vectorData
    SetParameterOutputVectorData("out", m_ProjectedVectorData);

  };

  VectorDataType::Pointer m_ProjectedVectorData;
  VectorDataType::Pointer m_OutVectorData;

};
276

277 278
}
}
279

280
OTB_APPLICATION_EXPORT(otb::Wrapper::ComputePolylineFeatureFromImage)
281 282