[vtkusers] vtkContourFilter issue
Francois Bertel
francois.bertel at kitware.com
Wed Dec 9 13:46:10 EST 2009
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
More information about the vtkusers
mailing list