An error occurred while loading the file. Please try again.
-
Victor Poughon authored6b3198ce
/*
* 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.
*/
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
#include "otbOGRDataSourceWrapper.h"
#include "otbOGRFeatureWrapper.h"
#include "itkVariableLengthVector.h"
#include "otbStatisticsXMLFileReader.h"
#include "itkListSample.h"
#include "otbShiftScaleSampleListFilter.h"
#include "otbDimensionalityReductionModelFactory.h"
#include <time.h>
namespace otb
{
namespace Wrapper
{
/** Utility function to negate std::isalnum */
bool IsNotAlphaNum(char c)
{
return !std::isalnum(c);
}
/**
* \class VectorDimensionalityReduction
*
* Apply a dimensionality reduction model on a vector file
*/
class VectorDimensionalityReduction : public Application
{
public:
/** Standard class typedefs. */
typedef VectorDimensionalityReduction Self;
typedef Application Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Standard macro */
itkNewMacro(Self);
itkTypeMacro(Self, Application)
/** Filters typedef */
typedef float ValueType;
typedef itk::VariableLengthVector<ValueType> InputSampleType;
typedef itk::Statistics::ListSample<InputSampleType> ListSampleType;
typedef MachineLearningModel<
itk::VariableLengthVector<ValueType>,
itk::VariableLengthVector<ValueType> > DimensionalityReductionModelType;
typedef DimensionalityReductionModelFactory<
ValueType,ValueType> DimensionalityReductionModelFactoryType;
typedef DimensionalityReductionModelType::Pointer ModelPointerType;
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
/** Statistics Filters typedef */
typedef itk::VariableLengthVector<ValueType> MeasurementType;
typedef otb::StatisticsXMLFileReader<MeasurementType> StatisticsReader;
typedef otb::Statistics::ShiftScaleSampleListFilter<
ListSampleType, ListSampleType> ShiftScaleFilterType;
protected:
~VectorDimensionalityReduction() override
{
DimensionalityReductionModelFactoryType::CleanFactories();
}
private:
void DoInit() override
{
SetName("VectorDimensionalityReduction");
SetDescription("Performs dimensionality reduction of the input vector data "
"according to a model file.");
SetDocName("Vector Dimensionality Reduction");
SetDocAuthors("OTB-Team");
SetDocLongDescription("This application performs a vector data "
"dimensionality reduction based on a model file produced by the "
"TrainDimensionalityReduction application.");
SetDocSeeAlso("TrainDimensionalityReduction");
SetDocLimitations("None");
AddDocTag(Tags::Learning);
AddParameter(ParameterType_InputVectorData, "in", "Name of the input vector data");
SetParameterDescription("in","The input vector data to reduce.");
AddParameter(ParameterType_InputFilename, "instat", "Statistics file");
SetParameterDescription("instat", "A XML file containing mean and standard "
"deviation to center and reduce samples before dimensionality reduction "
"(produced by ComputeImagesStatistics application).");
MandatoryOff("instat");
AddParameter(ParameterType_InputFilename, "model", "Model file");
SetParameterDescription("model", "A model file (produced by the "
"TrainDimensionalityReduction application,");
AddParameter(ParameterType_OutputFilename, "out", "Output vector data file "
"containing the reduced vector");
SetParameterDescription("out","Output vector data file storing sample "
"values (OGR format). If not given, the input vector data file is used. "
"In overwrite mode, the original features will be lost.");
MandatoryOff("out");
AddParameter(ParameterType_ListView, "feat", "Input features to use for reduction."); //
SetParameterDescription("feat","List of field names in the input vector "
"data used as features for reduction."); //
AddParameter(ParameterType_Choice, "featout", "Output feature"); //
SetParameterDescription("featout", "Naming of output features");
AddChoice("featout.prefix", "Prefix");
SetParameterDescription("featout.prefix", "Use a name prefix");
AddParameter(ParameterType_String, "featout.prefix.name", "Feature name prefix");
SetParameterDescription("featout.prefix.name","Name prefix for output "
"features. This prefix is followed by the numeric index of each output feature.");
SetParameterString("featout.prefix.name","reduced_", false);
AddChoice("featout.list","List");
SetParameterDescription("featout.list", "Use a list with all names");
AddParameter(ParameterType_StringList, "featout.list.names", "Feature name list");
SetParameterDescription("featout.list.names","List of field names for the output "
"features which result from the reduction."); //
AddParameter(ParameterType_Int, "pcadim", "Principal component dimension"); //
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
SetParameterDescription("pcadim","This optional parameter can be set to "
"reduce the number of eignevectors used in the PCA model file. This "
"parameter can't be used for other models"); //
MandatoryOff("pcadim");
AddParameter(ParameterType_Choice, "mode", "Writing mode"); //
SetParameterDescription("mode", "This parameter determines if the output "
"file is overwritten or updated [overwrite/update]. If an output file "
"name is given, the original file is copied before creating the new features.");
AddChoice("mode.overwrite", "Overwrite");
SetParameterDescription("mode.overwrite","Overwrite mode"); //
AddChoice("mode.update", "Update");
SetParameterDescription("mode.update", "Update mode");
// Doc example parameter settings
SetDocExampleParameterValue("in", "vectorData.shp");
SetDocExampleParameterValue("instat", "meanVar.xml");
SetDocExampleParameterValue("model", "model.txt");
SetDocExampleParameterValue("out", "vectorDataOut.shp");
SetDocExampleParameterValue("feat", "perimeter area width");
//SetOfficialDocLink();
}
void DoUpdateParameters() override
{
if ( HasValue("in") )
{
std::string shapefile = GetParameterString("in");
otb::ogr::DataSource::Pointer ogrDS;
OGRSpatialReference oSRS("");
std::vector<std::string> options;
ogrDS = otb::ogr::DataSource::New(shapefile, otb::ogr::DataSource::Modes::Read);
otb::ogr::Layer layer = ogrDS->GetLayer(0);
OGRFeatureDefn &layerDefn = layer.GetLayerDefn();
ClearChoices("feat");
for(int iField=0; iField< layerDefn.GetFieldCount(); iField++)
{
std::string item = layerDefn.GetFieldDefn(iField)->GetNameRef();
std::string key(item);
std::string::iterator end = std::remove_if( key.begin(), key.end(), IsNotAlphaNum );
std::transform( key.begin(), end, key.begin(), tolower );
std::string tmpKey = "feat." + key.substr( 0, static_cast<unsigned long>( end - key.begin() ) );
AddChoice(tmpKey,item);
}
}
}
void DoExecute() override
{
clock_t tic = clock();
std::string shapefile = GetParameterString("in");
otb::ogr::DataSource::Pointer source = otb::ogr::DataSource::New(
shapefile, otb::ogr::DataSource::Modes::Read);
otb::ogr::Layer layer = source->GetLayer(0);
ListSampleType::Pointer input = ListSampleType::New();
std::vector<int> inputIndexes = GetSelectedItems("feat");
int nbFeatures = inputIndexes.size();
input->SetMeasurementVectorSize(nbFeatures);
otb::ogr::Layer::const_iterator it = layer.cbegin();
otb::ogr::Layer::const_iterator itEnd = layer.cend();
// Get the list of non-selected field indexes
// /!\ The 'feat' is assumed to expose all available fields, hence the
// mapping between GetSelectedItems() and OGR field indexes
OGRFeatureDefn &inLayerDefn = layer.GetLayerDefn();
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
std::set<int> otherInputFields;
for (int i=0 ; i < inLayerDefn.GetFieldCount() ; i++)
otherInputFields.insert(i);
for (int k=0 ; k < nbFeatures ; k++)
otherInputFields.erase(inputIndexes[k]);
for( ; it!=itEnd ; ++it)
{
MeasurementType mv;
mv.SetSize(nbFeatures);
for(int idx=0; idx < nbFeatures; ++idx)
{
mv[idx] = static_cast<float>( (*it)[inputIndexes[idx]].GetValue<double>() );
}
input->PushBack(mv);
}
/** Statistics for shift/scale */
MeasurementType meanMeasurementVector;
MeasurementType stddevMeasurementVector;
if (HasValue("instat") && IsParameterEnabled("instat"))
{
StatisticsReader::Pointer statisticsReader = StatisticsReader::New();
std::string XMLfile = GetParameterString("instat");
statisticsReader->SetFileName(XMLfile);
meanMeasurementVector = statisticsReader->GetStatisticVectorByName("mean");
stddevMeasurementVector = statisticsReader->GetStatisticVectorByName("stddev");
otbAppLogINFO("Mean used: " << meanMeasurementVector);
otbAppLogINFO("Standard deviation used: " << stddevMeasurementVector);
}
else
{
meanMeasurementVector.SetSize(nbFeatures);
meanMeasurementVector.Fill(0.);
stddevMeasurementVector.SetSize(nbFeatures);
stddevMeasurementVector.Fill(1.);
}
ShiftScaleFilterType::Pointer trainingShiftScaleFilter = ShiftScaleFilterType::New();
trainingShiftScaleFilter->SetInput(input);
trainingShiftScaleFilter->SetShifts(meanMeasurementVector);
trainingShiftScaleFilter->SetScales(stddevMeasurementVector);
trainingShiftScaleFilter->Update();
otbAppLogINFO("Loading model");
/** Read the model */
m_Model = DimensionalityReductionModelFactoryType::CreateDimensionalityReductionModel(
GetParameterString("model"),
DimensionalityReductionModelFactoryType::ReadMode);
if (m_Model.IsNull())
{
otbAppLogFATAL(<< "Error when loading model " << GetParameterString("model")
<< " : unsupported model type");
}
m_Model->Load(GetParameterString("model"));
if (HasValue("pcadim") && IsParameterEnabled("pcadim"))
{
std::string modelName(m_Model->GetNameOfClass());
if (modelName != "PCAModel")
{
otbAppLogFATAL(<< "Can't set 'pcadim' on a model : "<< modelName);
}
m_Model->SetDimension( GetParameterInt("pcadim") );
}
otbAppLogINFO("Model loaded, dimension : "<< m_Model->GetDimension());
/** Perform Dimensionality Reduction */
ListSampleType::Pointer listSample = trainingShiftScaleFilter->GetOutput();
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
ListSampleType::Pointer target = m_Model->PredictBatch(listSample);
/** Create/Update Output Shape file */
ogr::DataSource::Pointer output;
ogr::DataSource::Pointer buffer = ogr::DataSource::New();
bool updateMode = false;
if (IsParameterEnabled("out") && HasValue("out"))
{
// Create new OGRDataSource
if (GetParameterString("mode")=="overwrite")
{
output = ogr::DataSource::New(GetParameterString("out"), ogr::DataSource::Modes::Overwrite);
otb::ogr::Layer newLayer = output->CreateLayer(
GetParameterString("out"),
const_cast<OGRSpatialReference*>(layer.GetSpatialRef()),
layer.GetGeomType());
// Copy existing fields except the ones selected for reduction
for (const int& k : otherInputFields)
{
OGRFieldDefn fieldDefn(inLayerDefn.GetFieldDefn(k));
newLayer.CreateField(fieldDefn);
}
}
else if (GetParameterString("mode")=="update")
{
//output = ogr::DataSource::New(GetParameterString("out"), ogr::DataSource::Modes::Update_LayerCreateOnly);
// Update mode
otb::ogr::DataSource::Pointer source_output =
otb::ogr::DataSource::New(GetParameterString("out"), otb::ogr::DataSource::Modes::Read);
layer = source_output->GetLayer(0);
updateMode = true;
otbAppLogINFO("Update input vector data.");
// fill temporary buffer for the transfer
otb::ogr::Layer inputLayer = layer;
layer = buffer->CopyLayer(inputLayer, std::string("Buffer"));
// close input data source
source_output->Clear();
// Re-open input data source in update mode
output = otb::ogr::DataSource::New(
GetParameterString("out"),
otb::ogr::DataSource::Modes::Update_LayerUpdate);
}
else
{
otbAppLogFATAL(<< "Error when creating the output file" <<
GetParameterString("mode") << " : unsupported writing mode type");
}
}
otb::ogr::Layer outLayer = output->GetLayer(0);
OGRErr errStart = outLayer.ogr().StartTransaction();
if (errStart != OGRERR_NONE)
{
otbAppLogFATAL(<< "Unable to start transaction for OGR layer " << outLayer.ogr().GetName() << ".");
}
// Build the list of output fields
std::vector<std::string> outFields;
if(GetParameterString("featout") == "prefix")
{
std::string prefix = GetParameterString("featout.prefix.name");
std::ostringstream oss;
for (unsigned int i=0 ; i < m_Model->GetDimension() ; i++)
{
oss.str(prefix);
oss.seekp(0,std::ios_base::end);
oss << i;
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
outFields.push_back(oss.str());
}
}
else if(GetParameterString("featout") == "list")
{
outFields = GetParameterStringList("featout.list.names");
if (outFields.size() != m_Model->GetDimension())
{
otbAppLogFATAL( << "Wrong number of output field names, expected "
<< m_Model->GetDimension() << " , got "<< outFields.size());
}
}
else
{
otbAppLogFATAL( << "Unsupported output feature mode : "
<< GetParameterString("featout"));
}
// Add the field of prediction in the output layer if field not exist
for (unsigned int i=0; i<outFields.size() ;i++)
{
OGRFeatureDefn &layerDefn = outLayer.GetLayerDefn();
int idx = layerDefn.GetFieldIndex(outFields[i].c_str());
if (idx >= 0)
{
if (layerDefn.GetFieldDefn(idx)->GetType() != OFTReal)
otbAppLogFATAL("Field name "<< outFields[i]
<< " already exists with a different type!");
}
else
{
OGRFieldDefn predictedField(outFields[i].c_str(), OFTReal);
ogr::FieldDefn predictedFieldDef(predictedField);
outLayer.CreateField(predictedFieldDef);
}
}
// Fill output layer
unsigned int count=0;
it = layer.cbegin();
itEnd = layer.cend();
for( ; it!=itEnd ; ++it, ++count)
{
ogr::Feature dstFeature(outLayer.GetLayerDefn());
dstFeature.SetFrom( *it , TRUE);
dstFeature.SetFID(it->GetFID());
for (std::size_t i=0; i<outFields.size(); ++i)
{
dstFeature[outFields[i]].SetValue<double>(target->GetMeasurementVector(count)[i]);
}
if (updateMode)
{
outLayer.SetFeature(dstFeature);
}
else
{
outLayer.CreateFeature(dstFeature);
}
}
if(outLayer.ogr().TestCapability("Transactions"))
{
const OGRErr errCommitX = outLayer.ogr().CommitTransaction();
if (errCommitX != OGRERR_NONE)
{
otbAppLogFATAL(<< "Unable to commit transaction for OGR layer " <<
outLayer.ogr().GetName() << ".");
421422423424425426427428429430431432433434435
}
}
output->SyncToDisk();
clock_t toc = clock();
otbAppLogINFO( "Elapsed: "<< ((double)(toc - tic) / CLOCKS_PER_SEC)<<" seconds.");
}
ModelPointerType m_Model;
};
} // end of namespace Wrapper
} // end of namespace otb
OTB_APPLICATION_EXPORT(otb::Wrapper::VectorDimensionalityReduction)