[vtk-developers] Re: [vtkusers] Problem with using vtkArrayCalculator to extract fields for vtkImageFFT Processing

Mathieu Malaterre mathieu.malaterre at kitware.com
Wed Sep 14 13:51:56 EDT 2005


Ok this is a bug:

i ========== End of Data Statistics ===============
==30679==
==30679== ERROR SUMMARY: 1 errors from 1 contexts (suppressed: 9 from 4)
==30679==
==30679== 1 errors in context 1 of 1:
==30679== Invalid read of size 8
==30679==    at 0x12A4B8F5: void 
vtkImageFFTExecute<double>(vtkImageFFT*, vtkImageData*, int*, double*, 
vtkImageData*, int*, double*, int) (vtkImageFFT.cxx:137)
==30679==    by 0x12A4A645: vtkImageFFT::ThreadedExecute(vtkImageData*, 
vtkImageData*, int*, int) (vtkImageFFT.cxx:203)
==30679==    by 0x1353AEA2: 
vtkThreadedImageAlgorithm::ThreadedRequestData(vtkInformation*, 
vtkInformationVector**, vtkInformationVector*, vtkImageData***, 
vtkImageData**, int*, int) (vtkThreadedImageAlgorithm.cxx:327)
==30679==    by 0x1353AA8E: 
vtkThreadedImageAlgorithmThreadedExecute(void*) 
(vtkThreadedImageAlgorithm.cxx:200)
==30679==    by 0x13BA30F9: start_thread (in /lib64/tls/libpthread-2.3.5.so)
==30679==    by 0x1561AA62: clone (in /lib64/tls/libc-2.3.5.so)
==30679==  Address 0x19992030 is 0 bytes after a block of size 16777216 
alloc'd
==30679==    at 0x11B1A8E2: operator new[](unsigned long) (in 
/usr/lib64/valgrind/vgpreload_memcheck.so)
==30679==    by 0x1380B136: vtkDataArrayTemplate<double>::Allocate(long 
long, long long) (vtkDataArrayTemplate.txx:99)
==30679==    by 0x1380BF2E: 
vtkDataArrayTemplate<double>::SetNumberOfValues(long long) 
(vtkDataArrayTemplate.txx:536)
==30679==    by 0x1380B409: 
vtkDataArrayTemplate<double>::SetNumberOfTuples(long long) 
(vtkDataArrayTemplate.txx:330)
==30679==    by 0x12154566: 
vtkArrayCalculator::RequestData(vtkInformation*, vtkInformationVector**, 
vtkInformationVector*) (vtkArrayCalculator.cxx:279)
==30679==    by 0x1341A5AA: 
vtkDataSetAlgorithm::ProcessRequest(vtkInformation*, 
vtkInformationVector**, vtkInformationVector*) 
(vtkDataSetAlgorithm.cxx:175)==30679==    by 0x1342CBC9: 
vtkExecutive::CallAlgorithm(vtkInformation*, int, 
vtkInformationVector**, vtkInformationVector*) (vtkExecutive.cxx:679)
==30679==    by 0x13427AF7: 
vtkDemandDrivenPipeline::ExecuteData(vtkInformation*, 
vtkInformationVector**, vtkInformationVector*) 
(vtkDemandDrivenPipeline.cxx:464)
==30679==    by 0x1342732B: 
vtkDemandDrivenPipeline::ProcessRequest(vtkInformation*, int, 
vtkInformationVector**, vtkInformationVector*) 
(vtkDemandDrivenPipeline.cxx:244)
==30679==    by 0x1352487E: 
vtkStreamingDemandDrivenPipeline::ProcessRequest(vtkInformation*, int, 
vtkInformationVector**, vtkInformationVector*) 
(vtkStreamingDemandDrivenPipeline.cxx:144)
==30679==    by 0x13427927: vtkDemandDrivenPipeline::UpdateData(int) 
(vtkDemandDrivenPipeline.cxx:406)
==30679==    by 0x13524A04: 
vtkStreamingDemandDrivenPipeline::Update(int) 
(vtkStreamingDemandDrivenPipeline.cxx:190)


Could you open a bug in the bugtracker(*) for this. Add your program as 
an attachment.

Thanks
Mathieu
(*) http://vtk.org/Bug



