diff --git a/.hgsigs b/.hgsigs new file mode 100644 index 0000000000000000000000000000000000000000..752dbc56655223c129380da99fb0d17cccd046b3 --- /dev/null +++ b/.hgsigs @@ -0,0 +1 @@ +68eebc1b170a636794feda66523e6c0135ebfef2 0 iEYEABECAAYFAkmvX3cACgkQwRJnCg+r8KFSCgCfVYSnx2ev+hIlpbM/arxzLVf3KvsAn3nqaJXaj62RcNUIjv+qcNWlJ5WG diff --git a/Code/BasicFilters/otbAmplitudePhaseToRGBFunctor.h b/Code/BasicFilters/otbAmplitudePhaseToRGBFunctor.h index 41f51633a49879aeb1e7602ec125a669ada0e719..e74c42dd5bce5785ec4750441a17242756142809 100644 --- a/Code/BasicFilters/otbAmplitudePhaseToRGBFunctor.h +++ b/Code/BasicFilters/otbAmplitudePhaseToRGBFunctor.h @@ -42,7 +42,7 @@ namespace otb * and itk::ComplexToPhaseImageFilter for example. * */ - template< class TInput1, class TInput2=TInput1, class TOutput=TInput1> + template< class TInput1, class TInput2=TInput1, class TInput3=TInput1, class TOutput=TInput1> class ITK_EXPORT AmplitudePhaseToRGBFunctor { public: @@ -68,7 +68,7 @@ namespace otb this->m_Minimum = min; } - inline TOutput operator()( const TInput1 & amplitude, const TInput2 & phase) const + inline TOutput operator()( const TInput1 & amplitude, const TInput2 & coherence, const TInput3 & phase) const { // std::cout << amplitude << " - " << phase << std::endl; double hinc, sinc, vinc; @@ -79,8 +79,9 @@ namespace otb double hue, sat, val; hue = 0.6 - (phase+M_PI)*hinc; - sat = 0.99; - val = itk::NumericTraits<RGBComponentType>::max()/(m_Maximum-m_Minimum) * (amplitude-m_Minimum); + sat = 0.6*coherence+0.3; + val = itk::NumericTraits<RGBComponentType>::max()/2 + * ( (amplitude-m_Minimum)/(m_Maximum-m_Minimum)+1.0); if (amplitude < m_Minimum) { diff --git a/Code/BasicFilters/otbStreamingResampleImageFilter.txx b/Code/BasicFilters/otbStreamingResampleImageFilter.txx index ad1ea471e1fb7baffe9ab57c6870ef12c9a7dbc8..500d998505c40a53b9d0cdea7c278930ea929eb0 100644 --- a/Code/BasicFilters/otbStreamingResampleImageFilter.txx +++ b/Code/BasicFilters/otbStreamingResampleImageFilter.txx @@ -88,7 +88,8 @@ StreamingResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType vPoints.push_back(indexTmp); //otbGenericMsgDebugMacro(<< "indexLL : (" << indexTmp[0] << "," << indexTmp[1] << ")"); - typedef itk::ContinuousIndex<TInterpolatorPrecisionType, 2> ContinuousIndexType; + typedef typename Superclass::CoordRepType CoordRepType; + typedef itk::ContinuousIndex<CoordRepType, 2> ContinuousIndexType; typename ContinuousIndexType::ValueType minX = itk::NumericTraits<typename ContinuousIndexType::ValueType>::max(); typename ContinuousIndexType::ValueType maxX = 0; typename ContinuousIndexType::ValueType minY = itk::NumericTraits<typename ContinuousIndexType::ValueType>::max(); diff --git a/Code/ChangeDetection/otbCBAMIChangeDetector.h b/Code/ChangeDetection/otbCBAMIChangeDetector.h index afe95fbae66bf5412e83c2c24ae56da85813e89d..b5e012c05643a64ac10153c701a67002d815cf2d 100644 --- a/Code/ChangeDetection/otbCBAMIChangeDetector.h +++ b/Code/ChangeDetection/otbCBAMIChangeDetector.h @@ -73,7 +73,7 @@ public: inline TOutput operator()( const TInput1 & itA, const TInput2 & itB) { - + const double epsilon = 0.01; VectorType vecA; VectorType vecB; @@ -93,8 +93,6 @@ public: protected: - static const double epsilon = 0.01; - inline void normalizeInPlace(VectorType vx) { diff --git a/Code/Common/otbVectorDataExtractROI.txx b/Code/Common/otbVectorDataExtractROI.txx index d7ec133c361818dbd5dc7ed4ab67328f98bd844a..703ba6f978b22a555ab833fa5b7e20ce031d2bb8 100644 --- a/Code/Common/otbVectorDataExtractROI.txx +++ b/Code/Common/otbVectorDataExtractROI.txx @@ -64,10 +64,11 @@ VectorDataExtractROI<TVectorData> typename VectorDataType::Pointer output = this->GetOutput(); /** put this here*/ - output->SetProjectionRef(input->GetProjectionRef()); + if(!input->GetProjectionRef().empty()) + output->SetProjectionRef(input->GetProjectionRef()); if(!input) - std::cout << " Probleme avec la recuperation du input"<<std::endl; + return; /** Need to check if it is necessary to project the roi*/ this->CompareInputAndRegionProjection(); @@ -76,6 +77,8 @@ VectorDataExtractROI<TVectorData> /** If Projection of the region is needed, we project on the vectorData coordinate axis*/ if(m_ProjectionNeeded) this->ProjectRegionToInputVectorProjection(); + else + m_GeoROI = m_ROI; /** Loop in the vectorData file @@ -87,49 +90,95 @@ VectorDataExtractROI<TVectorData> typename DataNodeType::Pointer root = input->GetDataTree()->GetRoot()->Get(); typename VectorDataType::DataTreeType::Pointer tree = output->GetDataTree(); - DataNodePointerType currentContainer; + DataNodePointerType currentContainer; + DataNodePointerType newDataNodeFolder = DataNodeType::New(); + DataNodePointerType newDataNodeMultiPolygon = DataNodeType::New(); + DataNodePointerType newDataNodeMultiLine = DataNodeType::New(); + DataNodePointerType newDataNodeMultiFeature = DataNodeType::New(); + /** Walking trough the input vector data */ typedef itk::PreOrderTreeIterator<DataTreeType> TreeIteratorType; TreeIteratorType it(input->GetDataTree()); it.GoToBegin(); + bool newFolder = false ; + bool newMultiFeature = false; + + while (!it.IsAtEnd()) { - DataNodePointerType dataNode = it.Get(); - DataNodePointerType newDataNode = DataNodeType::New(); - newDataNode->SetNodeType(dataNode->GetNodeType()); - newDataNode->SetNodeId(dataNode->GetNodeId()); + DataNodePointerType dataNode = it.Get(); + DataNodePointerType newDataNode = DataNodeType::New(); + + switch (dataNode->GetNodeType()) { case ROOT: { + newDataNode->SetNodeType(dataNode->GetNodeType()); + newDataNode->SetNodeId(dataNode->GetNodeId()); tree->SetRoot(newDataNode); currentContainer = newDataNode; break; } case DOCUMENT: { + newDataNode->SetNodeType(dataNode->GetNodeType()); + newDataNode->SetNodeId(dataNode->GetNodeId()); tree->Add(newDataNode,currentContainer); currentContainer = newDataNode; break; } case FOLDER: { - tree->Add(newDataNode,currentContainer); - currentContainer = newDataNode; + newDataNodeFolder->SetNodeType(dataNode->GetNodeType()); + newDataNodeFolder->SetNodeId(dataNode->GetNodeId()); + newFolder = true; break; } case FEATURE_POINT: { - newDataNode->SetPoint(dataNode->GetPoint()); - tree->Add(newDataNode,currentContainer); + if(m_GeoROI.IsInside(this->PointToContinuousIndex(dataNode->GetPoint()))) + { + if(newFolder) + { + tree->Add(newDataNodeFolder,currentContainer); + currentContainer = newDataNodeFolder; + newFolder = false; + } + if(newMultiFeature) + { + tree->Add(newDataNodeMultiFeature,currentContainer); + currentContainer = newDataNodeMultiFeature; + newMultiFeature = false; + } + newDataNode->SetNodeType(dataNode->GetNodeType()); + newDataNode->SetNodeId(dataNode->GetNodeId()); + newDataNode->SetPoint(dataNode->GetPoint()); + tree->Add(newDataNode,currentContainer); + } + break; } case FEATURE_LINE: { if(this->IsLineIntersectionNotNull(dataNode->GetLine())) { + if(newFolder) + { + tree->Add(newDataNodeFolder,currentContainer); + currentContainer = newDataNodeFolder; + newFolder = false; + } + if(newMultiFeature) + { + tree->Add(newDataNodeMultiFeature ,currentContainer); + currentContainer = newDataNodeMultiFeature ; + newMultiFeature = false; + } + newDataNode->SetNodeType(dataNode->GetNodeType()); + newDataNode->SetNodeId(dataNode->GetNodeId()); newDataNode->SetLine(dataNode->GetLine()); tree->Add(newDataNode,currentContainer); } @@ -139,35 +188,55 @@ VectorDataExtractROI<TVectorData> { if(this->IsPolygonIntersectionNotNull(dataNode->GetPolygonExteriorRing())) { + if(newFolder) + { + tree->Add(newDataNodeFolder,currentContainer); + currentContainer = newDataNodeFolder; + newFolder = false; + } + + if(newMultiFeature) + { + tree->Add(newDataNodeMultiFeature,currentContainer); + currentContainer = newDataNodeMultiFeature ; + newMultiFeature = false; + } + + newDataNode->SetNodeType(dataNode->GetNodeType()); + newDataNode->SetNodeId(dataNode->GetNodeId()); newDataNode->SetPolygonExteriorRing(dataNode->GetPolygonExteriorRing()); newDataNode->SetPolygonInteriorRings(dataNode->GetPolygonInteriorRings()); tree->Add(newDataNode,currentContainer); } - else - std::cout << " OUTSIDE The region" <<std::endl; - break; } case FEATURE_MULTIPOINT: { - tree->Add(newDataNode,currentContainer); - currentContainer = newDataNode; + newDataNodeMultiFeature->SetNodeType(dataNode->GetNodeType()); + newDataNodeMultiFeature ->SetNodeId(dataNode->GetNodeId()); + newMultiFeature = true; + + break; } case FEATURE_MULTILINE: { - tree->Add(newDataNode,currentContainer); - currentContainer = newDataNode; + newDataNodeMultiFeature ->SetNodeType(dataNode->GetNodeType()); + newDataNodeMultiFeature ->SetNodeId(dataNode->GetNodeId()); + newMultiFeature = true; break; } case FEATURE_MULTIPOLYGON: { - tree->Add(newDataNode,currentContainer); - currentContainer = newDataNode; + newDataNodeMultiFeature ->SetNodeType(dataNode->GetNodeType()); + newDataNodeMultiFeature->SetNodeId(dataNode->GetNodeId()); + newMultiFeature = true; break; } case FEATURE_COLLECTION: { + newDataNode->SetNodeType(dataNode->GetNodeType()); + newDataNode->SetNodeId(dataNode->GetNodeId()); tree->Add(newDataNode,currentContainer); currentContainer = newDataNode; break; @@ -211,10 +280,9 @@ void VectorDataExtractROI<TVectorData> ::CompareInputAndRegionProjection() { - /**Traces*/ std::string regionProjection = m_ROI.GetRegionProjection(); std::string inputVectorProjection = this->GetInput()->GetProjectionRef(); - + if(regionProjection == inputVectorProjection) m_ProjectionNeeded = false; else @@ -229,10 +297,10 @@ void VectorDataExtractROI<TVectorData> ::ProjectRegionToInputVectorProjection() { - typedef otb::GenericMapProjection<otb::FORWARD> ForwardMapProjectionType; + typedef otb::GenericMapProjection<otb::INVERSE> ForwardMapProjectionType; ForwardMapProjectionType::Pointer mapTransform = ForwardMapProjectionType::New(); mapTransform->SetWkt(m_ROI.GetRegionProjection()); - + /*FORWARD ; From RegionProjection To long/lat Projection */ typename VertexListType::Pointer regionCorners = VertexListType::New(); ProjPointType point1, point2 , point3, point4; @@ -255,30 +323,39 @@ VectorDataExtractROI<TVectorData> ProjPointType pGeo2 = mapTransform->TransformPoint(point2); ProjPointType pGeo3 = mapTransform->TransformPoint(point3); ProjPointType pGeo4 = mapTransform->TransformPoint(point4); - - /** Inverse : From long/lat to InputVectorData projection*/ - typedef otb::GenericMapProjection<otb::INVERSE> InverseMapProjectionType; + + /** INVERSE : From long/lat to InputVectorData projection*/ + typedef otb::GenericMapProjection<otb::FORWARD> InverseMapProjectionType; InverseMapProjectionType::Pointer mapInverseTransform = InverseMapProjectionType::New(); - if(this->GetInput()->GetProjectionRef().empty()) + + typedef itk::Transform<double, 2, 2> GenericTransformType; + GenericTransformType::Pointer outputTransform ; + + if(!this->GetInput()->GetProjectionRef().empty()) { - std::string inputProjectionRef = "GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]]"; - mapInverseTransform->SetWkt(inputProjectionRef); + mapInverseTransform->SetWkt(this->GetInput()->GetProjectionRef()); + outputTransform = mapInverseTransform.GetPointer(); + + //Case of kml + if(mapInverseTransform->GetMapProjection() == NULL) + outputTransform = itk::IdentityTransform< double, 2 >::New(); } else - mapInverseTransform->SetWkt(this->GetInput()->GetProjectionRef()); + outputTransform = itk::IdentityTransform< double, 2 >::New(); + /** Finally Project the four corners in the InputDataVector Projection*/ - ProjPointType pCarto1 = mapInverseTransform->TransformPoint(pGeo1); - ProjPointType pCarto2 = mapInverseTransform->TransformPoint(pGeo2); - ProjPointType pCarto3 = mapInverseTransform->TransformPoint(pGeo3); - ProjPointType pCarto4 = mapInverseTransform->TransformPoint(pGeo4); - + ProjPointType pCarto1 = /*mapInverseTransform*/ outputTransform->TransformPoint(pGeo1); + ProjPointType pCarto2 = /*mapInverseTransform*/ outputTransform->TransformPoint(pGeo2); + ProjPointType pCarto3 = /*mapInverseTransform*/ outputTransform ->TransformPoint(pGeo3); + ProjPointType pCarto4 = /*mapInverseTransform*/ outputTransform->TransformPoint(pGeo4); + /** Fill the vertex List : First Convert Point To*/ regionCorners->InsertElement(regionCorners->Size(),this->PointToContinuousIndex(pCarto1)); regionCorners->InsertElement(regionCorners->Size(),this->PointToContinuousIndex(pCarto2)); regionCorners->InsertElement(regionCorners->Size(),this->PointToContinuousIndex(pCarto3)); regionCorners->InsertElement(regionCorners->Size(),this->PointToContinuousIndex(pCarto4)); - + /** Due to The projection : the Projected ROI can be rotated */ m_GeoROI = this->ComputeVertexListBoudingRegion(regionCorners.GetPointer()); @@ -292,11 +369,12 @@ typename VectorDataExtractROI<TVectorData> VectorDataExtractROI<TVectorData> ::PointToContinuousIndex(ProjPointType point) { + VertexType vertex; vertex[0] = point[0]; vertex[1] = point[1]; - + return vertex; } @@ -326,13 +404,13 @@ VectorDataExtractROI<TVectorData> y = static_cast<double>(it.Value()[1]); index[0] = x; index[1] = y; - + ++it; while (it != vertexlist->End()) { x = static_cast<double>(it.Value()[0]); y = static_cast<double>(it.Value()[1]); - + // Index search if ( x < index[0] ) { diff --git a/Code/FeatureExtraction/otbLineSpatialObjectListToRightAnglePointSetFilter.txx b/Code/FeatureExtraction/otbLineSpatialObjectListToRightAnglePointSetFilter.txx index 6e78364cc11ce7929c4067de93a86059c4371ff6..79d3e3988b7cdaf04b0c06c416459311fde141d7 100644 --- a/Code/FeatureExtraction/otbLineSpatialObjectListToRightAnglePointSetFilter.txx +++ b/Code/FeatureExtraction/otbLineSpatialObjectListToRightAnglePointSetFilter.txx @@ -189,7 +189,7 @@ LineSpatialObjectListToRightAnglePointSetFilter<TImage,TLinesList ,TPointSet> * d(P,[Q1Q2]) = |PQ1^PQ2|/|Q1Q2| */ - float SegmentLength = (Xq1-Xq2)* (Xq1-Xq2) + (Yq1-Yq2) *(Yq1-Yq2); + float SegmentLength = vcl_sqrt((Xq1-Xq2)* (Xq1-Xq2) + (Yq1-Yq2) *(Yq1-Yq2)); float CrossProduct = Xq1*Yq2 - Xq2*Yq1 ; /** Define a line iterator */ diff --git a/Code/FeatureExtraction/otbSFSTexturesImageFilter.h b/Code/FeatureExtraction/otbSFSTexturesImageFilter.h index f381f5a548dd3ac99c73e2304d40341c9b96c8f2..faa9f3e7c6c0bc2fafeff19a24f9dbe021d1af3d 100644 --- a/Code/FeatureExtraction/otbSFSTexturesImageFilter.h +++ b/Code/FeatureExtraction/otbSFSTexturesImageFilter.h @@ -162,7 +162,8 @@ public: * 6: SD * Set to 1 means the texture will be computed. **/ - typedef enum {LENGTH, WIDTH, PSI, WMEAN, RATIO, SD} FeatureType; + typedef enum {LENGTH=1, WIDTH, PSI, WMEAN, RATIO, SD} FeatureType; + void SetFeatureStatus(FeatureType id, bool isSelected ) { if ( static_cast<unsigned int>(id) > this->GetTexturesStatus().size() || id == 0 ) @@ -180,16 +181,7 @@ public: return this->GetFunctor().GetTexturesStatus(); } - void InitTextureStatusFalse() - { - unsigned int id; - for (id=0;id<=1;id++) - { - this->SetFeatureStatus(static_cast<FeatureType>(id),false); - } - } - - + void InitFeatureStatus(bool status); /** Return output length image */ const OutputImageType * GetLengthOutput() const; diff --git a/Code/FeatureExtraction/otbSFSTexturesImageFilter.txx b/Code/FeatureExtraction/otbSFSTexturesImageFilter.txx index a04fbb48bdb3db5f66f3a427c3af785ce2a4c1d3..91993097e72e835ae1c5fd3b81dc63bbb2d735c0 100644 --- a/Code/FeatureExtraction/otbSFSTexturesImageFilter.txx +++ b/Code/FeatureExtraction/otbSFSTexturesImageFilter.txx @@ -46,6 +46,7 @@ SFSTexturesImageFilter<TInputImage,TOutputImage> m_Radius = this->GetSpatialThreshold(); m_FunctorList.clear(); + } /************************************************************ * @@ -264,6 +265,7 @@ SFSTexturesImageFilter<TInputImage, TOutputImage> { m_FunctorList.push_back(m_Functor); } + this->InitFeatureStatus(true); } template <class TInputImage, class TOutputImage> @@ -424,6 +426,19 @@ SFSTexturesImageFilter<TInputImage, TOutputImage> } } +template <class TInputImage, class TOutputImage> +void +SFSTexturesImageFilter<TInputImage, TOutputImage> +::InitFeatureStatus(bool status) + { + for (FeatureType id=LENGTH;id<=SD; + id=static_cast<FeatureType>(id+1)) + { + this->SetFeatureStatus(static_cast<FeatureType>(id),status); + } + } + + /** * Standard "PrintSelf" method diff --git a/Testing/Code/BasicFilters/otbAmplitudePhaseToRGBFunctor.cxx b/Testing/Code/BasicFilters/otbAmplitudePhaseToRGBFunctor.cxx index 0c32d4d5692ed5eb05c32d9d54e1a0630433e683..d351b5b453003e73860fd0972929da62388a899e 100644 --- a/Testing/Code/BasicFilters/otbAmplitudePhaseToRGBFunctor.cxx +++ b/Testing/Code/BasicFilters/otbAmplitudePhaseToRGBFunctor.cxx @@ -30,10 +30,11 @@ #include "otbImageFileReader.h" #include "otbStreamingImageFileWriter.h" -#include "itkBinaryFunctorImageFilter.h" +#include "itkTernaryFunctorImageFilter.h" #include "otbAmplitudePhaseToRGBFunctor.h" #include "itkComplexToModulusImageFilter.h" #include "itkComplexToPhaseImageFilter.h" +#include "itkShiftScaleImageFilter.h" int otbAmplitudePhaseToRGBFunctor(int argc, char * argv[]) { @@ -64,18 +65,23 @@ int otbAmplitudePhaseToRGBFunctor(int argc, char * argv[]) PhaseFilterType::Pointer phaseFilter = PhaseFilterType::New(); phaseFilter->SetInput(reader->GetOutput()); + typedef itk::ShiftScaleImageFilter<ImageType, ImageType> ConstFilterType; + ConstFilterType::Pointer constFilter = ConstFilterType::New(); + constFilter->SetScale(0.0); + constFilter->SetInput(modulusFilter->GetOutput()); typedef otb::Functor::AmplitudePhaseToRGBFunctor - <PixelType,PixelType,RGBPixelType> ColorMapFunctorType; - typedef itk::BinaryFunctorImageFilter - <ImageType, ImageType, RGBImageType, ColorMapFunctorType> ColorMapFilterType; + <PixelType,PixelType,PixelType,RGBPixelType> ColorMapFunctorType; + typedef itk::TernaryFunctorImageFilter + <ImageType, ImageType, ImageType, RGBImageType, ColorMapFunctorType> ColorMapFilterType; ColorMapFilterType::Pointer colormapper = ColorMapFilterType::New(); colormapper->GetFunctor().SetMaximum(4000); colormapper->GetFunctor().SetMinimum(0); colormapper->SetInput1(modulusFilter->GetOutput()); - colormapper->SetInput2(phaseFilter->GetOutput()); + colormapper->SetInput2(constFilter->GetOutput()); + colormapper->SetInput3(phaseFilter->GetOutput()); colormapper->SetNumberOfThreads(1); writer->SetInput(colormapper->GetOutput()); diff --git a/Testing/Code/BasicFilters/otbAmplitudePhaseToRGBFunctorNew.cxx b/Testing/Code/BasicFilters/otbAmplitudePhaseToRGBFunctorNew.cxx index 6c6fc24c05c604f6b98f3365db40cbe72b585a20..9532c4ce4c336feb175a9a904538cf0348991a01 100644 --- a/Testing/Code/BasicFilters/otbAmplitudePhaseToRGBFunctorNew.cxx +++ b/Testing/Code/BasicFilters/otbAmplitudePhaseToRGBFunctorNew.cxx @@ -28,7 +28,7 @@ #include "otbImage.h" -#include "itkBinaryFunctorImageFilter.h" +#include "itkTernaryFunctorImageFilter.h" #include "otbAmplitudePhaseToRGBFunctor.h" int otbAmplitudePhaseToRGBFunctorNew(int argc, char * argv[]) @@ -41,9 +41,9 @@ int otbAmplitudePhaseToRGBFunctorNew(int argc, char * argv[]) typedef otb::Image<RGBPixelType, 2> RGBImageType; typedef otb::Functor::AmplitudePhaseToRGBFunctor - <PixelType,PixelType,RGBPixelType> ColorMapFunctorType; - typedef itk::BinaryFunctorImageFilter - <ImageType, ImageType, RGBImageType, ColorMapFunctorType> ColorMapFilterType; + <PixelType,PixelType,PixelType,RGBPixelType> ColorMapFunctorType; + typedef itk::TernaryFunctorImageFilter + <ImageType, ImageType, ImageType, RGBImageType, ColorMapFunctorType> ColorMapFilterType; ColorMapFilterType::Pointer colormapper = ColorMapFilterType::New(); colormapper->GetFunctor().SetMaximum(150); colormapper->GetFunctor().SetMinimum(70); diff --git a/Testing/Code/Common/CMakeLists.txt b/Testing/Code/Common/CMakeLists.txt index 0cf97e9d779632035cd39f50e7cf81e2d4202c87..7170f7c439d2956f659d76a4d5828355f0292537 100644 --- a/Testing/Code/Common/CMakeLists.txt +++ b/Testing/Code/Common/CMakeLists.txt @@ -138,7 +138,8 @@ ADD_TEST(coTvVectorDataExtractROI ${COMMON_TESTS2} otbVectorDataExtractROI ${INPUTDATA}/ToulousePoints-examples.shp #${LARGEDATA}/TOULOUSE/QuickBird/GIS_FILES/000000128955_01_ORDER_SHAPE.shp ${TEMP}/coVectorDataExtractROIOutput.shp - 1.588 48.2544 1.2222 1.33001 + 374369.48850211215904 4828951.58612491376698 # Origin of the CartoRegion + 1000.25 25000.2 # Size of the Cartoregion ) # ------- otb::MultiChannelExtractROI ------------------------------ diff --git a/Utilities/ITK/Code/BasicFilters/itkResampleImageFilter.h b/Utilities/ITK/Code/BasicFilters/itkResampleImageFilter.h index 1e23b3de88986b24a8e9bfdad152ec0702e00ead..bf8b206f202e6c14d834cd58d914fa3156a7ab93 100644 --- a/Utilities/ITK/Code/BasicFilters/itkResampleImageFilter.h +++ b/Utilities/ITK/Code/BasicFilters/itkResampleImageFilter.h @@ -9,8 +9,8 @@ Copyright (c) Insight Software Consortium. All rights reserved. See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. - This software is distributed WITHOUT ANY WARRANTY; without even - the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ @@ -69,7 +69,7 @@ namespace itk * ProcessObject::GenerateInputRequestedRegion() and * ProcessObject::GenerateOutputInformation(). * - * This filter is implemented as a multithreaded filter. It provides a + * This filter is implemented as a multithreaded filter. It provides a * ThreadedGenerateData() method for its implementation. * * \ingroup GeometricTransforms @@ -94,7 +94,7 @@ public: typedef typename InputImageType::RegionType InputImageRegionType; /** Method for creation through the object factory. */ - itkNewMacro(Self); + itkNewMacro(Self); /** Run-time type information (and related methods). */ itkTypeMacro(ResampleImageFilter, ImageToImageFilter); @@ -107,13 +107,14 @@ public: /** Transform typedef. */ - typedef Transform<TInterpolatorPrecisionType, - itkGetStaticConstMacro(ImageDimension), + typedef double CoordRepType; + typedef Transform<CoordRepType, + itkGetStaticConstMacro(ImageDimension), itkGetStaticConstMacro(ImageDimension)> TransformType; typedef typename TransformType::ConstPointer TransformPointerType; /** Interpolator typedef. */ - typedef InterpolateImageFunction<InputImageType, TInterpolatorPrecisionType> InterpolatorType; + typedef InterpolateImageFunction<InputImageType,CoordRepType > InterpolatorType; typedef typename InterpolatorType::Pointer InterpolatorPointerType; /** Image size typedef. */ @@ -129,7 +130,7 @@ public: /** Image pixel value typedef. */ typedef typename TOutputImage::PixelType PixelType; typedef typename TInputImage::PixelType InputPixelType; - + /** Typedef to describe the output image region type. */ typedef typename TOutputImage::RegionType OutputImageRegionType; @@ -137,7 +138,7 @@ public: typedef typename TOutputImage::SpacingType SpacingType; typedef typename TOutputImage::PointType OriginPointType; typedef typename TOutputImage::DirectionType DirectionType; - + /** Set the coordinate transformation. * Set the coordinate transform to use for resampling. Note that this must * be in physical coordinates and it is the output-to-input transform, NOT @@ -145,7 +146,7 @@ public: * the filter uses an Identity transform. You must provide a different * transform here, before attempting to run the filter, if you do not want to * use the default Identity transform. */ - itkSetConstObjectMacro( Transform, TransformType ); + itkSetConstObjectMacro( Transform, TransformType ); /** Get a pointer to the coordinate transform. */ itkGetConstObjectMacro( Transform, TransformType ); @@ -166,7 +167,7 @@ public: /** Get the size of the output image. */ itkGetConstReferenceMacro( Size, SizeType ); - + /** Set the pixel value when a transformed pixel is outside of the * image. The default default pixel value is 0. */ itkSetMacro( DefaultPixelValue, PixelType ); @@ -211,7 +212,7 @@ public: this->SetOutputStartIndex ( Image->GetLargestPossibleRegion().GetIndex() ); this->SetSize ( Image->GetLargestPossibleRegion().GetSize() ); } - /** Set the start index of the output largest possible region. + /** Set the start index of the output largest possible region. * The default is an index of all zeros. */ itkSetMacro( OutputStartIndex, IndexType ); @@ -222,7 +223,7 @@ public: * the information is specified with the SetOutputSpacing, Origin, * and Direction methods. UseReferenceImage must be On and a * Reference image must be present to override the defaul behavior. - * NOTE: This function seems redundant with the + * NOTE: This function seems redundant with the * SetOutputParametersFromImage( image ) function */ void SetReferenceImage ( const TOutputImage *image ); const TOutputImage * GetReferenceImage( void ) const; @@ -245,11 +246,11 @@ public: * \sa ProcessObject::GenerateInputRequestedRegion() */ virtual void GenerateInputRequestedRegion( void ); - /** This method is used to set the state of the filter before + /** This method is used to set the state of the filter before * multi-threading. */ virtual void BeforeThreadedGenerateData( void ); - /** This method is used to set the state of the filter after + /** This method is used to set the state of the filter after * multi-threading. */ virtual void AfterThreadedGenerateData( void ); @@ -286,12 +287,12 @@ protected: int threadId ); /** Implementation for resampling that works for with linear - * transformation types. + * transformation types. */ void LinearThreadedGenerateData( const OutputImageRegionType& outputRegionForThread, int threadId ); - + private: ResampleImageFilter( const Self& ); //purposely not implemented @@ -312,13 +313,13 @@ private: }; - + } // end namespace itk - + #ifndef ITK_MANUAL_INSTANTIATION #include "itkResampleImageFilter.txx" #endif - + #endif - + #endif diff --git a/Utilities/ITK/Code/BasicFilters/itkResampleImageFilter.txx b/Utilities/ITK/Code/BasicFilters/itkResampleImageFilter.txx index d7c93968994a33fbd22ce87dde23d2fef3958729..17777312d88078bfdcee002d7167ae49e4b57a6b 100644 --- a/Utilities/ITK/Code/BasicFilters/itkResampleImageFilter.txx +++ b/Utilities/ITK/Code/BasicFilters/itkResampleImageFilter.txx @@ -9,8 +9,8 @@ Copyright (c) Insight Software Consortium. All rights reserved. See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. - This software is distributed WITHOUT ANY WARRANTY; without even - the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ @@ -54,9 +54,10 @@ ResampleImageFilter<TInputImage, TOutputImage,TInterpolatorPrecisionType> m_Size.Fill( 0 ); m_OutputStartIndex.Fill( 0 ); - - m_Transform = IdentityTransform<TInterpolatorPrecisionType, ImageDimension>::New(); - m_Interpolator = LinearInterpolateImageFunction<InputImageType, TInterpolatorPrecisionType>::New(); + + typedef double CoordRepType; + m_Transform = IdentityTransform<CoordRepType, ImageDimension>::New(); + m_Interpolator = LinearInterpolateImageFunction<InputImageType, CoordRepType>::New(); m_DefaultPixelValue = 0; } @@ -67,12 +68,12 @@ ResampleImageFilter<TInputImage, TOutputImage,TInterpolatorPrecisionType> * \todo Add details about this class */ template <class TInputImage, class TOutputImage, class TInterpolatorPrecisionType> -void +void ResampleImageFilter<TInputImage, TOutputImage,TInterpolatorPrecisionType> ::PrintSelf(std::ostream& os, Indent indent) const { Superclass::PrintSelf(os,indent); - + os << indent << "DefaultPixelValue: " << static_cast<typename NumericTraits<PixelType>::PrintType>(m_DefaultPixelValue) << std::endl; @@ -91,7 +92,7 @@ ResampleImageFilter<TInputImage, TOutputImage,TInterpolatorPrecisionType> * Set the output image spacing. */ template <class TInputImage, class TOutputImage, class TInterpolatorPrecisionType> -void +void ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> ::SetOutputSpacing( const double* spacing) @@ -105,7 +106,7 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> * Set the output image origin. */ template <class TInputImage, class TOutputImage, class TInterpolatorPrecisionType> -void +void ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> ::SetOutputOrigin( const double* origin) @@ -121,7 +122,7 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> * has to be set up before ThreadedGenerateData */ template <class TInputImage, class TOutputImage, class TInterpolatorPrecisionType> -void +void ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> ::BeforeThreadedGenerateData() { @@ -144,7 +145,7 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> * Set up state of filter after multi-threading. */ template <class TInputImage, class TOutputImage, class TInterpolatorPrecisionType> -void +void ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> ::AfterThreadedGenerateData() { @@ -157,7 +158,7 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> * ThreadedGenerateData */ template <class TInputImage, class TOutputImage, class TInterpolatorPrecisionType> -void +void ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> ::ThreadedGenerateData( const OutputImageRegionType& outputRegionForThread, @@ -183,7 +184,7 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> this->NonlinearThreadedGenerateData(outputRegionForThread, threadId); return; } - + // Check whether we can use a fast path for resampling. Fast path // can be used if the transformation is linear. Transform respond // to the IsLinear() call. @@ -196,12 +197,12 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> // Otherwise, we use the normal method where the transform is called // for computing the transformation of every point. this->NonlinearThreadedGenerateData(outputRegionForThread, threadId); - + } template <class TInputImage, class TOutputImage, class TInterpolatorPrecisionType> -void +void ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> ::NonlinearThreadedGenerateData( const OutputImageRegionType& outputRegionForThread, @@ -223,12 +224,13 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> PointType outputPoint; // Coordinates of current output pixel PointType inputPoint; // Coordinates of current input pixel - typedef ContinuousIndex<TInterpolatorPrecisionType, ImageDimension> ContinuousIndexType; + typedef double CoordRepType; + typedef ContinuousIndex<CoordRepType, ImageDimension> ContinuousIndexType; ContinuousIndexType inputIndex; // Support for progress methods/callbacks ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); - + typedef typename InterpolatorType::OutputType OutputType; // Min/max values of the output pixel type AND these values @@ -238,7 +240,7 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> const OutputType minOutputValue = static_cast<OutputType>(minValue); const OutputType maxOutputValue = static_cast<OutputType>(maxValue); - + // Walk the output region outIt.GoToBegin(); @@ -258,7 +260,7 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> inputPtr->TransformPhysicalPointToContinuousIndex(inputPoint, inputIndex); // The inputIndex is precise to many decimal points, but this precision - // involves some error in the last bits. + // involves some error in the last bits. // Sometimes, when an index should be inside of the image, the // index will be slightly // greater than the largest index in the image, like 255.00000000002 @@ -274,25 +276,16 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> double newInputIndexFrac = vcl_floor(precisionConstant * inputIndexFrac)/precisionConstant; inputIndex[i] = roundedInputIndex + newInputIndexFrac; } - + // Evaluate input at right position and copy to the output if( m_Interpolator->IsInsideBuffer(inputIndex) ) { PixelType pixval; const OutputType value = m_Interpolator->EvaluateAtContinuousIndex(inputIndex); - if( value < minOutputValue ) - { - pixval = minValue; - } - else if( value > maxOutputValue ) - { - pixval = maxValue; - } - else - { - pixval = static_cast<PixelType>( value ); - } + pixval = static_cast<PixelType>( + NumericTraits<OutputType>::Clamp(value,minOutputValue,maxOutputValue) + ); outIt.Set( pixval ); } else @@ -308,7 +301,7 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> } template <class TInputImage, class TOutputImage, class TInterpolatorPrecisionType> -void +void ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> ::LinearThreadedGenerateData( const OutputImageRegionType& outputRegionForThread, @@ -333,17 +326,18 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> PointType tmpOutputPoint; PointType tmpInputPoint; - typedef ContinuousIndex<TInterpolatorPrecisionType, ImageDimension> ContinuousIndexType; + typedef double CoordRepType; + typedef ContinuousIndex<CoordRepType, ImageDimension> ContinuousIndexType; ContinuousIndexType inputIndex; ContinuousIndexType tmpInputIndex; - + typedef typename PointType::VectorType VectorType; VectorType delta; // delta in input continuous index coordinate frame IndexType index; // Support for progress methods/callbacks ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); - + typedef typename InterpolatorType::OutputType OutputType; // Cache information from the superclass @@ -356,34 +350,34 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> const OutputType minOutputValue = static_cast<OutputType>(minValue); const OutputType maxOutputValue = static_cast<OutputType>(maxValue); - - + + // Determine the position of the first pixel in the scanline index = outIt.GetIndex(); outputPtr->TransformIndexToPhysicalPoint( index, outputPoint ); - + // Compute corresponding input pixel position inputPoint = m_Transform->TransformPoint(outputPoint); inputPtr->TransformPhysicalPointToContinuousIndex(inputPoint, inputIndex); - + // As we walk across a scan line in the output image, we trace // an oriented/scaled/translated line in the input image. Cache // the delta along this line in continuous index space of the input // image. This allows us to use vector addition to model the // transformation. // - // To find delta, we take two pixels adjacent in a scanline - // and determine the continuous indices of these pixels when + // To find delta, we take two pixels adjacent in a scanline + // and determine the continuous indices of these pixels when // mapped to the input coordinate frame. We use the difference // between these two continuous indices as the delta to apply - // to an index to trace line in the input image as we move + // to an index to trace line in the input image as we move // across the scanline of the output image. // // We determine delta in this manner so that Images and // OrientedImages are both handled properly (with the delta for // OrientedImages taking into account the direction cosines). - // + // ++index[0]; outputPtr->TransformIndexToPhysicalPoint( index, tmpOutputPoint ); tmpInputPoint = m_Transform->TransformPoint( tmpOutputPoint ); @@ -432,7 +426,7 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> // will incremented in the scanline loop inputPoint = m_Transform->TransformPoint(outputPoint); inputPtr->TransformPhysicalPointToContinuousIndex(inputPoint, inputIndex); - + // The inputIndex is precise to many decimal points, but this precision // involves some error in the last bits. // Sometimes, when an index should be inside of the image, the @@ -460,25 +454,17 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> PixelType pixval; const OutputType value = m_Interpolator->EvaluateAtContinuousIndex(inputIndex); - if( value < minOutputValue ) - { - pixval = minValue; - } - else if( value > maxOutputValue ) - { - pixval = maxValue; - } - else - { - pixval = static_cast<PixelType>( value ); - } + + pixval = static_cast<PixelType>( + NumericTraits<OutputType>::Clamp(value,minOutputValue,maxOutputValue) + ); outIt.Set( pixval ); } else { outIt.Set(defaultValue); // default background value } - + progress.CompletedPixel(); ++outIt; inputIndex += delta; @@ -490,7 +476,7 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> } -/** +/** * Inform pipeline of necessary input image region * * Determining the actual input region is non-trivial, especially @@ -498,7 +484,7 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> * So we do the easy thing and request the entire input image. */ template <class TInputImage, class TOutputImage, class TInterpolatorPrecisionType> -void +void ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> ::GenerateInputRequestedRegion() { @@ -511,7 +497,7 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> } // get pointers to the input and output - InputImagePointer inputPtr = + InputImagePointer inputPtr = const_cast< TInputImage *>( this->GetInput() ); // Request the entire input image @@ -523,7 +509,7 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> } -/** +/** * Set the smart pointer to the reference image that will provide * the grid parameters for the output image. */ @@ -533,18 +519,18 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> ::GetReferenceImage() const { Self * surrogate = const_cast< Self * >( this ); - const OutputImageType * referenceImage = + const OutputImageType * referenceImage = static_cast<const OutputImageType *>(surrogate->ProcessObject::GetInput(1)); return referenceImage; } -/** +/** * Set the smart pointer to the reference image that will provide * the grid parameters for the output image. */ template <class TInputImage, class TOutputImage, class TInterpolatorPrecisionType> -void +void ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> ::SetReferenceImage( const TOutputImage *image ) { @@ -556,11 +542,11 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> } } -/** +/** * Inform pipeline of required output region */ template <class TInputImage, class TOutputImage, class TInterpolatorPrecisionType> -void +void ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> ::GenerateOutputInformation() { @@ -605,15 +591,15 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> return; } -/** +/** * Verify if any of the components has been modified. */ template <class TInputImage, class TOutputImage, class TInterpolatorPrecisionType> -unsigned long +unsigned long ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> ::GetMTime( void ) const { - unsigned long latestTime = Object::GetMTime(); + unsigned long latestTime = Object::GetMTime(); if( m_Transform ) { diff --git a/Utilities/ITK/Code/Common/itkNumericTraits.h b/Utilities/ITK/Code/Common/itkNumericTraits.h index f0d95a38addd52526dcc400be81cbae95502fcab..8691ed1112f06c98b9ae39507dfd85aa096b55c0 100644 --- a/Utilities/ITK/Code/Common/itkNumericTraits.h +++ b/Utilities/ITK/Code/Common/itkNumericTraits.h @@ -183,6 +183,9 @@ public: static bool IsNonnegative(char val) {return val >= Zero; } static char ZeroValue() { return Zero; } static char OneValue() { return One; } + static ValueType Clamp(ValueType val,ValueType minVal, ValueType maxVal) + {return val<minVal?minVal:(val>maxVal?maxVal:val);} + }; /** \class NumericTraits<char> @@ -214,6 +217,9 @@ public: static bool IsNonnegative(signed char val) {return val >= Zero; } static signed char ZeroValue() { return Zero; } static signed char OneValue() { return One; } + static ValueType Clamp(ValueType val,ValueType minVal, ValueType maxVal) + {return val<minVal?minVal:(val>maxVal?maxVal:val);} + }; /** \class NumericTraits<unsigned char> @@ -243,6 +249,9 @@ public: static bool IsNonnegative(unsigned char /*val */) {return true; } static unsigned char ZeroValue() { return Zero; } static unsigned char OneValue() { return One; } + static ValueType Clamp(ValueType val,ValueType minVal, ValueType maxVal) + {return val<minVal?minVal:(val>maxVal?maxVal:val);} + }; /** \class NumericTraits<short> @@ -270,6 +279,9 @@ public: static bool IsNonnegative(short val) {return val >= Zero; } static short ZeroValue() { return Zero; } static short OneValue() { return One; } + static ValueType Clamp(ValueType val,ValueType minVal, ValueType maxVal) + {return val<minVal?minVal:(val>maxVal?maxVal:val);} + }; /** \class NumericTraits<unsigned short> @@ -298,6 +310,9 @@ public: static bool IsNonnegative(unsigned short /*val*/) {return true; } static unsigned short ZeroValue() { return Zero; } static unsigned short OneValue() { return One; } + static ValueType Clamp(ValueType val,ValueType minVal, ValueType maxVal) + {return val<minVal?minVal:(val>maxVal?maxVal:val);} + }; /** \class NumericTraits<int> @@ -325,6 +340,9 @@ public: static bool IsNonnegative(int val) {return val >= Zero; } static int ZeroValue() { return Zero; } static int OneValue() { return One; } + static ValueType Clamp(ValueType val,ValueType minVal, ValueType maxVal) + {return val<minVal?minVal:(val>maxVal?maxVal:val);} + }; /** \class NumericTraits<unsigned int> @@ -356,6 +374,9 @@ public: static bool IsNonnegative(unsigned int /*val*/) {return true; } static unsigned int ZeroValue() { return Zero; } static unsigned int OneValue() { return One; } + static ValueType Clamp(ValueType val,ValueType minVal, ValueType maxVal) + {return val<minVal?minVal:(val>maxVal?maxVal:val);} + }; /** \class NumericTraits<long> @@ -384,6 +405,9 @@ public: static bool IsNonnegative(long val) {return val >= Zero; } static long ZeroValue() { return Zero; } static long OneValue() { return One; } + static ValueType Clamp(ValueType val,ValueType minVal, ValueType maxVal) + {return val<minVal?minVal:(val>maxVal?maxVal:val);} + }; /** \class NumericTraits<unsigned long> @@ -412,6 +436,9 @@ public: static bool IsNonnegative(unsigned long) {return true; } static unsigned long ZeroValue() { return Zero; } static unsigned long OneValue() { return One; } + static ValueType Clamp(ValueType val,ValueType minVal, ValueType maxVal) + {return val<minVal?minVal:(val>maxVal?maxVal:val);} + }; /** \class NumericTraits<float> @@ -440,6 +467,9 @@ public: static bool IsNonnegative(float val) {return val >= Zero; } static float ZeroValue() { return Zero; } static float OneValue() { return One; } + static ValueType Clamp(ValueType val,ValueType minVal, ValueType maxVal) + {return val<minVal?minVal:(val>maxVal?maxVal:val);} + }; /** \class NumericTraits<double> @@ -468,6 +498,9 @@ public: static bool IsNonnegative(double val) {return val >= Zero; } static double ZeroValue() { return Zero; } static double OneValue() { return One; } + static ValueType Clamp(ValueType val,ValueType minVal, ValueType maxVal) + {return val<minVal?minVal:(val>maxVal?maxVal:val);} + }; /** \class NumericTraits<long double> @@ -496,6 +529,9 @@ public: static bool IsNonnegative(long double val) {return val >= Zero; } static long double ZeroValue() { return Zero; } static long double OneValue() { return One; } + static ValueType Clamp(ValueType val,ValueType minVal, ValueType maxVal) + {return val<minVal?minVal:(val>maxVal?maxVal:val);} + }; /** \class NumericTraits< std::complex<float> > @@ -517,6 +553,8 @@ public: static const TheType ITKCommon_EXPORT Zero; static const TheType ITKCommon_EXPORT One; + static TheType min() { return vcl_numeric_limits<ValueType>::min(); } + static TheType max() { return vcl_numeric_limits<ValueType>::max(); } static TheType min( TheType ) { return vcl_numeric_limits<ValueType>::min(); } static TheType max( TheType ) { return vcl_numeric_limits<ValueType>::max(); } static TheType NonpositiveMin() { @@ -527,6 +565,12 @@ public: static bool IsNonnegative(TheType val) {return val.real() >= 0.0; } static TheType ZeroValue() { return Zero; } static TheType OneValue() { return One; } + static TheType Clamp(TheType val,TheType minVal, TheType maxVal) + { + return TheType(NumericTraits<ValueType>::Clamp(val.real(),minVal.real(),maxVal.real()), + NumericTraits<ValueType>::Clamp(val.imag(),minVal.imag(),maxVal.imag())); + } + }; @@ -549,6 +593,8 @@ public: static const TheType ITKCommon_EXPORT Zero; static const TheType ITKCommon_EXPORT One; + static TheType min() { return vcl_numeric_limits<ValueType>::min(); } + static TheType max() { return vcl_numeric_limits<ValueType>::max(); } static TheType min( TheType ) { return vcl_numeric_limits<ValueType>::min(); } static TheType max( TheType ) { return vcl_numeric_limits<ValueType>::max(); } static TheType NonpositiveMin() { @@ -559,6 +605,11 @@ public: static bool IsNonnegative(TheType val) {return val.real() >= 0.0; } static TheType ZeroValue() { return Zero; } static TheType OneValue() { return One; } + static TheType Clamp(TheType val,TheType minVal, TheType maxVal) + { + return TheType(NumericTraits<ValueType>::Clamp(val.real(),minVal.real(),maxVal.real()), + NumericTraits<ValueType>::Clamp(val.imag(),minVal.imag(),maxVal.imag())); + } }; diff --git a/Utilities/ITK/Code/Review/itkOptResampleImageFilter.h b/Utilities/ITK/Code/Review/itkOptResampleImageFilter.h index ae1f32d0bd08877df3cb24df9d2f79f1a7517c57..ef306d1d1f96dd0c93ae683d1bf9bab66946bd6a 100755 --- a/Utilities/ITK/Code/Review/itkOptResampleImageFilter.h +++ b/Utilities/ITK/Code/Review/itkOptResampleImageFilter.h @@ -101,26 +101,27 @@ public: TInputImage::ImageDimension); - /** + /** * Transform typedef. */ - typedef Transform<TInterpolatorPrecisionType, + typedef double CoordRepType; + typedef Transform<CoordRepType, itkGetStaticConstMacro(ImageDimension), itkGetStaticConstMacro(ImageDimension)> TransformType; typedef typename TransformType::ConstPointer TransformPointerType; /** Interpolator typedef. */ typedef InterpolateImageFunction<InputImageType, - TInterpolatorPrecisionType> InterpolatorType; + CoordRepType> InterpolatorType; typedef typename InterpolatorType::Pointer InterpolatorPointerType; typedef LinearInterpolateImageFunction<InputImageType, - TInterpolatorPrecisionType> LinearInterpolatorType; + CoordRepType> LinearInterpolatorType; typedef typename LinearInterpolatorType::Pointer LinearInterpolatorPointerType; typedef BSplineInterpolateImageFunction<InputImageType, - TInterpolatorPrecisionType> BSplineInterpolatorType; + CoordRepType> BSplineInterpolatorType; typedef typename BSplineInterpolatorType::Pointer BSplineInterpolatorPointerType; diff --git a/Utilities/ITK/Code/Review/itkOptResampleImageFilter.txx b/Utilities/ITK/Code/Review/itkOptResampleImageFilter.txx index 65aa6a4cec3147361007a7804dbf4c768cc3a518..e51ad8e881ff305f8a78d3ab7942c77e310fdd6a 100755 --- a/Utilities/ITK/Code/Review/itkOptResampleImageFilter.txx +++ b/Utilities/ITK/Code/Review/itkOptResampleImageFilter.txx @@ -49,15 +49,15 @@ ResampleImageFilter<TInputImage, TOutputImage,TInterpolatorPrecisionType> m_Size.Fill( 0 ); m_OutputStartIndex.Fill( 0 ); - - m_Transform = IdentityTransform<TInterpolatorPrecisionType, ImageDimension>::New(); + typedef double CoordRepType; + m_Transform = IdentityTransform<CoordRepType, ImageDimension>::New(); m_InterpolatorIsBSpline = false; m_BSplineInterpolator = NULL; m_InterpolatorIsLinear = true; m_LinearInterpolator = LinearInterpolateImageFunction<InputImageType, - TInterpolatorPrecisionType>::New(); + CoordRepType>::New(); m_Interpolator = dynamic_cast<InterpolatorType *> ( m_LinearInterpolator.GetPointer() ); @@ -273,7 +273,8 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> PointType outputPoint; // Coordinates of current output pixel PointType inputPoint; // Coordinates of current input pixel - typedef ContinuousIndex<TInterpolatorPrecisionType, ImageDimension> + typedef double CoordRepType; + typedef ContinuousIndex<CoordRepType, ImageDimension> ContinuousIndexType; ContinuousIndexType inputIndex; @@ -343,18 +344,10 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> const OutputType value = m_BSplineInterpolator ->EvaluateAtContinuousIndex(inputIndex, threadId); - if( value < minOutputValue ) - { - pixval = minValue; - } - else if( value > maxOutputValue ) - { - pixval = maxValue; - } - else - { - pixval = static_cast<PixelType>( value ); - } + pixval = static_cast<PixelType>( + NumericTraits<OutputType>::Clamp(value,minOutputValue,maxOutputValue) + ); + outIt.Set( pixval ); } else @@ -407,18 +400,10 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> PixelType pixval; const OutputType value = m_LinearInterpolator ->EvaluateAtContinuousIndex(inputIndex); - if( value < minOutputValue ) - { - pixval = minValue; - } - else if( value > maxOutputValue ) - { - pixval = maxValue; - } - else - { - pixval = static_cast<PixelType>( value ); - } + + pixval = static_cast<PixelType>( + NumericTraits<OutputType>::Clamp(value,minOutputValue,maxOutputValue) + ); outIt.Set( pixval ); } else @@ -471,18 +456,10 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> PixelType pixval; const OutputType value = m_Interpolator ->EvaluateAtContinuousIndex(inputIndex); - if( value < minOutputValue ) - { - pixval = minValue; - } - else if( value > maxOutputValue ) - { - pixval = maxValue; - } - else - { - pixval = static_cast<PixelType>( value ); - } + + pixval = static_cast<PixelType>( + NumericTraits<OutputType>::Clamp(value,minOutputValue,maxOutputValue) + ); outIt.Set( pixval ); } else @@ -526,7 +503,8 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> PointType tmpOutputPoint; PointType tmpInputPoint; - typedef ContinuousIndex<TInterpolatorPrecisionType, ImageDimension> + typedef double CoordRepType; + typedef ContinuousIndex<CoordRepType, ImageDimension> ContinuousIndexType; ContinuousIndexType inputIndex; ContinuousIndexType tmpInputIndex; @@ -686,18 +664,10 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> { value = m_Interpolator->EvaluateAtContinuousIndex(inputIndex); } - if( value < minOutputValue ) - { - pixval = minValue; - } - else if( value > maxOutputValue ) - { - pixval = maxValue; - } - else - { - pixval = static_cast<PixelType>( value ); - } + + pixval = static_cast<PixelType>( + NumericTraits<OutputType>::Clamp(value,minOutputValue,maxOutputValue) + ); outIt.Set( pixval ); progress.CompletedPixel(); @@ -727,18 +697,10 @@ ResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType> { value = m_Interpolator->EvaluateAtContinuousIndex(inputIndex); } - if( value < minOutputValue ) - { - pixval = minValue; - } - else if( value > maxOutputValue ) - { - pixval = maxValue; - } - else - { - pixval = static_cast<PixelType>( value ); - } + + pixval = static_cast<PixelType>( + NumericTraits<OutputType>::Clamp(value,minOutputValue,maxOutputValue) + ); outIt.Set( pixval ); } else