[Insight-users] Calculating least squares straight-line
Matthias Dodt
matthias.dodt at mdc-berlin.de
Thu Dec 3 08:41:54 EST 2009
Hi Sajendra!
Thanks for the advice- i have binary images containing "tubes" and just
want to compute the line for determining the orientation of those. So
its not based on intensities-
Thanks for the links!
Best,
mat
Sajendra schrieb:
> Hi Matthias,
>
> I'm not familiar with the itk FEM classes nor using them for data
> fitting; however, If you just need simple linear least solution to an
> overdetermined system can you not use the closed form solution?:
> http://en.wikipedia.org/wiki/Linear_least_squares#Computation
> http://mathworld.wolfram.com/LeastSquaresFitting.html
>
> It is quite simple to compute. I recently hacked an ITK image filter
> to perform that computation on input ITK image intensities so I could
> use ITKs multithreading and IO. Let me know if that would help with
> your problem.
>
> Sajendra
>
>
>
> On Tue, Nov 24, 2009 at 4:10 AM, Matthias Dodt
> <matthias.dodt at mdc-berlin.de> wrote:
>
>> Thanks for the hint - but i think the problem is that the matrix A has to be
>> a square matrix. My system of equations is overdetermined, so maybe i should
>> use a different solver...?
>>
>> I would need a MultipleValuedLinearOptimizer- but i am not sure which class
>> to use-
>>
>> best,
>>
>> mat
>>
>> C.Cagatay Bilgin schrieb:
>>
>>> 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
>>>>
>> --
>> ------------------------------------------------
>> Matthias Dodt
>>
>> Scientific Programmer at Bioinformaitcs platform AG Dieterich
>>
>> Berlin Institute for Medical Systems Biology at the
>> Max-Delbrueck-Center for Molecular Medicine
>> Robert-Roessle-Strasse 10, 13125 Berlin, Germany
>>
>> fon: +49 30 9406 4261
>> email: matthias.dodt at mdc-berlin.de
>>
>> _____________________________________
>> 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
>>
>>
--
------------------------------------------------
Matthias Dodt
Scientific Programmer at Bioinformatics platform AG Dieterich
Berlin Institute for Medical Systems Biology at the
Max-Delbrueck-Center for Molecular Medicine
Robert-Roessle-Strasse 10, 13125 Berlin, Germany
fon: +49 30 9406 4261
email: matthias.dodt at mdc-berlin.de
More information about the Insight-users
mailing list