[vtkusers] vtkMath::Random vs vtkMinimalStandardRandomSequence
Bill Lorensen
bill.lorensen at gmail.com
Wed Dec 9 13:54:50 EST 2009
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
>
More information about the vtkusers
mailing list