<html><head></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><br><div><br><div>Begin forwarded message:</div><br class="Apple-interchange-newline"><blockquote type="cite"><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px;"><span style="font-family:'Helvetica'; font-size:medium; color:rgba(0, 0, 0, 1);"><b>From: </b></span><span style="font-family:'Helvetica'; font-size:medium;">Nicholas Tustison <<a href="mailto:ntustison@gmail.com">ntustison@gmail.com</a>><br></span></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px;"><span style="font-family:'Helvetica'; font-size:medium; color:rgba(0, 0, 0, 1);"><b>Date: </b></span><span style="font-family:'Helvetica'; font-size:medium;">November 24, 2009 12:49:28 PM EST<br></span></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px;"><span style="font-family:'Helvetica'; font-size:medium; color:rgba(0, 0, 0, 1);"><b>To: </b></span><span style="font-family:'Helvetica'; font-size:medium;">Matthias Dodt <<a href="mailto:matthias.dodt@mdc-berlin.de">matthias.dodt@mdc-berlin.de</a>><br></span></div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px;"><span style="font-family:'Helvetica'; font-size:medium; color:rgba(0, 0, 0, 1);"><b>Subject: </b></span><span style="font-family:'Helvetica'; font-size:medium;"><b>Re: [Insight-users] Calculating least squares straight-line</b><br></span></div><br><div>Hi Matthias,<br><br>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<br><br><a href="http://www.itk.org/Doxygen/html/classBSplineScatteredDataPointSetToImageFilter_1_1h.html#_details">http://www.itk.org/Doxygen/html/classBSplineScatteredDataPointSetToImageFilter_1_1h.html#_details</a><br><br>with the following parameter settings.<br><br>1. Set the spline order to 1 (to specify a linear B-spline basis)<br>2. Set the number of levels to 1.<br>3. Set the number of control points to 2.<br><br>The resulting control point image will define the endpoints of your approximating line. <br><br>Nick<br><br><br><br>On Nov 24, 2009, at 4:10 AM, Matthias Dodt wrote:<br><br><blockquote type="cite">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...?<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite">I would need a MultipleValuedLinearOptimizer- but i am not sure which class to use-<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite">best,<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite">mat<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite">C.Cagatay Bilgin schrieb:<br></blockquote><blockquote type="cite"><blockquote type="cite">I did not know this functionality existed in ITK and implemented<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">my own. Can you try running this on a smaller set of points and<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">provide those points to us so that I can give it a try as well.<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">Also, I see the following lines in the source code that looks related<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">to your problem. If you have zero elements, just setting them to<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">something really small like 1.0e-16 might help, maybe.<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">/* *******************************************************************<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00541 * FIX ME: itpack does not allow for any non-zero diagonal elements<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00542 * so "very small" numbers are inserted to allow for a solution<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00543 *<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00544 int i;<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00545 doublereal fakeZero = 1.0e-16;<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00546<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00547 //insert "fake" zeros<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00548 for (i=0; i<static_cast<int>(m_Order); i++)<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00549 {<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00550 if ( (*m_Matrices)[0].Get(i,i) == 0.0)<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00551 {<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00552 (*m_Matrices)[0].Set(i,i,fakeZero);<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00553 }<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00554 }<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00555 // END FIXME<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">00556 *********************************************************************/<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite">On Nov 23, 2009, at 10:22 AM, Matthias Dodt wrote:<br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">Hi!<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">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:<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">void CalculateLeastSquares(){<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> using namespace fem;<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> //Least squares: b + A * x = y<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> //given Pi(xi,yi): 1b + xi * x = yi<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> fem::LinearSystemWrapperItpack lsw;<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> lsw.SetNumberOfMatrices(1u);<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> lsw.SetNumberOfVectors(1u);<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> lsw.SetNumberOfSolutions(1u);<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> lsw.SetSystemOrder(areaPoints->size());<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> lsw.SetMaximumNonZeroValuesInMatrix(areaPoints->size() *4);<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> lsw.InitializeMatrix(0);<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> lsw.InitializeSolution(0);<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> lsw.InitializeVector(0);<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> //Matrix<double,2> coeff;<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> //Vector<double,1> y;<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> for(unsigned int t=0;t<areaPoints->size();++t){<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> //coeff(0,i) = 1;<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> std::cout << "Setting: " << t << "," << areaPoints->at(t)[0] << "," << areaPoints->at(t)[1] << std::endl;<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> lsw.SetMatrixValue(t,0u,1u,0);<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> lsw.SetMatrixValue(t,1,static_cast<float>(areaPoints->at(t)[0]),0);<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> //coeff(1,i) = static_cast<double>(areaPoints->at(i)[0]);<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> lsw.SetSolutionValue(t,static_cast<float>(areaPoints->at(t)[1]),0);<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> //y(i) = static_cast<double>(areaPoints->at(i)[1]);<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> }<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> lsw.Solve();<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> for(unsigned int i=0;i<lsw.GetNumberOfSolutions();++i){<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> std::cout << "Coefficient: " << lsw.GetSolutionValue(i,0);<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> }<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"> }<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">Anyway the .Solve() crashes with:<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">itk::FEMExceptionItpackSolver (0x66c1a0)<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">Location: "LinearSystemWrapperItpack::Solve"<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">File: /Users/mat/Documents/InsightToolkit-3.14.0/Code/Numerics/FEM/itkFEMLinearSystemWrapperItpack.cxx <br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">Line: 659<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">Description: Error: A diagonal element is not positive<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">Somebody has a clue?<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">thanks!<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">_____________________________________<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">Powered by www.kitware.com<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">Visit other Kitware open-source projects at<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">http://www.kitware.com/opensource/opensource.html<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">Kitware offers ITK Training Courses, for more information visit:<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">http://www.kitware.com/products/protraining.html<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">Please keep messages on-topic and check the ITK FAQ at:<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">http://www.itk.org/Wiki/ITK_FAQ<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">Follow this link to subscribe/unsubscribe:<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><blockquote type="cite">http://www.itk.org/mailman/listinfo/insight-users<br></blockquote></blockquote></blockquote><blockquote type="cite"><blockquote type="cite"><br></blockquote></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite">-- <br></blockquote><blockquote type="cite">------------------------------------------------<br></blockquote><blockquote type="cite">Matthias Dodt<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite">Scientific Programmer at Bioinformaitcs platform AG Dieterich<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite">Berlin Institute for Medical Systems Biology at the<br></blockquote><blockquote type="cite">Max-Delbrueck-Center for Molecular Medicine<br></blockquote><blockquote type="cite">Robert-Roessle-Strasse 10, 13125 Berlin, Germany<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite">fon: +49 30 9406 4261<br></blockquote><blockquote type="cite">email: matthias.dodt@mdc-berlin.de<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite">_____________________________________<br></blockquote><blockquote type="cite">Powered by www.kitware.com<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite">Visit other Kitware open-source projects at<br></blockquote><blockquote type="cite">http://www.kitware.com/opensource/opensource.html<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite">Kitware offers ITK Training Courses, for more information visit:<br></blockquote><blockquote type="cite">http://www.kitware.com/products/protraining.html<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite">Please keep messages on-topic and check the ITK FAQ at:<br></blockquote><blockquote type="cite">http://www.itk.org/Wiki/ITK_FAQ<br></blockquote><blockquote type="cite"><br></blockquote><blockquote type="cite">Follow this link to subscribe/unsubscribe:<br></blockquote><blockquote type="cite">http://www.itk.org/mailman/listinfo/insight-users<br></blockquote><br></div></blockquote></div><br></body></html>