[Insight-users] Calculating least squares straight-line
C.Cagatay Bilgin
bilgic at cs.rpi.edu
Mon Nov 23 13:47:39 EST 2009
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
More information about the Insight-users
mailing list