Randall Hand wrote:
> Ok, I posted this once before, but here it is again with some more 
> information.
> 
> I've written a piece of code to create a synthetic dataset.  It's a 
> 128^3 dataset, each point with a 3-component double value consisting of 
> (sin(x), junk, morejunk).  I extract just the first component, and do an 
> FFT.  You'll find the code at the end of this message, along with my 
> "PrintStatistics" function attached.
> 
> You'll see from the #if/#endif blocks that there's two possible ways to 
> do it.  The first is with a vtkArrayCalculator, where I set a very 
> trivial function to simply copy out the one value.  The other way is 
> with a vtkImageExtractComponents.  When doing it with the 
> vtkArrayCalculator, the filter seems to work just fine but then the FFT 
> yields values like :
> p       Array [0]: vtkDoubleArray "result", 2,097,152 points, 8 bytes 
> per point
> p                       Magnitude Range:        0.0000  1482912.5793
> p                       Component 0 Range:      -1048575.9848   1048575.8760
> p                       Component 1 Range:      -1048572.9243   1048579.0968
> With a simple sin wave, the result should be almost purely imaginary, so 
> the Component 0 range should be almost 0 to 0. 
> 
> When I switch to using the vtkImageExtractComponents, the results look 
> better:
> p       Array [0]: vtkDoubleArray "image", 2,097,152 points, 8 bytes per 
> point
> p                       Magnitude Range:        0.0000  1048576.0091
> p                       Component 0 Range:      -0.0557 0.0011
> p                       Component 1 Range:      -1048576.0089   1048576.0091
> 
> 
> Unfortunately, when working with the vtkImageExtractComponents filter, I 
> lose the other two fields if I had planned on using them for coloring or 
> something.  Also, my results aren't quite what I expected.  Much better, 
> but on more complex datasets I seem to get results that are a little off 
> from the expected result.  I don't know if the cause of all this is 
> related, but I'm chasing down what I can.
> 
> So, any ideas?
> 
> ----- Code starts here -------
> #include <stdio.h>
> #include <string.h>
> #include <stdlib.h>
> #include <unistd.h>
> #include "stat.h"
> #include <sys/types.h>
> #include <sys/stat.h>
> #include <sys/mman.h>
> #include <fcntl.h>
> #include <vtkObject.h>
> #include "vtkImageData.h"
> #include "vtkPointData.h"
> #include <vtkDataArray.h>
> #include <vtkUnsignedCharArray.h>
> #include <vtkDoubleArray.h>
> #include <vtkImageData.h>
> #include <vtkImageFFT.h>
> #include <vtkArrayCalculator.h>
> #include <vtkImageExtractComponents.h>
> 
> int main (int argc, char *argv[])
> {
>     double *data;
>     int xdim, ydim, zdim;
>     double kx;
>     int x,y,z;
>     long index;
>     xdim = 128;
>     ydim = 128;
>     zdim = 128;
> 
>     if ((data = (double*)malloc(sizeof(double)*xdim*ydim*zdim*3)) == NULL) {
>         perror("Unable to allocate memory!\n\t");
>         return -1;
>     }
> 
>     printf("Constructing sin(kx * x) dataset...\n");
>     kx = 2.0 * 3.1415926 / (double)xdim;
>     for(x=0; x<xdim; x++) {
>         for(y=0; y<ydim; y++) {
>             for(z=0; z<zdim; z++) {
>                 index = z + (y*zdim) + (x * ydim * zdim);
>                 data[(index*3)+0] = sin(kx * (double)x);
>                 data[(index*3)+1] = cos(kx * (double)x);
>                 data[(index*3)+2] = tan(kx * (double)x);
>             }
>         }
>     }
> 
>     vtkDoubleArray *ucPointer = vtkDoubleArray::New();
>     ucPointer->SetNumberOfComponents(3);
>     ucPointer->SetArray(data, xdim*ydim*zdim*3, 1);
>     ucPointer->SetName("image");
> 
>     vtkImageData *image = vtkImageData::New();
>     image->Initialize();
>     image->SetDimensions(xdim, ydim, zdim);
>     image->SetExtent(1, xdim, 1, ydim, 1, zdim);
>     image->SetScalarTypeToDouble();
>     image->SetNumberOfScalarComponents(3);
>     image->GetPointData()->SetScalars(ucPointer);
>     PrintStatistics(image);
>     printf("Extracting Velocity X-Component...\n");
> 
> #if 0
>     vtkArrayCalculator *extract = vtkArrayCalculator::New();
>     extract->SetInput(image);
>     extract->SetResultArrayName("result");
>     extract->AddScalarVariable("x", "image", 0);
>     extract->SetFunction("(x*1.0)+0.0");
>     extract->ReplaceInvalidValuesOff();
>     extract->Update();
>     extract->GetOutput()->GetPointData()->RemoveArray("image");
>     extract->GetOutput()->GetPointData()->SetActiveScalars("result");
> #else
>     vtkImageExtractComponents *extract = vtkImageExtractComponents::New();
>     extract->SetInput(image);
>     extract->SetComponents(0);
>     extract->Update();
> #endif
>     PrintStatistics(extract->GetOutput());
> 
>     printf("Computing FFT...\n");
>     vtkImageFFT *fft = vtkImageFFT::New();
>     fft->SetInput(extract->GetOutput());
>     fft->Update();
>     PrintStatistics(fft->GetOutput());
> 
> }
> 
> 
> -- 
> Randall Hand
> Visualization Scientist,
> ERDC-MSRC Vicksburg, MS
> Homepage: http://www.yeraze.com
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> This is the private VTK discussion list. 
> Please keep messages on-topic. Check the 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 vtk-developers mailing list