[Insight-users] randomness in MRIRegistration example

[Luis Ibanez]

Hi Peter,



Yes, your observation is correct.

The MRIRegistration example uses a MutualInformationMetric.
This metric takes a population of pixels from the image
in order to evaluate MutulaInformation. Pixel sampling
is performed by the ImageRandomIteratorWithIndex class.

This iterator takes a random walk on N-Dimensions over a
selected region of the N-D image. The user define how many
samples should be extracted from the image region.

The randomness of the iterator walk is generated by the
"vnl_sample_uniform( a, b )" function which produces samples
with a uniform distribution in the interval [a,b].

The actual implementation of vnl_sample_uniform can be
found in

     Insight/Code/Numerics/vxl/vnl/vnl_sample.cxx

It is basically:


   double vnl_sample_uniform(double a, double b)
   {
   #if VXL_STDLIB_HAS_DRAND48
     // it's your lucky day.
     double u = drand48(); // uniform on [0, 1)
   #else
     // unlucky! perhaps your luck will change if you try
     // a different random number generator.
     vnl_sample_seed = (vnl_sample_seed*16807)%2147483647L;
     double u = double(vnl_sample_seed)/2147483711L;
   #endif
     return (1.0 - u)*a + u*b;
   }



As you can see, if you are lucky (in vnl terms ) the
drand48() function will generate the random number for you.
Otherwise an inlined method is used.

In order to prevent the sequence of random numbers to be
repeated at each run of the program, you can reinitialize
the seed using the function:

           vnl_sample_reseed();

The implementation is found in the same vnl_sample.cxx file.
It is basically:


   void vnl_sample_reseed()
   {
   # if VXL_STDLIB_HAS_DRAND48
     srand48( vcl_time(0) );
   # else
     vnl_sample_seed = vcl_time(0);
   # endif
   }


This is exactly what as you suggested: To reinitialize the
seed using the current time.


To summarize,
In order prevent the runs from being deterministics
you can call :

             vnl_sample_reseed();

before invoking:

            ->StartRegistration();

in the ImageRegistrationMethod.


It could make sense to invoke this method in the
RandomIterator on each call to GoToBegin(). Or
even in the constructor.  It could also be invoked
in the MutualInformationMetric before every metric
computation... but as you pointed out, this will
prevent users from repeating experiments.



One possibility for making everybody happy is to
set a boolean in the MutualInformation metric
that will enable/disable the reinitialization of
the seed.



Please let us know which of these options you
may find appropriate.




   Thanks


     Luis



====================================================
 >
 > I've been playing with the MRIRegistration example
 > code, and I notice that the random number seeding
 > appears to be constant (i.e. fixed starting seed
 > rather than time-based).  Multiple runs on the same
 > sample give identical results.  This is good for
 > testing purposes but seems not so good for obtaining
 > robust results, particularly given that the example
 > code involves running several iterations of the
 > starting program.
 > Is this observation indeed correct?
 > Thanks,
 > --Peter