[vtkusers] vtkContourFilter issue
Bill Lorensen
bill.lorensen at gmail.com
Wed Dec 9 13:38:20 EST 2009
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
>>
>
More information about the vtkusers
mailing list