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

[Berk Geveci]

This is a guess but if the array calculator outputs a data type than
the input's scalars, no imaging filter can be used down the pipeline.
This is a known bug. Essentially no (non-imaging) filter that changes
the data type of the scalars is allowed in a pipeline with imaging
filters. There is work (slow) going on to fix this issue. I would not
hold my breath for it's completion.

On 9/14/05, Randall Hand <randall.hand at gmail.com> 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 
> _______________________________________________
> vtk-developers mailing list
> vtk-developers at vtk.org
> http://www.vtk.org/mailman/listinfo/vtk-developers
> 
> 
> 
>