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

[Randall Hand]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtk-developers/attachments/20050914/217f3dd2/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: stat.h
Type: application/octet-stream
Size: 7396 bytes
Desc: not available
URL: <http://public.kitware.com/pipermail/vtk-developers/attachments/20050914/217f3dd2/attachment-0001.obj>