[vtkusers] vtkMath::Random vs vtkMinimalStandardRandomSequence
Francois Bertel
francois.bertel at kitware.com
Wed Dec 9 14:10:09 EST 2009
BTW, ITK has the same defect that vtkMath::Random() regarding the
side-effect behavior: GetVariate() go the next value in the random
sequence AND return the new value. So if in ITK you have code like
f(seq->GetVariate(), seq->GetVariate()) you are in trouble even if seq
is an object.
On Wed, Dec 9, 2009 at 2:00 PM, Francois Bertel
<francois.bertel at kitware.com> wrote:
> ... and I just saw a spelling mistake that I have to fix :-)
>
> On Wed, Dec 9, 2009 at 1:58 PM, Francois Bertel
> <francois.bertel at kitware.com> wrote:
>> It has not be announced but it is documented in vtkMath::Random() and
>> others with:
>>
>> // DON'T USE Random(), RandomSeed(), GetSeed(), Gaussian()
>> // THIS IS STATIC SO THIS IS PRONE TO ERRORS (SPECIAL FOR REGRESSION TESTS)
>> // THIS IS HERE FOR BACKWARD COMPATIBILITY ONLY.
>> // Instead, for a sequence of random numbers with a uniform distribution
>> // create a vtkMinimalStandardRandomSequence object.
>> // For a sequence of random numbers with a gaussian/normal distribution
>> // create a vtkBoxMuellerRandomSequence object.
>>
>>
>> On Wed, Dec 9, 2009 at 1:54 PM, Bill Lorensen <bill.lorensen at gmail.com> wrote:
>>> OK, so where is all this documented/announced?
>>>
>>> On Wed, Dec 9, 2009 at 1:46 PM, Francois Bertel
>>> <francois.bertel at kitware.com> wrote:
>>>> vtkMath::Random() has been refactored to use
>>>> vtkMinimalStandardRandomSequence, a static reference a
>>>> vtkMinimalStandardRandomSequence object.
>>>>
>>>> vtkMath::Random() is still there for backward compatibility reason.
>>>>
>>>>
>>>> But it does not matter how Random() is implemented: this is just a
>>>> function, not a class so *by construction* it has to use a static
>>>> storage, (or you have to change Random to take a seed in argument, but
>>>> you lose the only reason to keep Random() around, the backward
>>>> compatibility).
>>>>
>>>> Regarding a different algorithm for random generation,
>>>> vtkMinimalStandardRandomSequence
>>>> is actually a subclass of vtkRandomSequence (an abstract random
>>>> sequence), so it is possible to have another subclass of
>>>> vtkRandomSequence, vtkMersenneTwisterRandomVariateSequence with the a
>>>> new implementation of the algorithm. In this case, the implementation
>>>> does not matter much. The discussion was about the architecture
>>>> (static storage vs sequence object).
>>>>
>>>> On Wed, Dec 9, 2009 at 1:38 PM, Bill Lorensen <bill.lorensen at gmail.com> wrote:
>>>>> BTW, here is the random number generator we use in itk:
>>>>>
>>>>> /** \class MersenneTwisterRandomVariateGenerator
>>>>> * \brief MersenneTwisterRandom random variate generator
>>>>> *
>>>>> * \warning This class is NOT thread-safe.
>>>>> *
>>>>> * This notice was included with the original implementation.
>>>>> * The only changes made were to obfuscate the author's email addresses.
>>>>> *
>>>>> * MersenneTwister.h
>>>>> * Mersenne Twister random number generator -- a C++ class MTRand
>>>>> * Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
>>>>> * Richard J. Wagner v1.0 15 May 2003 rjwagner at writeme dot com
>>>>> *
>>>>> * The Mersenne Twister is an algorithm for generating random numbers. It
>>>>> * was designed with consideration of the flaws in various other generators.
>>>>> * The period, 2^19937-1, and the order of equidistribution, 623 dimensions,
>>>>> * are far greater. The generator is also fast; it avoids multiplication and
>>>>> * division, and it benefits from caches and pipelines. For more information
>>>>> * see the inventors' web page at http:*www.math.keio.ac.jp/~matumoto/emt.html
>>>>>
>>>>> * Reference
>>>>> * M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally
>>>>> * Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions on
>>>>> * Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3-30.
>>>>> *
>>>>>
>>>>>
>>>>> On Wed, Dec 9, 2009 at 1:35 PM, Bill Lorensen <bill.lorensen at gmail.com> wrote:
>>>>>> Why dpn't we change vtkMath::Random's implementation so that it uses
>>>>>> vtkMinimalStandardRandomSequence.
>>>>>>
>>>>>> On Wed, Dec 9, 2009 at 1:31 PM, Francois Bertel
>>>>>> <francois.bertel at kitware.com> wrote:
>>>>>>> Don't use drand48() but DON'T USE vtkMath::Random() ANYMORE. Use a
>>>>>>> vtkMinimalStandardRandomSequence object.
>>>>>>>
>>>>>>>
>>>>>>> Even if you use vtkMath::Random() , the following line is not portable
>>>>>>> across compilers because
>>>>>>> the random number generator uses a static variable and the evaluation
>>>>>>> order of the argument of a method/function is compiler dependent:
>>>>>>>
>>>>>>> points->InsertNextPoint(drand48(), drand48(), drand48());
>>>>>>>
>>>>>>> Imaging that drand48() or vtkMath::Random() generates the following
>>>>>>> sequence of three numbers: 0.1 0.2 0.3 (randomly chosen by myself :-)
>>>>>>> )
>>>>>>>
>>>>>>> If the compiler evaluates the last argument first, the line will be
>>>>>>> equivalent to
>>>>>>> points->InsertNextPoint(0.3,0.2,0.1);
>>>>>>>
>>>>>>> Another compiler can decide to evaluate the first argument first, the
>>>>>>> line will be equivalent to
>>>>>>> points->InsertNextPoint(0.1,0.2,0.3);
>>>>>>>
>>>>>>> Notably, gcc and Visual Studio uses different order.
>>>>>>>
>>>>>>> It can also differ from the same compiler if the optimization flags
>>>>>>> are activated or not.
>>>>>>>
>>>>>>> One way to ensure the order of evaluation is to force a sequence of calls:
>>>>>>> double x=vtkMath::Random(); // 0.1
>>>>>>> double y=vtkMath::Random(); // 0.2
>>>>>>> double z=vtkMath::Random(); // 0.3
>>>>>>>
>>>>>>> points->InsertNextPoint(x,y,z);
>>>>>>>
>>>>>>> But I just said to not use vtkMath::Random() ANYMORE. The replacement
>>>>>>> for vtkMath::Random() and others is the class
>>>>>>> Common/vtkMinimalStandardRandomSequence:
>>>>>>>
>>>>>>> This class as no nasty side-effect: GetValue() just get the current
>>>>>>> value of the random sequence, it does not go the next random value.
>>>>>>> Next() just move the the next number in the random sequence but does
>>>>>>> not return the new value. (Command-query separation principle).
>>>>>>>
>>>>>>> Therefore the way the class is written force you to break down the
>>>>>>> calls so you magically avoid the compiler specific behavior:
>>>>>>>
>>>>>>> vtkMinimalStandardRandomSequence *seq=vtkMinimalStandardRandomSequence::New();
>>>>>>> seq->SetSeed(1); // initialize the sequence
>>>>>>> double x=seq->GetValue(); // or seq->GetRangeValue(-1.0,1.0); in your case
>>>>>>> seq->Next();
>>>>>>> double y=seq->GetValue();
>>>>>>> seq->Next();
>>>>>>> double z=seq->GetValue();
>>>>>>> points->InsertNextPoint(x,y,z);
>>>>>>> ...
>>>>>>> ...
>>>>>>>
>>>>>>> seq->Delete();
>>>>>>>
>>>>>>> The other great thing with sequence objects is you can have multiple
>>>>>>> independent random sequences. Some classes use random numbers
>>>>>>> internally for diverse reasons (it can be sampling/jittering, cheery
>>>>>>> picking). With vtkMath::Random() using a static seed, your random
>>>>>>> number generation could interfere with the random number generation of
>>>>>>> an other component. (and all the other drawbacks of static storage
>>>>>>> you can image, like issues with multithreading)
>>>>>>>
>>>>>>>
>>>>>>> On Wed, Dec 9, 2009 at 1:11 PM, David Doria <daviddoria+vtk at gmail.com> wrote:
>>>>>>>> On Wed, Dec 9, 2009 at 12:38 PM, Bill Lorensen <bill.lorensen at gmail.com> wrote:
>>>>>>>>> David,
>>>>>>>>>
>>>>>>>>> 1) Please do not use drand48. It's not portable. Use
>>>>>>>>> vtkMath::Random(0.0,1.0) for equivalent.
>>>>>>>>> 2) drand48 generates numbers between 0,1. You cannot extract an
>>>>>>>>> isosurface of 0. Try vtkMath::Random(-1.0,1.0)
>>>>>>>>>
>>>>>>>>> Bill
>>>>>>>>
>>>>>>>> Bill,
>>>>>>>>
>>>>>>>> 1) Good point. Will do.
>>>>>>>>
>>>>>>>> 2) You are correct, that was an artifact of simplifying my code into a
>>>>>>>> demo. In the demo, when this was changed, it worked as expected. In my
>>>>>>>> real code, the values were indeed on both sides of zero.
>>>>>>>>
>>>>>>>> I have isolated the problem: I was setting up the ImageData in a
>>>>>>>> separate class that derives from vtkImageAlgorithm (aside: why is it
>>>>>>>> not called vtkImageDataAlgorithm?)
>>>>>>>>
>>>>>>>> When I cut and paste the contents of the RequestData function from
>>>>>>>> that class into the function we have been discussing, it works fine.
>>>>>>>>
>>>>>>>> I checked the input pointer in RequestData and it has the correct
>>>>>>>> number of points, so the input of the pipeline has been setup
>>>>>>>> correctly. The output of the voxelizeFilter seems to be setup properly
>>>>>>>> back in the calling function. It seems to me that these two should
>>>>>>>> behave identically, but one produces valid output while the other
>>>>>>>> produces empty output.
>>>>>>>>
>>>>>>>> Here is the VoxelizePolyData class:
>>>>>>>> http://www.rpi.edu/~doriad/VTK_List/vtkVoxelizePolyData/
>>>>>>>>
>>>>>>>> Does anyone see why these two things would be any different:
>>>>>>>>
>>>>>>>> A)
>>>>>>>> vtkSmartPointer<vtkVoxelizePolyData> voxelizeFilter =
>>>>>>>> vtkSmartPointer<vtkVoxelizePolyData>::New();
>>>>>>>> voxelizeFilter->SetInput(input);
>>>>>>>> voxelizeFilter->Update();
>>>>>>>> vtkImageData* grid = voxelizeFilter->GetOutput();
>>>>>>>>
>>>>>>>> B)
>>>>>>>>
>>>>>>>> //get the bounds of the input PolyData
>>>>>>>> double bounds[6];
>>>>>>>> input->GetBounds(bounds);
>>>>>>>> //xmin, xmax, ymin, ymax, zmin, zmax
>>>>>>>> double xmin = bounds[0];
>>>>>>>> double xmax = bounds[1];
>>>>>>>> double ymin = bounds[2];
>>>>>>>> double ymax = bounds[3];
>>>>>>>> double zmin = bounds[4];
>>>>>>>> double zmax = bounds[5];
>>>>>>>>
>>>>>>>> //create a grid the size of the point set
>>>>>>>> vtkSmartPointer<vtkImageData> grid = vtkSmartPointer<vtkImageData>::New();
>>>>>>>>
>>>>>>>> //set the bottom left corner of the grid to the bottom left corner
>>>>>>>> of the data bounding box
>>>>>>>> grid->SetOrigin(xmin, ymin, zmin);
>>>>>>>>
>>>>>>>> //create a grid of voxels of size (NumberOfCellsX x NumberOfCellsY x
>>>>>>>> NumberOfCellsZ)
>>>>>>>> //there are (NumberOfCellsX+1 x NumberOfCellsY+1 x NumberOfCellsZ+1)
>>>>>>>> points to make this size voxel grid
>>>>>>>> grid->SetExtent(0, 10, 0, 10, 0, 10);
>>>>>>>>
>>>>>>>> //set the size of each element/cell/voxel so that the grid spans the
>>>>>>>> entire input point set volume
>>>>>>>> //grid->SetSpacing((xmax-xmin)/static_cast<double>(this->NumberOfCellsX),
>>>>>>>> (ymax-ymin)/static_cast<double>(this->NumberOfCellsY),
>>>>>>>> (zmax-zmin)/static_cast<double>(this->NumberOfCellsZ));
>>>>>>>> grid->SetSpacing((xmax-xmin)/static_cast<double>(10),
>>>>>>>> (ymax-ymin)/static_cast<double>(10),
>>>>>>>> (zmax-zmin)/static_cast<double>(10));
>>>>>>>>
>>>>>>>> Have I done something wrong in setting up the pipeline connections of
>>>>>>>> the voxelizePolyData filter?
>>>>>>>>
>>>>>>>> Thanks,
>>>>>>>>
>>>>>>>> David
>>>>>>>> _______________________________________________
>>>>>>>> Powered by www.kitware.com
>>>>>>>>
>>>>>>>> Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html
>>>>>>>>
>>>>>>>> Please keep messages on-topic and check the VTK FAQ at: http://www.vtk.org/Wiki/VTK_FAQ
>>>>>>>>
>>>>>>>> Follow this link to subscribe/unsubscribe:
>>>>>>>> http://www.vtk.org/mailman/listinfo/vtkusers
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> --
>>>>>>> François Bertel, PhD | Kitware Inc. Suite 204
>>>>>>> 1 (518) 371 3971 x113 | 28 Corporate Drive
>>>>>>> | Clifton Park NY 12065, USA
>>>>>>> _______________________________________________
>>>>>>> Powered by www.kitware.com
>>>>>>>
>>>>>>> Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html
>>>>>>>
>>>>>>> Please keep messages on-topic and check the VTK FAQ at: http://www.vtk.org/Wiki/VTK_FAQ
>>>>>>>
>>>>>>> Follow this link to subscribe/unsubscribe:
>>>>>>> http://www.vtk.org/mailman/listinfo/vtkusers
>>>>>>>
>>>>>>
>>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> François Bertel, PhD | Kitware Inc. Suite 204
>>>> 1 (518) 371 3971 x113 | 28 Corporate Drive
>>>> | Clifton Park NY 12065, USA
>>>> _______________________________________________
>>>> Powered by www.kitware.com
>>>>
>>>> Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html
>>>>
>>>> Please keep messages on-topic and check the VTK FAQ at: http://www.vtk.org/Wiki/VTK_FAQ
>>>>
>>>> Follow this link to subscribe/unsubscribe:
>>>> http://www.vtk.org/mailman/listinfo/vtkusers
>>>>
>>>
>>
>>
>>
>> --
>> François Bertel, PhD | Kitware Inc. Suite 204
>> 1 (518) 371 3971 x113 | 28 Corporate Drive
>> | Clifton Park NY 12065, USA
>>
>
>
>
> --
> François Bertel, PhD | Kitware Inc. Suite 204
> 1 (518) 371 3971 x113 | 28 Corporate Drive
> | Clifton Park NY 12065, USA
>
--
François Bertel, PhD | Kitware Inc. Suite 204
1 (518) 371 3971 x113 | 28 Corporate Drive
| Clifton Park NY 12065, USA
More information about the vtkusers
mailing list