[Insight-users] Calculating least squares straight-line

[C.Cagatay Bilgin]

I did not know this functionality existed in ITK and implemented
my own. Can you try running this on a smaller set of points and
provide those points to us so that I can give it a try as well.

Also, I see the following lines in the source code that looks related
to your problem. If you have zero elements, just setting them to
something really small like 1.0e-16 might help, maybe.

/* *******************************************************************
00541    * FIX ME: itpack does not allow for any non-zero diagonal  
elements
00542    * so "very small" numbers are inserted to allow for a solution
00543    *
00544   int i;
00545   doublereal fakeZero = 1.0e-16;
00546
00547   //insert "fake" zeros
00548   for (i=0; i<static_cast<int>(m_Order); i++)
00549   {
00550     if ( (*m_Matrices)[0].Get(i,i) == 0.0)
00551     {
00552       (*m_Matrices)[0].Set(i,i,fakeZero);
00553     }
00554   }
00555   // END FIXME
00556     
*********************************************************************/


On Nov 23, 2009, at 10:22 AM, Matthias Dodt wrote:

> 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!
>
> _____________________________________
> 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.html
>
> 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://www.itk.org/mailman/listinfo/insight-users