[Insight-users] Fwd: Calculating least squares straight-line

[Nicholas Tustison]

Begin forwarded message:

> From: Nicholas Tustison <ntustison at gmail.com>
> Date: November 24, 2009 12:49:28 PM EST
> To: Matthias Dodt <matthias.dodt at mdc-berlin.de>
> Subject: Re: [Insight-users] Calculating least squares straight-line
> 
> Hi Matthias,
> 
> It seems a bit of an overkill to use the FEM functionality for least-squares fitting to a line.  One option is to use the B-spline approximation filter
> 
> http://www.itk.org/Doxygen/html/classBSplineScatteredDataPointSetToImageFilter_1_1h.html#_details
> 
> with the following parameter settings.
> 
> 1. Set the spline order to 1 (to specify a linear B-spline basis)
> 2. Set the number of levels to 1.
> 3. Set the number of control points to 2.
> 
> The resulting control point image will define the endpoints of your approximating line.  
> 
> Nick
> 
> 
> 
> On Nov 24, 2009, at 4:10 AM, Matthias Dodt 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
> 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091124/e86f14b2/attachment-0001.htm>