[ITK-users] vnl_svd Segfaults on 3xn matrix, but not on nx3?

[DVigneault]

All--

I'm trying to fit a plane through an image (treating indices as x and y and
pixel value as z) by reading an image into a 3xn matrix and computing the
singular value decomposition with vnl.  However, I'm getting a segfault when
I call the SVD function.  I do not get a segfault when I call the SVD on the
transpose of the same matrix.  Could anyone suggest (a) why this is
segfaulting or (b) a better way to fit a plane to an image in ITK?

Below is a minimal example of the problem.  (Apologies for the many
questions--I very much appreciate Brad's and other's help.)

#include "itkImageFileReader.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkVariableSizeMatrix.h"
#include "vnl/vnl_matrix.h"

const unsigned int             Dimension = 2;
typedef double                 WorkPixelType;
typedef itk::Image< WorkPixelType, Dimension >    WorkImageType;
typedef itk::ImageFileReader< WorkImageType > ReaderType;
typedef itk::VariableSizeMatrix< WorkPixelType > MatrixType;
typedef itk::ImageRegionIteratorWithIndex< WorkImageType > ItType;

int main( int argc, char * argv[] )
{

	// Read in file
    ReaderType::Pointer reader = ReaderType::New();
	reader->SetFileName( "../data/input/wrapped_WT_square.png" );
    reader->Update();
    
    WorkImageType::Pointer input = reader->GetOutput();
    
    MatrixType matrix;
    matrix.SetSize(3,
input->GetLargestPossibleRegion().GetSize()[0]*input->GetLargestPossibleRegion().GetSize()[0]);
    
    ItType inputIt( input, input->GetLargestPossibleRegion() );
    inputIt.GoToBegin();
    unsigned int count = 0;
    while(!inputIt.IsAtEnd()) {
    	
    	
    	matrix(0,count) = inputIt.GetIndex()[0];
    	matrix(1,count) = inputIt.GetIndex()[1];
    	matrix(2,count) = inputIt.Get();
    
    	++count;
    	++inputIt;
    	
    }

    vnl_matrix< double > VNLMatrix = matrix.GetVnlMatrix();

    // n x 3 matrix works:
    vnl_svd<double> svdTranspose(VNLMatrix.transpose()); // Segmentation
fault: 11
    std::cout << svdTranspose << std::endl;
    
    // 3 x n matrix segfaults:
    vnl_svd<double> svd(VNLMatrix); // Segmentation fault: 11
    std::cout << svd << std::endl; // Doesn't make it this far

  return EXIT_SUCCESS;
  
}




--
View this message in context: http://itk-insight-users.2283740.n2.nabble.com/vnl-svd-Segfaults-on-3xn-matrix-but-not-on-nx3-tp7585974.html
Sent from the ITK Insight Users mailing list archive at Nabble.com.