[Insight-users] Re: Using KalmanLinearEstimator: patch

[Mathieu Malaterre]

Hi Luis,

  Indeed giving a better starting point make the convergence go a
little faster (it looks like this a pretty bad convergence still).

  Anyway the patch is totally backward compatible.

HTH
-Mathieu

On 5/7/07, Mathieu Malaterre <mathieu.malaterre at gmail.com> wrote:
> As usual Luis thanks a bunch for your help. I will look into your suggestion.
>
> just for reference I was able to get a correct approximate solution by doing:
>
>   for(unsigned int duplicate = 0; duplicate < 100000; ++duplicate)
>     for(unsigned int ax=0; ax < N; ax++)
>       {
>       predictor(0) = points[ax][0];
>       measure = points[ax][1];
>       //std::cout << predictor(0) << "," << measure << std::endl;
>       filter.UpdateWithNewMeasure(measure,predictor);
>       }
>
> If you plot the values using excel or gnuplot you can see the linear
> regression is pretty obvious. Anyway I do not know enough about Kalman
> algorithm.
>
> Thanks
> -Mathieu
>
> On 5/7/07, Luis Ibanez <luis.ibanez at kitware.com> wrote:
> >
> > Hi Mathieu,
> >
> >
> > 1) Do you have only 6 data points in your real problem ?
> >
> >     They may not be enough data for bringing the
> >     Kalman Estimator to convergence to a good
> >     estimation.
> >
> >     I was assuming that you had in the order
> >     of hundreds to thousands of values.
> >
> >
> > 2) You may find useful to print out the estimations
> >     as you go adding new points. One of the interesting
> >     advantages of the Kalman estimator is that it
> >     provides at estimation using the currently provided
> >     samples.
> >
> >     E.g., move the line:
> >
> >  >  VectorType estimation = filter.GetEstimator();
> >
> >     inside the for loop and print out the value of
> >     the estimation. Then you will see how it get
> >     close to the [A,B] values.
> >
> >     Note that the estimator is initialized with [0,0]
> >     so it may take a while (many measurement) to shift
> >     the estimation to your real A,B values.
> >
> >     Maybe we should add a method to initialize the
> >     Estimator to a known value. In this way, the
> >     Kalman estimator would focus on refining the
> >     estimation.
> >
> >
> >
> >     Regards,
> >
> >
> >       Luis
> >
> >
> > -------------------------
> > Mathieu Malaterre wrote:
> > > Hello,
> > >
> > >  I must be something obviously really wrong when using the
> > > KalmanLinearEstimator class, because the results are way off. Could
> > > someone please help find the obvious mistake I must have done.
> > >
> > > thanks !
> > > -Mathieu
> > >
> > > Code is:
> > >
> > > #include "itkKalmanLinearEstimator.h"
> > > #include <iostream>
> > >
> > > int main(int, char* [] )
> > > {
> > >  typedef itk::KalmanLinearEstimator<double,2> KalmanFilterType;
> > >
> > >  typedef KalmanFilterType::VectorType    VectorType;
> > >  typedef KalmanFilterType::MatrixType    MatrixType;
> > >  typedef KalmanFilterType::ValueType     ValueType;
> > >
> > >  KalmanFilterType filter;
> > >
> > >  filter.ClearEstimation();
> > >  filter.SetVariance(1.0);
> > >
> > >  ValueType     measure;
> > >  VectorType    predictor;
> > >
> > >  // A ~ 0.323
> > >  // B ~ 172.4
> > >  const double points[][2] = {
> > >    {-16.6817,167.095},
> > >    {-44.4257,158.118},
> > >    {-68.0884,150.462},
> > >    {-87.7176,144.111},
> > >    {-103.447,139.021},
> > >    {-115.437,135.142},
> > >  };
> > >  const unsigned int N = sizeof(points) / ( 2 * sizeof( double ));
> > >
> > >  predictor(1)  =  1.0;
> > >  for(unsigned int ax=0; ax < N; ax++)
> > >    {
> > >    predictor(0) = points[ax][0];
> > >    measure = points[ax][1];
> > >    //std::cout << predictor(0) << "," << measure << std::endl;
> > >    filter.UpdateWithNewMeasure(measure,predictor);
> > >    }
> > >
> > >  VectorType estimation = filter.GetEstimator();
> > >
> > >  std::cout << std::endl << "The Estimation is : " << std::endl;
> > >  std::cout << estimation;
> > >
> > >  std::cout << std::endl << std::endl;
> > >  return EXIT_SUCCESS;
> > > }
> > >
> > >
> > > On 5/6/07, Luis Ibanez <luis.ibanez at kitware.com> wrote:
> > >
> > >>
> > >> Hi Mathieu,
> > >>
> > >> You may want to try the KalmanLinearEstimator.
> > >>
> > >> http://www.itk.org/Insight/Doxygen/html/classitk_1_1KalmanLinearEstimator.html
> > >>
> > >>
> > >>
> > >> Your estimator will be the vector [a,b],
> > >>
> > >> and your measures will be the vector [1,Y]
> > >>
> > >>
> > >> The Kalman Linear Estimator should converge to the vector [a,b] that
> > >> minimizes the sum of squared errors between X and (a+bY).
> > >>
> > >> For an example, you may want to look at:
> > >>
> > >>
> > >>
> > >>      Insight/Testing/Code/Algorithms/
> > >>           itkKalmanLinearEstimatorTest.cxx
> > >>
> > >>
> > >>
> > >>     Regards,
> > >>
> > >>
> > >>         Luis
> > >>
> > >>
> > >> ---------------------------
> > >> Mathieu Malaterre wrote:
> > >> > Hello,
> > >> >
> > >> >  I have a pretty simple problem to solve and I was wondering if I
> > >> > could reuse any brick already existing in ITK.
> > >> >
> > >> >  I am working with a 4D dataset (x,y,y and theta), where I can
> > >> > express for each pixel that:
> > >> >
> > >> > I(x,y,z)/sin(theta) = F(x,y,z) + b*I(x,y,z)/tan(theta)
> > >> >
> > >> >  My goal here is to find out the b constant, I do not need the
> > >> > F(x,y,z) part. Obvisouly all I need to do is draw the line
> > >> >
> > >> >   X = a + b Y
> > >> >
> > >> > where X = I(x,y,z)/sin(theta) and Y = I(x,y,z)/tan(theta) and the
> > >> > slope will give me 'b'.
> > >> >
> > >> > Any suggestions on which class to reuse ?
> > >> > thanks !
> > >>
> > >
> > >
> >
>
>
> --
> Mathieu
> Tel: +33 6 32 13 33 40
>


-- 
Mathieu
Tel: +33 6 32 13 33 40
-------------- next part --------------
A non-text attachment was scrubbed...
Name: kalman.patch
Type: application/octet-stream
Size: 894 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20070507/7b3745e0/kalman.obj