[ITK] [ITK-users] ITK and VNL - SVD decomposition of image - Seg fault
Tejas
tejas.rox at gmail.com
Sat Nov 22 20:04:31 EST 2014
Hi,
I am trying to decompose my affine-registered input ITK image using the SVD
routine in VNL. I have followed prior source code available from the
following links:
http://www.itk.org/pipermail/insight-users/2014-July/050694.html
http://public.kitware.com/vxl/doc/development/books/core/book_6.html
http://public.kitware.com/pipermail/insight-users/2007-February/020883.html
A part of my code, which is the root cause of the two problems below, is
pasted after a short explanation of what I am facing:
There are two problems that occur:
Problem 1:
While using the standard VNL SVD routine (vnl_svd) to decompose my input
image into the singular matrix (W) and two orthogonal matrices (U and V), I
witness segmentation faults when running my code from the command prompt.
This is similar to what one user had posted in the first link above.
Elaborating further:
My code seg faults when I use this line: vnl_svd<double> svd(VNLMatrix);
Yet, it runs fine when I use this line instead: vnl_svd<double>
svdTranspose(VNLMatrix.transpose());
Is there an inherent issue with the vnl_svd routine that has not been fixed
yet? If so, what are the other possible ways to run svd on a single 200x200
image?
Problem 2:
Recomposing the U,W,and V component matrices yields an abnormal result. The
recomposition of the matrices is returned as a vnl_matrix of a specific data
type; in my case, a double valued matrix. When visualizing the resulting
.MHA file in ITK-SNAP, the image appears all white. At first, I thought I
was seeing a scaled intensity of the image. So, I divided the image by it's
maximum value, and the result was pretty much the same, so scaling wasn't
the issue here. I believe that problem 2 is connected to problem 1, and I am
not sure about how to solve problem 1 itself. Does anybody have thoughts on
why this is the case?
Code:
typedef itk::VariableSizeMatrix< double > MatrixType;
// create a matrix structure in ITK
MatrixType matrix;
matrix.SetSize(rows, cols); // set the size of the ITK matrix
// create an image iterator
// dAlignedIm is a double data-type, ITK affine-registered image
passed to vnl_svd, via an ITK VariableSizeMatrix and vnl_matrix.
IterType inputIter( dAlignedIm, dAlignedIm->GetLargestPossibleRegion() );
inputIter.GoToBegin();
while(!inputIter.IsAtEnd()) {
// copy the input image intensity value into the corresponding matrix
unsigned int c = inputIter.GetIndex()[0];
unsigned int r = inputIter.GetIndex()[1];
matrix(r, c) = inputIter.Get();
++inputIter;
}
// create a vnl matrix
vnl_matrix< double > VNLMatrix = matrix.GetVnlMatrix();
//// perform SVD decomposition
vnl_svd<double> svdTranspose(VNLMatrix.transpose());
//svdTranspose.recompose();
//svdTranspose.W(0,0) = 0;
typedef itk::ImportImageFilter< double, Dimension > ImportFilterType;
ImportFilterType::Pointer importFilter = ImportFilterType::New();
ImportFilterType::SizeType size;
size[0] = cols; // size along X
size[1] = rows; // size along Y
ImportFilterType::IndexType start;
start.Fill( 0 );
ImportFilterType::RegionType region;
region.SetIndex( start );
region.SetSize( size );
importFilter->SetRegion( region );
double origin[ Dimension ];
origin[0] = 0.0; // X coordinate
origin[1] = 0.0; // Y coordinate
importFilter->SetOrigin( origin );
double spacing[ Dimension ];
spacing[0] = 1.0; // along X direction
spacing[1] = 1.0; // along Y direction
importFilter->SetSpacing( spacing );
const bool importImageFilterWillOwnTheBuffer = false;
vnl_matrix <double> reconstruct = (svdTranspose.recompose());
//reconstruct = reconstruct.transpose();
importFilter->SetImportPointer( reconstruct.data_block(), rows*cols,
importImageFilterWillOwnTheBuffer );
typedef itk::ImageFileWriter< doubleImageType > WriterType2;
//typedef itk::CastImageFilter< doubleImageType, ImageType >
CastFilterType3;
WriterType2::Pointer writernew = WriterType2::New();
//CastFilterType3::Pointer casternew = CastFilterType3::New();
writernew->SetFileName( "reconstruct.mha" );
//casternew->SetInput( importFilter->GetOutput() );
writernew->SetInput( importFilter->GetOutput() );
try
{
writernew->Update();
}
catch( itk::ExceptionObject & exp )
{
std::cerr << "Exception caught !" << std::endl;
std::cerr << exp << std::endl;
}
--
View this message in context: http://itk-insight-users.2283740.n2.nabble.com/ITK-and-VNL-SVD-decomposition-of-image-Seg-fault-tp7586604.html
Sent from the ITK Insight Users mailing list archive at Nabble.com.
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php
Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ
Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/insight-users
More information about the Community
mailing list