[Insight-users] Fwd: Calculating least squares straight-line
Karthik Krishnan
karthik.krishnan at kitware.com
Tue Nov 24 13:55:46 EST 2009
Yes..
- total overkill to use FEM...
- The BSplineScatteredApprox filter might be an overkill as well...
The orientation of LS fit direction to a point cloud is simply the
eigen vector of the covariance matrix of the points. The origin of
that line would be the centroid of the points.. something like the
following excerpt, from code I've written as part of a larger app....
void AddPoint( const Point p )
{
++m_Size;
// Update the centroid...
m_Sum += p;
const double n = static_cast< double >(m_Size);
for (unsigned int i = 0; i < PointType::PointDimension; i++)
{
m_Centroid[i] = m_Sum[i] / n;
}
// Update the convariance...
m_CovarianceSum[5] += (p[0] * p[0]);
m_CovarianceSum[4] += (p[0] * p[1]);
m_CovarianceSum[3] += (p[1] * p[1]);
m_CovarianceSum[2] += (p[0] * p[2]);
m_CovarianceSum[1] += (p[1] * p[2]);
m_CovarianceSum[0] += (p[2] * p[2]);
m_Covariance[5] = m_CovarianceSum[5] - ( n * m_Centroid[0] *
m_Centroid[0] );
m_Covariance[4] = m_CovarianceSum[4] - ( n * m_Centroid[1] *
m_Centroid[0] );
m_Covariance[3] = m_CovarianceSum[3] - ( n * m_Centroid[1] *
m_Centroid[1] );
m_Covariance[2] = m_CovarianceSum[2] - ( n * m_Centroid[2] *
m_Centroid[0] );
m_Covariance[1] = m_CovarianceSum[1] - ( n * m_Centroid[2] *
m_Centroid[1] );
m_Covariance[0] = m_CovarianceSum[0] - ( n * m_Centroid[2] *
m_Centroid[2] );
// Now fit a LS line to these points..
// Compute the line that best approximates the existing points in a least
// squares sense.
this->FitLeastSquaresLineToPoints();
}
void FitLeastSquaresLineToNeedlePoints()
{
CovarianceMatrixType::EigenValuesArrayType eigenValues;
CovarianceMatrixType::EigenVectorsMatrixType eigenVectors;
// EigenValues are sorted in ascending order by value. (sign is taken into
// account). EigenVectors are sorted in accordance with their eigen values.
m_Covariance.ComputeEigenAnalysis( eigenValues, eigenVectors );
this->m_LineEndPoints[0] = m_Centroid;
if ( eigenValues[2] != eigenValues[1] )
{
// regular case
PointType::VectorType direction;
direction[0] = eigenVectors[2][2];
direction[1] = eigenVectors[2][1];
direction[2] = eigenVectors[2][0];
// Point 1 is the centroid.
// Point 2 is the centroid + direction.
this->m_LineEndPoints[1] = m_Centroid + direction;
}
else
{
// isotropic case (infinite number of directions)
// by default: assemble a horizontal plane that goes through the centroid.
// cannot find LS line...
}
}
On Tue, Nov 24, 2009 at 12:53 PM, Nicholas Tustison <ntustison at gmail.com> wrote:
>
>
> 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
>
>
>
> _____________________________________
> 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
>
>
--
_________________________________
Karthik Krishnan
R&D Engineer,
Kitware Inc.
Ph: +1 5188814919, +91 9538477060
More information about the Insight-users
mailing list