Commit aa7100db authored by Julien Malik's avatar Julien Malik
Browse files

WIP: merge Hyperspectral code

parent 89860f88
......@@ -38,6 +38,7 @@ SET(OTB_INCLUDE_DIRS_BUILD_TREE ${OTB_INCLUDE_DIRS_BUILD_TREE}
${OTB_SOURCE_DIR}/Code/Fuzzy
${OTB_SOURCE_DIR}/Code/GeospatialAnalysis
${OTB_SOURCE_DIR}/Code/Gui
${OTB_SOURCE_DIR}/Code/Hyperspectral
${OTB_SOURCE_DIR}/Code/IO
${OTB_SOURCE_DIR}/Code/Learning
${OTB_SOURCE_DIR}/Code/Markov
......@@ -302,6 +303,7 @@ SET(OTB_INCLUDE_RELATIVE_DIRS ${OTB_INCLUDE_RELATIVE_DIRS}
Markov
Fusion
GeospatialAnalysis
Hyperspectral
Testing
UtilitiesAdapters/CurlAdapters
UtilitiesAdapters/OssimAdapters
......
......@@ -19,7 +19,7 @@ ADD_SUBDIRECTORY(ObjectDetection)
ADD_SUBDIRECTORY(Fuzzy)
ADD_SUBDIRECTORY(MultiTemporal)
ADD_SUBDIRECTORY(Simulation)
ADD_SUBDIRECTORY(Hyperspectral)
IF(OTB_USE_VISU_GUI)
ADD_SUBDIRECTORY(Visu)
......
# Sources of non-templated classes.
FILE(GLOB OTBHyperspectral_SRCS "*.cxx" )
ADD_LIBRARY(OTBHyperspectral ${OTBHyperspectral_SRCS})
TARGET_LINK_LIBRARIES (OTBHyperspectral OTBCommon)
IF(OTB_LIBRARY_PROPERTIES)
SET_TARGET_PROPERTIES(OTBHyperspectral PROPERTIES ${OTB_LIBRARY_PROPERTIES})
ENDIF(OTB_LIBRARY_PROPERTIES)
IF(NOT OTB_INSTALL_NO_LIBRARIES)
INSTALL(TARGETS OTBHyperspectral
RUNTIME DESTINATION ${OTB_INSTALL_BIN_DIR_CM24} COMPONENT RuntimeLibraries
LIBRARY DESTINATION ${OTB_INSTALL_LIB_DIR_CM24} COMPONENT RuntimeLibraries
ARCHIVE DESTINATION ${OTB_INSTALL_LIB_DIR_CM24} COMPONENT Development)
ENDIF(NOT OTB_INSTALL_NO_LIBRARIES)
IF(NOT OTB_INSTALL_NO_DEVELOPMENT)
FILE(GLOB __files1 "${CMAKE_CURRENT_SOURCE_DIR}/*.h")
FILE(GLOB __files2 "${CMAKE_CURRENT_SOURCE_DIR}/*.txx")
INSTALL(FILES ${__files1} ${__files2}
DESTINATION ${OTB_INSTALL_INCLUDE_DIR_CM24}/Hyperspectral
COMPONENT Development)
ENDIF(NOT OTB_INSTALL_NO_DEVELOPMENT)
......@@ -1962,6 +1962,92 @@ ADD_TEST(bfTvMaskMuParserFilterTest ${BASICFILTERS_TESTS13}
"(b1>100)*(b2>120)"
)
ADD_TEST(bfTuHorizontalSobelVectorImageFilterNew
${BASICFILTERS_TESTS13}
otbHorizontalSobelVectorImageFilterNewTest)
ADD_TEST(bfTvHorizontalSobelVectorImageFilter
${BASICFILTERS_TESTS13}
--compare-image ${NOTOL}
${BASELINE}/bfTvHorizontalSobelVectorImageFilter.tif
${TEMP}/bfTvHorizontalSobelVectorImageFilter.tif
otbHorizontalSobelVectorImageFilterTest
-in ${INPUTDATA}/cupriteSubHsi.tif
-out ${TEMP}/bfTvHorizontalSobelVectorImageFilter.tif)
ADD_TEST(bfTuVerticalSobelVectorImageFilterNew
${BASICFILTERS_TESTS13}
otbVerticalSobelVectorImageFilterNewTest)
ADD_TEST(bfTvVerticalSobelVectorImageFilter
${BASICFILTERS_TESTS13}
--compare-image ${NOTOL}
${BASELINE}/bfTvVerticalSobelVectorImageFilter.tif
${TEMP}/bfTvVerticalSobelVectorImageFilter.tif
otbVerticalSobelVectorImageFilterTest
-in ${INPUTDATA}/cupriteSubHsi.tif
-out ${TEMP}/bfTvVerticalSobelVectorImageFilter.tif)
ADD_TEST(bfTuSobelVectorImageFilterNew
${BASICFILTERS_TESTS13}
otbSobelVectorImageFilterNewTest)
ADD_TEST(bfTvSobelVectorImageFilter
${BASICFILTERS_TESTS13}
--compare-image ${NOTOL}
${BASELINE}/bfTvSobelVectorImageFilter.tif
${TEMP}/bfTvSobelVectorImageFilter.tif
otbSobelVectorImageFilterTest
-in ${INPUTDATA}/cupriteSubHsi.tif
-out ${TEMP}/bfTvSobelVectorImageFilter.tif)
ADD_TEST(bfTuLocalGradientVectorImageFilterNew
${BASICFILTERS_TESTS13}
otbLocalGradientVectorImageFilterNewTest)
ADD_TEST(bfTvLocalGradientVectorImageFilter
${BASICFILTERS_TESTS13}
--compare-image ${NOTOL}
${BASELINE}/bfTvLocalGradientVectorImageFilter.tif
${TEMP}/bfTvLocalGradientVectorImageFilter.tif
otbLocalGradientVectorImageFilterTest
-in ${INPUTDATA}/cupriteSubHsi.tif
-out ${TEMP}/bfTvLocalGradientVectorImageFilter.tif)
ADD_TEST(bfTuLocalActivityVectorImageFilterNew
${BASICFILTERS_TESTS13}
otbLocalActivityVectorImageFilterNewTest)
ADD_TEST(bfTvLocalActivityVectorImageFilter
${BASICFILTERS_TESTS13}
--compare-image ${NOTOL}
${BASELINE}/bfTvLocalActivityVectorImageFilter.tif
${TEMP}/bfTvLocalActivityVectorImageFilter.tif
otbLocalActivityVectorImageFilterTest
-in ${INPUTDATA}/cupriteSubHsi.tif
-out ${TEMP}/bfTvLocalActivityVectorImageFilter.tif)
ADD_TEST(bfTuNormalizeVectorImageFilter
${BASICFILTERS_TESTS13}
otbNormalizeVectorImageFilterNewTest)
ADD_TEST(bfTvNormalizeVectorImageFilter
${BASICFILTERS_TESTS13}
--compare-image ${NOTOL}
${BASELINE}/bfTvNormalizeVectorImageFilter.tif
${TEMP}/bfTvNormalizeVectorImageFilter.tif
otbNormalizeVectorImageFilterTest
-in ${INPUTDATA}/cupriteSubHsi.tif
-out ${TEMP}/bfTvNormalizeVectorImageFilter.tif)
ADD_TEST(bfTuVectorImageToMatrixNew
${BASICFILTERS_TESTS13}
otbVectorImageToMatrixNewTest)
ADD_TEST(bfTvVectorImageToMatrix
${BASICFILTERS_TESTS13}
otbVectorImageToMatrixTest)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ otbBasicFiltersTests14 ~~~~~~~~~~~~~~~~~~~~~~~~~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
......@@ -2425,6 +2511,18 @@ otbLog10ThresholdedImageFilterTest.cxx
otbMatrixImageFilterTest.cxx
otbMaskMuParserFilterTest.cxx
otbMaskMuParserFilterNew.cxx
otbProjectiveProjection.cxx
otbNormalizeVectorImageFilter.cxx
otbHorizontalSobelVectorImageFilter.cxx
otbVerticalSobelVectorImageFilter.cxx
otbSobelVectorImageFilter.cxx
otbLocalActivityVectorImageFilter.cxx
otbLocalGradientVectorImageFilter.cxx
otbPCAImageFilter.cxx
otbMNFImageFilter.cxx
otbNAPCAImageFilter.cxx
otbFastICAImageFilter.cxx
otbVectorImageToMatrixImageFilter.cxx
)
SET(BasicFilters_SRCS14
......
......@@ -38,4 +38,29 @@ void RegisterTests()
REGISTER_TEST(otbMatrixImageFilterTest);
REGISTER_TEST(otbMaskMuParserFilterNew);
REGISTER_TEST(otbMaskMuParserFilterTest);
REGISTER_TEST(otbProjectiveProjectionTestHighSNR);
REGISTER_TEST(otbPCAImageFilterNewTest);
REGISTER_TEST(otbPCAImageFilterTest);
REGISTER_TEST(otbMNFImageFilterNewTest);
REGISTER_TEST(otbMNFImageFilterTest);
REGISTER_TEST(otbNAPCAImageFilterNewTest);
REGISTER_TEST(otbNAPCAImageFilterTest);
REGISTER_TEST(otbNormalizeVectorImageFilterNewTest);
REGISTER_TEST(otbNormalizeVectorImageFilterTest);
REGISTER_TEST(otbLocalActivityVectorImageFilterNewTest);
REGISTER_TEST(otbLocalActivityVectorImageFilterTest);
REGISTER_TEST(otbLocalGradientVectorImageFilterNewTest);
REGISTER_TEST(otbLocalGradientVectorImageFilterTest);
REGISTER_TEST(otbHorizontalSobelVectorImageFilterNewTest);
REGISTER_TEST(otbHorizontalSobelVectorImageFilterTest);
REGISTER_TEST(otbVerticalSobelVectorImageFilterNewTest);
REGISTER_TEST(otbVerticalSobelVectorImageFilterTest);
REGISTER_TEST(otbSobelVectorImageFilterNewTest);
REGISTER_TEST(otbSobelVectorImageFilterTest);
REGISTER_TEST(otbFastICAInternalOptimizerVectorImageFilterNewTest);
REGISTER_TEST(otbFastICAImageFilterNewTest);
REGISTER_TEST(otbFastICAImageFilterTest);
REGISTER_TEST(otbVectorImageToMatrixNewTest);
REGISTER_TEST(otbVectorImageToMatrixTest);
REGISTER_TEST(otbProjectiveProjectionTestHighSNR);
}
......@@ -20,8 +20,7 @@ ADD_SUBDIRECTORY(ObjectDetection)
ADD_SUBDIRECTORY(Fuzzy)
ADD_SUBDIRECTORY(MultiTemporal)
ADD_SUBDIRECTORY(Simulation)
ADD_SUBDIRECTORY(Hyperspectral)
IF(OTB_USE_VISU_GUI)
ADD_SUBDIRECTORY(Visu)
......
IF( NOT OTB_DISABLE_CXX_TESTING AND BUILD_TESTING )
SET(BASELINE ${OTB_DATA_ROOT}/Baseline/OTB/Images)
SET(BASELINE_FILES ${OTB_DATA_ROOT}/Baseline/OTB/Files)
SET(INPUTDATA ${OTB_DATA_ROOT}/Input)
SET(TEMP ${OTBTesting_BINARY_DIR}/Temporary)
IF(OTB_DATA_USE_LARGEINPUT)
SET(LARGEINPUT ${OTB_DATA_LARGEINPUT_ROOT} )
ENDIF(OTB_DATA_USE_LARGEINPUT)
SET(NOTOL 0.0)
SET(EPSILON_9 0.000000001)
SET(EPSILON_3 0.001)
SET(Hyperspectral_TESTS1 ${CXX_TEST_PATH}/otbHyperspectralTests1)
SET(Hyperspectral_TESTS2 ${CXX_TEST_PATH}/otbHyperspectralTests2)
# HyperspectralTests1
ADD_TEST(hyTuHyperspectralVariableNew ${Hyperspectral_TESTS1}
otbHyperspectralVariableNew)
SET(Hyperspectral_SRCS1
otbHyperspectralTests1.cxx
otbEigenvalueLikelihoodMaximization.cxx
otbCLSPSTOUnmixingImageFilter.cxx
otbFCLSUnmixingImageFilter.cxx
otbFullyConstrainedLeastSquareImageFilter.cxx
otbISRAUnmixingImageFilter.cxx
otbNCLSUnmixingImageFilter.cxx
#otbSparseUnmixingImageFilter.cxx
otbUnConstrainedLeastSquareImageFilter.cxx
otbVCAImageFilter.cxx
)
SET(Hyperspectral_SRCS2
otbHyperspectralTests2.cxx
otbLocalRxDetectorTest.cxx
otbLocalRxDetectorRoiTest.cxx
)
ADD_EXECUTABLE(otbHyperspectralTests1 otbHyperspectralTests1.cxx ${Hyperspectral_SRCS1})
TARGET_LINK_LIBRARIES(otbHyperspectralTests1 OTBIO OTBHyperspectral OTBTesting)
ADD_EXECUTABLE(otbHyperspectralTests2 otbHyperspectralTests2.cxx ${Hyperspectral_SRCS2})
TARGET_LINK_LIBRARIES(otbHyperspectralTests2 OTBIO OTBHyperspectral OTBTesting)
ENDIF( NOT OTB_DISABLE_CXX_TESTING AND BUILD_TESTING )
......@@ -25,11 +25,6 @@ void RegisterTests()
{
REGISTER_TEST(otbEigenvalueLikelihoodMaximizationNewTest);
REGISTER_TEST(otbEigenvalueLikelihoodMaximizationTest);
REGISTER_TEST(otbStreamingStatisticsVectorImageFilterTest);
// REGISTER_TEST(otbStreamingStatisticsVectorImageFilterNewTest);
// REGISTER_TEST(otbStreamingStatisticsVectorImageFilterTest);
REGISTER_TEST(otbMatrixMultiplyImageFilterNewTest);
REGISTER_TEST(otbMatrixMultiplyImageFilterTest);
REGISTER_TEST(otbUnConstrainedLeastSquareImageFilterNewTest);
REGISTER_TEST(otbUnConstrainedLeastSquareImageFilterTest);
REGISTER_TEST(otbISRAUnmixingImageFilterNewTest);
......@@ -42,34 +37,7 @@ void RegisterTests()
REGISTER_TEST(otbCLSPSTOUnmixingImageFilterTest);
REGISTER_TEST(otbFullyConstrainedLeastSquareImageFilterNewTest);
REGISTER_TEST(otbFullyConstrainedLeastSquareImageFilterTest);
REGISTER_TEST(otbVectorImageToMatrixNewTest);
REGISTER_TEST(otbVectorImageToMatrixTest);
REGISTER_TEST(vahineVCA);
REGISTER_TEST(vahineElmVCA);
REGISTER_TEST(otbProjectiveProjectionTestHighSNR);
REGISTER_TEST(otbVCATestHighSNR);
REGISTER_TEST(otbVCAImageFilterTestHighSNR);
REGISTER_TEST(otbHorizontalSobelVectorImageFilterNewTest);
REGISTER_TEST(otbHorizontalSobelVectorImageFilterTest);
REGISTER_TEST(otbVerticalSobelVectorImageFilterNewTest);
REGISTER_TEST(otbVerticalSobelVectorImageFilterTest);
REGISTER_TEST(otbSobelVectorImageFilterNewTest);
REGISTER_TEST(otbSobelVectorImageFilterTest);
REGISTER_TEST(otbLocalGradientVectorImageFilterNewTest);
REGISTER_TEST(otbLocalGradientVectorImageFilterTest);
REGISTER_TEST(otbLocalActivityVectorImageFilterNewTest);
REGISTER_TEST(otbLocalActivityVectorImageFilterTest);
REGISTER_TEST(otbNormalizeVectorImageFilterNewTest);
REGISTER_TEST(otbNormalizeVectorImageFilterTest);
REGISTER_TEST(otbPCAImageFilterNewTest);
REGISTER_TEST(otbPCAImageFilterTest);
REGISTER_TEST(otbMNFImageFilterNewTest);
REGISTER_TEST(otbMNFImageFilterTest);
REGISTER_TEST(otbNAPCAImageFilterNewTest);
REGISTER_TEST(otbNAPCAImageFilterTest);
REGISTER_TEST(otbFastICAInternalOptimizerVectorImageFilterNewTest);
REGISTER_TEST(otbFastICAImageFilterNewTest);
REGISTER_TEST(otbFastICAImageFilterTest);
// REGISTER_TEST(otbSparseWvltToAngleMapperListFilterNewTest);
// REGISTER_TEST(otbAngularProjectionBinaryImageFilterNewTest);
// REGISTER_TEST(otbAngularProjectionBinaryImageFilterTest);
......
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
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.
=========================================================================*/
#include "otbImage.h"
#include "otbVectorImage.h"
#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
#include "otbDotProductImageFilter.h"
#include "otbProjectiveProjectionImageFilter.h"
#include "otbMatrixMultiplyImageFilter.h"
#include "otbVectorImageToMatrixImageFilter.h"
#include "otbStreamingMinMaxImageFilter.h"
#include "otbStreamingStatisticsImageFilter.h"
#include "otbStreamingStatisticsVectorImageFilter.h"
#include "itkAbsImageFilter.h"
#include "itkMersenneTwisterRandomVariateGenerator.h"
#include "vnl/algo/vnl_svd.h"
const unsigned int Dimension = 2;
typedef double PixelType;
typedef double PrecisionType;
typedef otb::Image<PixelType, Dimension> ImageType;
typedef otb::VectorImage<PixelType, Dimension> VectorImageType;
typedef ImageType::IndexType IndexType;
typedef ImageType::SizeType SizeType;
typedef ImageType::RegionType RegionType;
typedef otb::ImageFileReader<VectorImageType> ReaderType;
typedef itk::AbsImageFilter<ImageType, ImageType> AbsImageFilterType;
typedef otb::ProjectiveProjectionImageFilter<VectorImageType,VectorImageType,PrecisionType> ProjectiveProjectionImageFilterType;
typedef otb::DotProductImageFilter<VectorImageType,ImageType> DotProductImageFilterType;
typedef otb::MatrixMultiplyImageFilter<VectorImageType,VectorImageType,PrecisionType> MatrixMultiplyImageFilterType;
typedef otb::VectorImageToMatrixImageFilter<VectorImageType> VectorImageToMatrixImageFilterType;
typedef otb::ImageFileWriter<VectorImageType> WriterType;
typedef otb::StreamingMinMaxImageFilter<ImageType> StreamingMinMaxImageFilterType;
typedef otb::StreamingStatisticsVectorImageFilter<VectorImageType> StreamingStatisticsVectorImageFilterType;
typedef otb::StreamingStatisticsImageFilter<ImageType> StreamingStatisticsImageFilterType;
typedef itk::Statistics::MersenneTwisterRandomVariateGenerator RandomVariateGeneratorType;
typedef StreamingStatisticsVectorImageFilterType::MatrixType MatrixType;
#define logmatrixsize(p) std::cout << #p << " rows : " << p.rows() << " cols : " << p.cols() << std::endl;
#define logvectorsize(p) std::cout << #p << " n : " << p.size() << std::endl;
#define logImagesize(p) p->UpdateOutputInformation(); std::cout << #p << " size : " << p->GetLargestPossibleRegion().GetSize() << " nbBands : " << p->GetNumberOfComponentsPerPixel() << std::endl;
int otbVCATestHighSNR(int argc, char * argv[])
{
const char * inputImage = argv[1];
const char * outputImage = argv[2];
const unsigned int nbEndmembers = atoi(argv[3]);
ReaderType::Pointer readerImage = ReaderType::New();
readerImage->SetFileName(inputImage);
readerImage->UpdateOutputInformation();
const unsigned int nbBands = readerImage->GetOutput()->GetNumberOfComponentsPerPixel();
std::cout << "Computing image stats" << std::endl;
StreamingStatisticsVectorImageFilterType::Pointer statsInput = \
StreamingStatisticsVectorImageFilterType::New();
statsInput->SetInput(readerImage->GetOutput());
statsInput->Update();
std::cout << "Computing SVD of correlation matrix" << std::endl;
// Take the correlation matrix
vnl_matrix<PrecisionType> R = statsInput->GetCorrelation().GetVnlMatrix();
// Apply SVD
vnl_svd<PrecisionType> svd(R);
vnl_matrix<PrecisionType> U = svd.U();
vnl_matrix<PrecisionType> Ud = U.get_n_columns(0, nbEndmembers);
vnl_matrix<PrecisionType> UdT = Ud.transpose();
logmatrixsize(U);
logmatrixsize(Ud);
std::cout << "Apply dimensionality reduction" << std::endl;
// Xd = Ud.'*M;
MatrixMultiplyImageFilterType::Pointer mulUd = MatrixMultiplyImageFilterType::New();
mulUd->SetInput(readerImage->GetOutput());
mulUd->SetMatrix(UdT);
mulUd->UpdateOutputInformation();
VectorImageType::Pointer Xd = mulUd->GetOutput();
logImagesize(Xd);
// Compute mean(Xd)
std::cout << "Compute mean(Xd)" << std::endl;
StreamingStatisticsVectorImageFilterType::Pointer statsXd = \
StreamingStatisticsVectorImageFilterType::New();
statsXd->SetInput(Xd);
statsXd->Update();
VectorImageType::PixelType Xdmean = statsXd->GetMean();
// Projective projection
// Xd ./ repmat( sum( Xd .* repmat(u,[1 N]) ) ,[d 1]);
std::cout << "Compute projective projection" << std::endl;
ProjectiveProjectionImageFilterType::Pointer proj = ProjectiveProjectionImageFilterType::New();
proj->SetInput(Xd);
proj->SetProjectionDirection(Xdmean);
VectorImageType::Pointer Y = proj->GetOutput();
logImagesize(Y);
// E : result, will contain the endmembers
vnl_matrix<PrecisionType> E(nbBands, nbEndmembers);
logmatrixsize(E);
// A = zeros(q,q)
// A(q,1) = 1
vnl_matrix<PrecisionType> A(nbEndmembers, nbEndmembers);
A.fill(0);
A(nbEndmembers - 1, 0) = 1;
RandomVariateGeneratorType::Pointer randomGen = RandomVariateGeneratorType::New();
for (unsigned int i = 0; i < nbEndmembers; ++i)
{
std::cout << "Iteration " << i << std::endl;
// w = rand(q,1)
std::cout << "Random vector generation " << std::endl;
vnl_vector<PrecisionType> w(nbEndmembers);
for (unsigned int j = 0; j < nbEndmembers; ++j)
{
w(j) = randomGen->GetVariateWithOpenRange();
}
// f = ((I - A*pinv(A))*w) / (norm(I - A*pinv(A))*w))
std::cout << "f = ((I - A*pinv(A))*w) /(norm(I - A*pinv(A))*w))" << std::endl;
vnl_matrix<PrecisionType> tmpMat(nbEndmembers,nbEndmembers);
tmpMat.set_identity();
std::cout << "A" << std::endl << A << std::endl;
vnl_svd<PrecisionType> Asvd(A);
tmpMat -= A * Asvd.inverse();
vnl_vector<PrecisionType> tmpNumerator = tmpMat * w;
vnl_vector<PrecisionType> f = tmpNumerator / tmpNumerator.two_norm();
// v = f.'*Y
std::cout << "v = f.'*Y" << std::endl;
DotProductImageFilterType::Pointer dotfY = DotProductImageFilterType::New();
dotfY->SetInput(Y);
VectorImageType::PixelType fV(f.data_block(), f.size());
dotfY->SetVector(VectorImageType::PixelType(fV));
ImageType::Pointer v = dotfY->GetOutput();
logImagesize(v);
// abs(v)
std::cout << "abs(v)" << std::endl;
AbsImageFilterType::Pointer absVFilter = AbsImageFilterType::New();
absVFilter->SetInput(v);
// max(abs(v))
std::cout << "max(abs(v))" << std::endl;
StreamingMinMaxImageFilterType::Pointer maxAbs = StreamingMinMaxImageFilterType::New();
maxAbs->SetInput(absVFilter->GetOutput());
maxAbs->Update();
// k = arg_max( max(abs(v)) )
std::cout << "k = arg_max( max(abs(v)) )" << std::endl;
IndexType maxIdx = maxAbs->GetMaximumIndex();
std::cout << "maxIdx : " << maxIdx << std::endl;
// extract Y(:,k)
std::cout << "Y(:,k)" << std::endl;
RegionType region;
region.SetIndex( maxIdx );
SizeType size;
size.Fill(1);
region.SetSize(size);
Y->SetRequestedRegion(region);
Y->Update();
// store new endmember in A
// A(:,i) = Y(:,k)
std::cout << "A(:,i) = Y(:,k)" << std::endl;
VectorImageType::PixelType e = Y->GetPixel(maxIdx);
A.set_column(i, e.GetDataPointer());
// reproject in original space
// u = Ud * Xd(:,k)
std::cout << "u = Ud * Xd(:,k)" << std::endl;
Xd->SetRequestedRegion(region);
Xd->Update();
VectorImageType::PixelType xd = Xd->GetPixel(maxIdx);
vnl_vector<PrecisionType> xdV(xd.GetDataPointer(), xd.GetSize());
logvectorsize(xdV);
logmatrixsize(Ud);
vnl_vector<PrecisionType> u = Ud * xdV;
logvectorsize(u);
// E(:, i) = u
std::cout << "E(:, i) = u" << std::endl;
E.set_column(i, u);
}
VectorImageType::Pointer output = VectorImageType::New();
RegionType region;
IndexType index;
SizeType size;
index.Fill(0);
size[0] = nbEndmembers;
size[1] = 1;
region.SetIndex(index);
region.SetSize(size);
output->SetRegions(region);
output->SetNumberOfComponentsPerPixel(nbBands);
output->Allocate();
itk::ImageRegionIteratorWithIndex<VectorImageType> it(output, output->GetLargestPossibleRegion());
unsigned int i;
for (it.GoToBegin(), i = 0; !it.IsAtEnd(); ++it, ++i)
{
vnl_vector<PrecisionType> e = E.get_column(i);
VectorImageType::PixelType pixel(e.data_block(), e.size());
it.Set(pixel);
}
WriterType::Pointer writer = WriterType::New();
writer->SetFileName(outputImage);
writer->SetInput(output);
writer->Update();
return EXIT_SUCCESS;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment