<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=iso-8859-1">
<META content="MSHTML 6.00.2800.1491" name=GENERATOR></HEAD>
<BODY>
<DIV><SPAN class=271203213-17032005><FONT face=Verdana color=#0000ff size=2>Ahh,
now I can sleep. Your analysis seems solid (although binomials were never
my thing).</FONT></SPAN></DIV>
<DIV><SPAN class=271203213-17032005><FONT face=Verdana color=#0000ff
size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=271203213-17032005><FONT face=Verdana color=#0000ff
size=2>Jim</FONT></SPAN></DIV>
<BLOCKQUOTE>
<DIV class=OutlookMessageHeader dir=ltr align=left><FONT face=Tahoma
size=2>-----Original Message-----<BR><B>From:</B>
insight-users-bounces@itk.org [mailto:insight-users-bounces@itk.org]<B>On
Behalf Of </B>Stefan Klein<BR><B>Sent:</B> Thursday, March 17, 2005 8:27
AM<BR><B>To:</B> Insight-users@itk.org<BR><B>Subject:</B> RE: [Insight-users]
random number generator<BR><BR></FONT></DIV><FONT size=3>Hi Jim,<BR><BR>You
observed well. I think this is the expected behaviour. I do not think the
truncation issues are the cause.<BR><BR>M = the number of pixels in the
image<BR>N = the number of loops over the image (200 in the example)<BR>R =
The number of pixels selected by the RandomIterator (= 100M in this
example)<BR><BR>image2:<BR>after one loop through the
image:<BR><X-TAB> </X-TAB>P(
pixel_i = 0x sampled ) = 1/2 for i =
1...M<BR><X-TAB> </X-TAB>P(
pixel_i = 1x sampled ) = 1/2<X-TAB> </X-TAB>for i = 1...M<BR>after N
loops through the image:
<BR><X-TAB> </X-TAB>P (pixel_i
= K times sampled) = \sum_{k=1...N} binom_coeff(N,K) (1/2)^K
(1/2)^{N-K}<X-TAB> </X-TAB>for i =
1...M<BR><X-TAB> </X-TAB>(which
is of course the binomial
distribution)<BR><X-TAB> </X-TAB>Expectation:
E( K ) = N/2 = 200/2 =
100<BR><X-TAB> </X-TAB>Variance:
Var( K ) = N/4 = 200/4 =
50<BR><X-TAB> </X-TAB>Standard
Deviation: SD( K ) = sqrt (N/4) = 7. as you observed.<BR><BR>image 4:<BR>after
one jump of the
RandomIterator<BR><X-TAB> </X-TAB>P(
pixel_i = 0x sampled ) = (M-1)/M for i =
1...M<BR><X-TAB> </X-TAB>P(
pixel_i = 1x sampled ) =
1/M<X-TAB> </X-TAB> for i =
1...M<BR>after R jumps of the
RandomIterator:<BR><X-TAB> </X-TAB>P(
pixel_i = K times sampled ) = \sum_{k=1...R} binom_coeff(R,K) (1/M)^K
((M-1)/M)^{R-K} for i =
1...M<BR><X-TAB> </X-TAB>(so
again a binomial
distribution)<BR><X-TAB> </X-TAB>Expectation:
E(K) = R/M = 100M / M =
100<BR><X-TAB> </X-TAB>Variance:
Var(K) = R/M * (M-1)/M = approximately R/M (if M is big) =
100.<BR><X-TAB> </X-TAB>Standard
Deviation: SD(K) = sqrt( var(K) ) = 10. as you observed.<BR><BR>My statistic
math classes were some time ago, so correct me if i made a mistake
:)<BR><BR>Regards,<BR>Stefan.<BR><BR>At 16:43 16/03/05, Miller, James V
(Research) wrote:<BR></FONT>
<BLOCKQUOTE class=cite cite="" type="cite"><FONT face=verdana color=#0000ff
size=2>An interesting observation:</FONT><FONT
size=3><BR> <BR></FONT><FONT face=verdana color=#0000ff size=2>The
standard deviation of the pixels in image2 and image4 are suprisingly
different. image2 has a standard deviation of 7 and image4 has a
standard deviation of 10. </FONT><FONT
size=3><BR> <BR></FONT><FONT face=verdana color=#0000ff size=2>If I
understand the original post correctly, both images were created using the
same random number generator. image2 was created by drawing a random
number 200 times at each pixel and counting the number of times the random
number was greater than 0.5. image4 was created using a
ImageRandomIteratorWithIndex to visit pixels in the image at random,
performing enough random visits so that each pixel was on average visited
100 times (the same expected value as the pixels in image2).</FONT><FONT
size=3><BR> <BR></FONT><FONT face=verdana color=#0000ff size=2>Naively,
I would expect the mean and standard deviations of these two images to be
very close.</FONT><FONT size=3><BR> <BR></FONT><FONT face=verdana
color=#0000ff size=2>The RandomIterator pretty much just draws a random
number between 0 and the number of pixels and unrolls the resulting linear
index into an index in the image. The truncation of the random number
to an integer and then to an index increases the standard deviation of how
many times a pixel is visited from 7 to 10.</FONT><FONT
size=3><BR> <BR></FONT><FONT face=verdana color=#0000ff size=2>Probably
not a useful observation.... but interesting.</FONT><FONT
size=3><BR> <BR></FONT><FONT face=verdana color=#0000ff
size=2>Jim</FONT><FONT size=3><BR></FONT>
<DL>
<DD><FONT face=tahoma size=2>-----Original Message-----<BR>
<DD>From:</B> insight-developers-bounces@itk.org [<A
href="mailto:insight-developers-bounces@itk.org"
eudora="autourl">mailto:insight-developers-bounces@itk.org</A>]On Behalf
Of </B>Blezek, Daniel J (Research)<BR>
<DD>Sent:</B> Wednesday, March 16, 2005 10:05 AM<BR>
<DD>To:</B> Stefan Klein; Insight-users@itk.org;
insight-developers@itk.org<BR>
<DD>Subject:</B> [Insight-developers] RE: [Insight-users] random number
generator<BR><BR></FONT>
<DD><FONT face=arial color=#0000ff size=2>Stefan, I started a thread
on this a while ago relating to non-deterministic behavior across
platforms. If the ITK design committee approves the move, it would
be great to standardize the random number generator to be sufficiently
random, fast and generate the same sequences (from a given seed) across
platforms. The discussion culminated with Brad's post: <A
href="http://www.itk.org/mailman/private/insight-developers/2005-January/006220.html">http://www.itk.org/mailman/private/insight-developers/2005-January/006220.html</A>
I'm not sure any action happened on the suggestions.</FONT><FONT
size=3><BR>
<DD><BR></FONT>
<DD><FONT face=arial color=#0000ff size=2>As you point out, many
registration algorithms depend on random sampling, exactly where I came
across the problem.</FONT><FONT size=3><BR>
<DD><BR></FONT>
<DD><FONT face=arial color=#0000ff size=2>Cheers,</FONT><FONT
size=3><BR></FONT>
<DD><FONT face=arial color=#0000ff size=2>-dan</FONT><FONT
size=3><BR></FONT>
<DD><FONT face=tahoma size=2>-----Original Message-----<BR>
<DD>From:</B> insight-users-bounces@itk.org [<A
href="mailto:insight-users-bounces@itk.org"
eudora="autourl">mailto:insight-users-bounces@itk.org</A>]On Behalf Of
</B>Stefan Klein<BR>
<DD>Sent:</B> Wednesday, March 16, 2005 9:16 AM<BR>
<DD>To:</B> Insight-users@itk.org<BR>
<DD>Subject:</B> [Insight-users] random number
generator<BR><BR></FONT><FONT size=3><BR>
<DD>Dear itk-users,<BR><BR>
<DD>I did some tests with the underlying random generator of the
itkImageRandomIteratorWithIndex and it seems that, in Windows, it is 'not
very random'. <BR><BR>
<DD>The itkImageRandomIterator uses the following random-generator:<BR>
<DD><ITKSOURCE>\Utilities\vxl\core\vnl\vnl_sample.<h/cxx><BR><BR>
<DD>In Linux the drand48 random-generator is used, which is a good
choice.<BR>
<DD>In Windows however a "simple congruential random number generator" is
implemented, since drand48 is not available in Windows. This gives
inferior results, in my experience.<BR><BR>
<DD>To get a feel for the result look at the image1.<mhd/raw>, which
is in the file: randomtestresults.zip, which you can download from: <A
href="http://www.isi.uu.nl/People/Stefan/"
eudora="autourl">http://www.isi.uu.nl/People/Stefan/</A><BR>
<DD>The gray-values in this image show how many times a pixel was sampled.
The sampling<BR>
<DD>process works was defined as follows:<BR>
<DD> "An ImageRegionIterator walks N=200 times through the
image and tests<BR>
<DD> at each voxel whether to sample it or not. The test is
performed by <BR>
<DD> drawing a number between 0 and 1 using the random
generator defined in<BR>
<DD> vnl_sample.h; a value >=0.5 means that the voxel is
sampled"<BR>
<DD>As you can see in the image the pixels are not really selected at
random...<BR><BR>
<DD>If we use the itkImageRandomIteratorWithIndex to sample the image, the
result (image3.<mhd,raw>) may look better at first sight, but the
histogram looks terrible.<BR><BR>
<DD>On the internet I found the following link:<BR>
<DD><A
href="http://paine.wiau.man.ac.uk/pub/doc_vxl/core/vnl/html/classvnl__random.html"
eudora="autourl">http://paine.wiau.man.ac.uk/pub/doc_vxl/core/vnl/html/classvnl__random.html</A><BR>
<DD>It describes the files vnl_random.<h,cxx>; those files are not
included in the vxl-version that comes with ITK. I tried to use these as a
random generator. The results are now much better. Image2, which is
generated in the same was as image1, but with the random generator defined
in vnl_random, does not show any structure anymore. Image4, which is
generated using an ImageRandomIterator that uses the vnl_random as its
underlying random number generator, has the histogram that you would
expect.<BR><BR>
<DD>For more details please look at the source of the test program, which
you can also download from <A href="http://www.isi.uu.nl/People/Stefan/"
eudora="autourl">http://www.isi.uu.nl/People/Stefan/</A> :
itkrandomtestsource.zip. You may reproduce the results with this program.
Note that only in Windows bad results will be obtained. In Linux there is
no problem.<BR><BR>
<DD>Algorithms that rely on a good random generator may fail if you use
the itkRandomImageIteratorWithIndex under Windows. For example, in my
research on nonrigid registration I tried to use a stochastic gradient
descent optimisation method for minimising the MattesMutualInformation in
combination with a B-spline transform (it may speed up your registration
algorithms; if you are interested: <BR>
<DD><A
href="http://www.isi.uu.nl/Research/Publications/publicationview.php?id=1011"
eudora="autourl">http://www.isi.uu.nl/Research/Publications/publicationview.php?id=1011</A>
). When I used the itkImageRandomConstIteratorWithIndex for selecting
spatial samples the registration results got significantly worse than when
I used the itkImageMoreRandomConstIteratorWithIndex (which uses the
vnl_random as underlying random number generator).<BR><BR>
<DD>Any comments on this would be appreciated!<BR><BR>
<DD>Stefan<BR><BR><BR><BR><BR><BR></FONT></DD></DL></BLOCKQUOTE><BR></BLOCKQUOTE></BODY></HTML>