[ITK-users] ITK and VNL - SVD decomposition of image - Seg fault

[Tejas]

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.