diff --git a/Applications/Segmentation/otbSegmentation.cxx b/Applications/Segmentation/otbSegmentation.cxx
index 1eb6a8a6a04258dec9a49e150043258877c93bbc..325713a231a1dfdf863ae6730e7805ffadae2109 100644
--- a/Applications/Segmentation/otbSegmentation.cxx
+++ b/Applications/Segmentation/otbSegmentation.cxx
@@ -29,6 +29,13 @@
#include "otbVectorImageToAmplitudeImageFilter.h"
#include "itkGradientMagnitudeImageFilter.h"
#include "otbWatershedSegmentationFilter.h"
+#include "otbMultiScaleConvexOrConcaveClassificationFilter.h"
+#include "itkBinaryBallStructuringElement.h"
+#include "otbMorphologicalOpeningProfileFilter.h"
+#include "otbMorphologicalClosingProfileFilter.h"
+#include "otbProfileToProfileDerivativeFilter.h"
+#include "otbProfileDerivativeToMultiScaleCharacteristicsFilter.h"
+#include "itkScalarConnectedComponentImageFilter.h"
// Large scale vectorization framework
#include "otbStreamingImageToOGRLayerSegmentationFilter.h"
@@ -38,7 +45,6 @@
#include "otbOGRLayerStreamStitchingFilter.h"
#include "otbGeoInformationConversion.h"
-
namespace otb
{
namespace Wrapper
@@ -79,6 +85,11 @@ public:
FunctorType,
MaskImageType > ConnectedComponentSegmentationFilterType;
+ typedef itk::ScalarConnectedComponentImageFilter
+ <LabelImageType,
+ LabelImageType> LabeledConnectedComponentSegmentationFilterType;
+
+
// Watershed
typedef otb::VectorImageToAmplitudeImageFilter
<FloatVectorImageType,
@@ -90,6 +101,18 @@ public:
typedef otb::WatershedSegmentationFilter
<FloatImageType,LabelImageType> WatershedSegmentationFilterType;
+ // Geodesic morphology multiscale segmentation
+ typedef itk::BinaryBallStructuringElement<FloatImageType::PixelType, 2> StructuringElementType;
+ typedef otb::MorphologicalOpeningProfileFilter<FloatImageType, FloatImageType, StructuringElementType>
+ OpeningProfileFilterType;
+ typedef otb::MorphologicalClosingProfileFilter<FloatImageType, FloatImageType, StructuringElementType>
+ ClosingProfileFilterType;
+ typedef otb::ProfileToProfileDerivativeFilter<FloatImageType, FloatImageType> DerivativeFilterType;
+ typedef otb::ProfileDerivativeToMultiScaleCharacteristicsFilter<FloatImageType,FloatImageType, LabelImageType>
+ MultiScaleCharacteristicsFilterType;
+ typedef otb::MultiScaleConvexOrConcaveClassificationFilter<FloatImageType,LabelImageType> MultiScaleClassificationFilterType;
+
+
// mask filter
typedef otb::MaskMuParserFilter
@@ -112,6 +135,12 @@ public:
<FloatVectorImageType,
ConnectedComponentSegmentationFilterType>
ConnectedComponentStreamingVectorizedSegmentationOGRType;
+ typedef otb::StreamingImageToOGRLayerSegmentationFilter
+ <LabelImageType,
+ LabeledConnectedComponentSegmentationFilterType>
+ LabeledConnectedComponentStreamingVectorizedSegmentationOGRType;
+
+
typedef otb::OGRLayerStreamStitchingFilter
<FloatVectorImageType> FusionFilterType;
@@ -134,7 +163,7 @@ private:
// Documentation
SetDocName("Segmentation");
- SetDocLongDescription("This application allows to perform various segmentation algorithms on an multispectral image."
+ SetDocLongDescription("This application allows to perform various segmentation algorithms on a multispectral image."
"Available segmentation algorithms are two different versions of Mean-Shift segmentation algorithm (one being multi-threaded),"
" simple pixel based connected components according to a user-defined criterion, and watershed from the gradient of the intensity"
" (norm of spectral bands vector). The application has two different modes that affects the nature of its output.\n\nIn raster mode,"
@@ -241,6 +270,22 @@ private:
AddChoice("mode.raster", "Standard segmentation with labeled raster output");
SetParameterDescription("mode.raster","In this mode, the application will output a standard labeled raster. This mode can not handle large data.");
+ // GeoMorpho
+ AddChoice("filter.geomorpho","Multiscale Geodesic Morphology");
+
+ AddParameter(ParameterType_Int,"filter.geomorpho.size","Profile Size");
+ SetParameterDescription("filter.geomorpho.size","TODO");
+ SetDefaultParameterInt("filter.geomorpho.size",5);
+ AddParameter(ParameterType_Int,"filter.geomorpho.step","Profile Step");
+ SetParameterDescription("filter.geomorpho.step","TODO");
+ SetDefaultParameterInt("filter.geomorpho.step",1);
+ AddParameter(ParameterType_Int,"filter.geomorpho.start","Profile Start");
+ SetParameterDescription("filter.geomorpho.start","TODO");
+ SetDefaultParameterInt("filter.geomorpho.start",1);
+ AddParameter(ParameterType_Float,"filter.geomorpho.sigma","Sigma for decision rule");
+ SetParameterDescription("filter.geomorpho.sigma","TODO");
+ SetDefaultParameterFloat("filter.geomorpho.sigma",1.);
+
//Raster mode parameters
AddParameter(ParameterType_OutputImage, "mode.raster.out", "Output labeled image");
SetParameterDescription( "mode.raster.out", "The output labeled image.");
@@ -345,7 +390,7 @@ private:
typedef TSegmentationFilter SegmentationFilterType;
typedef typename SegmentationFilterType::Pointer SegmentationFilterPointerType;
typedef otb::StreamingImageToOGRLayerSegmentationFilter
- <FloatVectorImageType,
+ <TInputImage,
SegmentationFilterType> StreamingVectorizedSegmentationOGRType;
// Retrieve tile size parameter
@@ -637,6 +682,73 @@ private:
gradientMagnitudeFilter->GetOutput(),
layer, 0);
}
+
+ else
+ if (segType == "geomorpho")
+ {
+ otbAppLogINFO(<<"Using multiscale geodesic morphology segmentation."<<std::endl);
+
+
+ unsigned int profileSize = GetParameterInt("filter.geomorpho.size");
+ unsigned int initialValue = GetParameterInt("filter.geomorpho.start");
+ unsigned int step = GetParameterInt("filter.geomorpho.step");
+ double sigma = GetParameterFloat("filter.geomorpho.sigma");
+
+
+ AmplitudeFilterType::Pointer amplitudeFilter = AmplitudeFilterType::New();
+
+ amplitudeFilter->SetInput(this->GetParameterFloatVectorImage("in"));
+
+ OpeningProfileFilterType::Pointer oprofileFilter = OpeningProfileFilterType::New();
+ oprofileFilter->SetInput(amplitudeFilter->GetOutput());
+ oprofileFilter->SetProfileSize(profileSize);
+ oprofileFilter->SetInitialValue(initialValue);
+ oprofileFilter->SetStep(step);
+
+ ClosingProfileFilterType::Pointer cprofileFilter = ClosingProfileFilterType::New();
+ cprofileFilter->SetInput(amplitudeFilter->GetOutput());
+ cprofileFilter->SetProfileSize(profileSize);
+ cprofileFilter->SetInitialValue(initialValue);
+ cprofileFilter->SetStep(step);
+
+ DerivativeFilterType::Pointer oderivativeFilter = DerivativeFilterType::New();
+ oderivativeFilter->SetInput(oprofileFilter->GetOutput());
+
+ DerivativeFilterType::Pointer cderivativeFilter = DerivativeFilterType::New();
+ cderivativeFilter->SetInput(cprofileFilter->GetOutput());
+
+ MultiScaleCharacteristicsFilterType::Pointer omsCharFilter = MultiScaleCharacteristicsFilterType::New();
+ omsCharFilter->SetInput(oderivativeFilter->GetOutput());
+ omsCharFilter->SetInitialValue(initialValue);
+ omsCharFilter->SetStep(step);
+
+ MultiScaleCharacteristicsFilterType::Pointer cmsCharFilter = MultiScaleCharacteristicsFilterType::New();
+ cmsCharFilter->SetInput(cderivativeFilter->GetOutput());
+ cmsCharFilter->SetInitialValue(initialValue);
+ cmsCharFilter->SetStep(step);
+
+ MultiScaleClassificationFilterType::Pointer classificationFilter = MultiScaleClassificationFilterType::New();
+ classificationFilter->SetOpeningProfileDerivativeMaxima(omsCharFilter->GetOutput());
+ classificationFilter->SetOpeningProfileCharacteristics(omsCharFilter->GetOutputCharacteristics());
+ classificationFilter->SetClosingProfileDerivativeMaxima(cmsCharFilter->GetOutput());
+ classificationFilter->SetClosingProfileCharacteristics(cmsCharFilter->GetOutputCharacteristics());
+ classificationFilter->SetSigma(sigma);
+ classificationFilter->SetLabelSeparator(initialValue + profileSize * step);
+
+ LabeledConnectedComponentStreamingVectorizedSegmentationOGRType::Pointer
+ ccVectorizationFilter = LabeledConnectedComponentStreamingVectorizedSegmentationOGRType::New();
+
+ if (HasValue("mode.vector.inmask"))
+ {
+ ccVectorizationFilter->GetSegmentationFilter()->SetMaskImage(
+ this->GetParameterUInt32Image("mode.vector.inmask"));
+ }
+
+ streamSize = GenericApplySegmentation<LabelImageType, LabeledConnectedComponentSegmentationFilterType> (
+ ccVectorizationFilter, classificationFilter->GetOutput(),
+ layer, 0);
+
+ }
else
{
otbAppLogFATAL(<<"non defined filtering method "<<GetParameterInt("filter")<<std::endl);
diff --git a/Code/UtilitiesAdapters/OssimAdapters/otbSensorModelAdapter.cxx b/Code/UtilitiesAdapters/OssimAdapters/otbSensorModelAdapter.cxx
index a8df4fbbb701be967a7c7f35f23cd209500b11bc..c7ef65c8fe4e731accbd5a833889d0e71de2f97f 100644
--- a/Code/UtilitiesAdapters/OssimAdapters/otbSensorModelAdapter.cxx
+++ b/Code/UtilitiesAdapters/OssimAdapters/otbSensorModelAdapter.cxx
@@ -120,7 +120,9 @@ void SensorModelAdapter::ForwardTransformPoint(double x, double y, double z,
otbMsgDevMacro(<< "USING DEM ! ");
- while ((diffHeight > m_Epsilon) && (nbIter < m_NbIter))
+ bool nanHeight = false;
+
+ while (!nanHeight && (diffHeight > m_Epsilon) && (nbIter < m_NbIter))
{
otbMsgDevMacro(<< "Iter " << nbIter);
@@ -128,10 +130,18 @@ void SensorModelAdapter::ForwardTransformPoint(double x, double y, double z,
heightTmp = this->m_DEMHandler->GetHeightAboveMSL(lon, lat);
- this->m_SensorModel->lineSampleHeightToWorld(ossimPoint, heightTmp, ossimGPointRef);
-
+ if(ossim::isnan(heightTmp))
+ {
+ nanHeight = true;
+ this->m_SensorModel->lineSampleToWorld(ossimPoint, ossimGPointRef);
+ }
+ else
+ {
+ this->m_SensorModel->lineSampleHeightToWorld(ossimPoint, heightTmp, ossimGPointRef);
+ }
+
diffHeight = fabs(heightTmp - height);
-
+
++nbIter;
}
ossimGPoint = ossimGPointRef;
@@ -245,4 +255,41 @@ double SensorModelAdapter::Optimize()
return precision;
}
+bool SensorModelAdapter::ReadGeomFile(const std::string & infile)
+{
+ ossimKeywordlist geom;
+
+ geom.add(infile.c_str());
+
+ m_SensorModel = ossimSensorModelFactory::instance()->createProjection(geom);
+
+ if (m_SensorModel == NULL)
+ {
+ m_SensorModel = ossimplugins::ossimPluginProjectionFactory::instance()->createProjection(geom);
+ }
+
+ return (m_SensorModel != NULL);
+}
+
+bool SensorModelAdapter::WriteGeomFile(const std::string & outfile)
+{
+ // If tie points and model are allocated
+ if(m_SensorModel != NULL)
+ {
+ // try to retrieve a sensor model
+ ossimSensorModel * sensorModel = NULL;
+ sensorModel = dynamic_cast<ossimSensorModel *>(m_SensorModel);
+
+ ossimKeywordlist geom;
+
+ bool success = sensorModel->saveState(geom);
+
+ if(success)
+ {
+ return geom.write(outfile.c_str());
+ }
+ }
+ return false;
+}
+
} // namespace otb
diff --git a/Code/UtilitiesAdapters/OssimAdapters/otbSensorModelAdapter.h b/Code/UtilitiesAdapters/OssimAdapters/otbSensorModelAdapter.h
index c5f67ee978161e16bd0bf948eb7c97bd41cf23fd..150feab485f9bf1fdbb2085fb01d393de9b8fe03 100644
--- a/Code/UtilitiesAdapters/OssimAdapters/otbSensorModelAdapter.h
+++ b/Code/UtilitiesAdapters/OssimAdapters/otbSensorModelAdapter.h
@@ -86,6 +86,12 @@ public:
/** Is sensor model valid method. return false if the m_SensorModel is null*/
bool IsValidSensorModel();
+ /** Read geom file and instanciate sensor model */
+ bool ReadGeomFile(const std::string & infile);
+
+ /** Write geom file corresponding to sensor model */
+ bool WriteGeomFile(const std::string& outfile);
+
protected:
SensorModelAdapter();
virtual ~SensorModelAdapter();