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