Ok, I posted this once before, but here it is again with some more information.<br>
<br>
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.<br>
<br>
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 :<br>
<div style="margin-left: 40px; font-family: courier new,monospace;">p Array [0]: vtkDoubleArray "result", 2,097,152 points, 8 bytes per point<br>
p
Magnitude Range: 0.0000
1482912.5793<br>
p
Component 0 Range:
-1048575.9848 1048575.8760<br>
p
Component 1 Range:
-1048572.9243 1048579.0968<br>
</div>
With a simple sin wave, the result should be almost purely imaginary, so the Component 0 range should be almost 0 to 0. <br>
<br>
When I switch to using the vtkImageExtractComponents, the results look better:<br>
<div style="font-family: courier new,monospace;">
<div style="margin-left: 40px; font-family: courier new,monospace;">p Array [0]: vtkDoubleArray "image", 2,097,152 points, 8 bytes per point<br>
p
Magnitude Range: 0.0000
1048576.0091<br>
p
Component 0 Range: -0.0557 0.0011<br>
p
Component 1 Range:
-1048576.0089 1048576.0091<br>
</div>
<br>
</div>
<br>
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.<br>
<br>
So, any ideas?<br>
<br>
<span style="font-family: courier new,monospace;">----- Code starts here -------</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"></span>
<div style="margin-left: 40px;"><span style="font-family: courier new,monospace;">#include <stdio.h></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include <string.h></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include <stdlib.h></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include <unistd.h></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include "stat.h"</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include <sys/types.h></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include <sys/stat.h></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include <sys/mman.h></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include <fcntl.h></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include <vtkObject.h></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include "vtkImageData.h"</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include "vtkPointData.h"</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include <vtkDataArray.h></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include <vtkUnsignedCharArray.h></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include <vtkDoubleArray.h></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include <vtkImageData.h></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include <vtkImageFFT.h></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include <vtkArrayCalculator.h></span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#include <vtkImageExtractComponents.h></span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">int main (int argc, char *argv[])</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">{</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> double *data;</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> int xdim, ydim, zdim;</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> double kx;</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> int x,y,z;</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> long index;</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> xdim = 128;</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> ydim = 128;</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> zdim = 128;</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> if ((data = (double*)malloc(sizeof(double)*xdim*ydim*zdim*3)) == NULL) {</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> perror("Unable to allocate memory!\n\t");</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> return -1;</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> }</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> printf("Constructing sin(kx * x) dataset...\n");</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> kx = 2.0 * 3.1415926 / (double)xdim;</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> for(x=0; x<xdim; x++) {</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> for(y=0; y<ydim; y++) {</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> for(z=0; z<zdim; z++) {</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">
index = z + (y*zdim) + (x * ydim * zdim);</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">
data[(index*3)+0] = sin(kx * (double)x);</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">
data[(index*3)+1] = cos(kx * (double)x);</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">
data[(index*3)+2] = tan(kx * (double)x);</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> }</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> }</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> }</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> vtkDoubleArray *ucPointer = vtkDoubleArray::New();</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> ucPointer->SetNumberOfComponents(3);</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> ucPointer->SetArray(data, xdim*ydim*zdim*3, 1);</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> ucPointer->SetName("image");</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> vtkImageData *image = vtkImageData::New();</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> image->Initialize();</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> image->SetDimensions(xdim, ydim, zdim);</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> image->SetExtent(1, xdim, 1, ydim, 1, zdim);</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> image->SetScalarTypeToDouble();</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> image->SetNumberOfScalarComponents(3);</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> image->GetPointData()->SetScalars(ucPointer);</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> PrintStatistics(image);</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> printf("Extracting Velocity X-Component...\n");</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#if 0</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> vtkArrayCalculator *extract = vtkArrayCalculator::New();</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> extract->SetInput(image);</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> extract->SetResultArrayName("result");</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> extract->AddScalarVariable("x", "image", 0);</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> extract->SetFunction("(x*1.0)+0.0");</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> extract->ReplaceInvalidValuesOff();</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> extract->Update();</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> extract->GetOutput()->GetPointData()->RemoveArray("image");</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> extract->GetOutput()->GetPointData()->SetActiveScalars("result");</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#else</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> vtkImageExtractComponents *extract = vtkImageExtractComponents::New();</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> extract->SetInput(image);</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> extract->SetComponents(0);</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> extract->Update();</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">#endif</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> PrintStatistics(extract->GetOutput());</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> printf("Computing FFT...\n");</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> vtkImageFFT *fft = vtkImageFFT::New();</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> fft->SetInput(extract->GetOutput());</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> fft->Update();</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> PrintStatistics(fft->GetOutput());</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">}</span><br>
<br>
</div>
<br>-- <br>Randall Hand<br>Visualization Scientist, <br>ERDC-MSRC Vicksburg, MS<br>Homepage: <a href="http://www.yeraze.com">http://www.yeraze.com</a>