[vtkusers] vtkContourFilter issue
David Doria
daviddoria+vtk at gmail.com
Wed Dec 9 13:41:25 EST 2009
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
>>>
Can we fork this random number discussion so this thread can remain
about the ImageData derived RequestData?
Thanks,
David
More information about the vtkusers
mailing list