[Insight-users] Calculating least squares straight-line

[Matthias Dodt]

Hi!

I have some problems using the LinearSystemWrapperItpack class. I want 
to compute the straight-line wich fits a set of points best. Till now i got:

void CalculateLeastSquares(){
        using namespace fem;
        //Least squares: b + A * x = y
        //given Pi(xi,yi): 1b + xi * x = yi

        fem::LinearSystemWrapperItpack lsw;
        lsw.SetNumberOfMatrices(1u);
        lsw.SetNumberOfVectors(1u);
        lsw.SetNumberOfSolutions(1u);
        lsw.SetSystemOrder(areaPoints->size());
        lsw.SetMaximumNonZeroValuesInMatrix(areaPoints->size() *4);
        lsw.InitializeMatrix(0);
        lsw.InitializeSolution(0);
        lsw.InitializeVector(0);


        //Matrix<double,2> coeff;
        //Vector<double,1> y;
        for(unsigned int t=0;t<areaPoints->size();++t){
            //coeff(0,i) = 1;
            std::cout << "Setting: " << t << "," << areaPoints->at(t)[0] 
<< "," << areaPoints->at(t)[1] << std::endl;
            lsw.SetMatrixValue(t,0u,1u,0);
            
lsw.SetMatrixValue(t,1,static_cast<float>(areaPoints->at(t)[0]),0);
            //coeff(1,i) = static_cast<double>(areaPoints->at(i)[0]);
            
lsw.SetSolutionValue(t,static_cast<float>(areaPoints->at(t)[1]),0);
            //y(i) = static_cast<double>(areaPoints->at(i)[1]);
        }

        lsw.Solve();

        for(unsigned int i=0;i<lsw.GetNumberOfSolutions();++i){
            std::cout << "Coefficient: " << lsw.GetSolutionValue(i,0);
        }



    }

Anyway the .Solve() crashes with:

itk::FEMExceptionItpackSolver (0x66c1a0)
Location: "LinearSystemWrapperItpack::Solve"
File: 
/Users/mat/Documents/InsightToolkit-3.14.0/Code/Numerics/FEM/itkFEMLinearSystemWrapperItpack.cxx
Line: 659
Description: Error: A diagonal element is not positive

Somebody has a clue?

thanks!