[vtkusers] vtkMath::Random vs vtkMinimalStandardRandomSequence
Francois Bertel
francois.bertel at kitware.com
Wed Dec 9 14:00:20 EST 2009
... 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
More information about the vtkusers
mailing list