[vtkusers] A Problem Viewing 3D vtkImageData

[Alex Southern]

  Hi,

Im trying to view a regular 3D grid of points using vtkImageData.   Each 
point is a single scalar that is mapped to a color.
I've included the complete source code and cmakelists.txt file of a 
simplified problem scenario.

Basically, can someone confirm that they also do not see a cuboid but 
instead a 2D plane, why is this?

In addition,  the render window is not updated after each iteration as I 
intended. What am I doing wrong?

(Note the source is called ImageSlicing.cpp, but I am not at the stage 
of the slicing yet, I want to view the whole 3D pressure field first)

Thanks
Alex

// ----------- CmakeLists.txt ------------------

cmake_minimum_required(VERSION 2.6)

PROJECT(ImageSlicing)

FIND_PACKAGE(VTK REQUIRED)
INCLUDE(${VTK_USE_FILE})

ADD_EXECUTABLE(ImageSlicing ImageSlicing.cpp)
TARGET_LINK_LIBRARIES(ImageSlicing vtkHybrid vtkWidgets)

//------------------------------ End Cmakelists.txt

// --------------------------Start ImageSlicing.cpp 
----------------------------------------
#include "vtkSmartPointer.h"
#include "vtkImageReslice.h"
#include "vtkLookupTable.h"
#include "vtkImageMapToColors.h"
#include "vtkImageActor.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkInteractorStyleImage.h"
#include "vtkCommand.h"
#include "vtkImageData.h"
#include "vtkPointData.h"
#include "vtkProperty.h"
#include "vtkFloatArray.h";
#include "vtkLookupTable.h";
#include "vtkScalarBarActor.h";
#include "vtkImageMapToColors.h";

int RunMockSimulationStep(int L, int W, int H, int offset1, int offset2, 
float* p, int numPoints)
{

         for(int z = 1; z < H-1;z++)
         {
             for(int y = 1; y < W-1;y++)
             {
                 for(int x = 1; x < L-1;x++)
                 {
                     p[offset1+(z*L*W)+(y*W)+x] = (1.0/6.0)*
                         (p[offset2+((z+1)*L*W)+(y*W)+x] +
                             p[offset2+((z-1)*L*W)+(y*W)+x] +
                                 p[offset2+(z*L*W)+((y+1)*W)+x] +
                                     p[offset2+(z*L*W)+((y-1)*W)+x] +
                                         p[offset2+(z*L*W)+(y*W)+x+1] +
                                             p[offset2+(z*L*W)+(y*W)+x-1]) -
                                                 
p[offset1+(z*L*W)+(y*W)+x];
                 }
             }
         }

     return 0;
}

int main (int argc, char *argv[])
{

     // Setup Scalar Data Memory
     int L = 100; // Set Length/Width/Height of Pressure grid
     int W = 50;
     int H = 100;
     int numPoints = L*W*H;
     float* pressure;
     pressure = new float [2*numPoints]; // There are two pressure 
grids; we only want to view the first one though...
     float** viewPressure;
     viewPressure = new float* [numPoints];

     //Create any required VTK objects
     vtkSmartPointer<vtkFloatArray> vtkPressure = 
vtkSmartPointer<vtkFloatArray>::New();
     vtkSmartPointer<vtkImageData> imageData = 
vtkSmartPointer<vtkImageData>::New();
     vtkSmartPointer<vtkLookupTable> table = 
vtkSmartPointer<vtkLookupTable>::New();
     vtkSmartPointer<vtkImageMapToColors> color = 
vtkSmartPointer<vtkImageMapToColors>::New();
     vtkSmartPointer<vtkImageActor> imActor = 
vtkSmartPointer<vtkImageActor>::New();
     vtkSmartPointer<vtkRenderer> ren = vtkSmartPointer<vtkRenderer>::New();
     vtkSmartPointer<vtkRenderWindow> renWin = 
vtkSmartPointer<vtkRenderWindow>::New();
     vtkSmartPointer<vtkRenderWindowInteractor> iren = 
vtkSmartPointer<vtkRenderWindowInteractor>::New();

     // Specify the size of the image data
     imageData->SetDimensions(W,L,H);
     imageData->SetOrigin(0,0,0);
     imageData->SetSpacing(1.0,1.0,1.0);
     imageData->SetNumberOfScalarComponents(1);
     imageData->SetScalarTypeToFloat();
     imageData->AllocateScalars();


     int* dims = imageData->GetDimensions();
     std::cout << "ImageData Dimensions: " << "x: " << dims[0] << " y: " 
<< dims[1] << " z: " << dims[2] << std::endl;

     //Initialize pressure across grid to zero
     for(int z = 0; z < H;z++)
     {
         for(int y = 0; y < W;y++)
         {
             for(int x = 0; x < L;x++)
             {
                 pressure[(z*L*W)+(y*W)+x] = 0.0;
                 pressure[numPoints+(z*L*W)+(y*W)+x] = 0.0;
                 viewPressure[(y*W)+x] = &pressure[(z*L*W)+(y*W)+x];
             }
         }
     }

     // Now point imageData to the pressure vector points
     vtkPressure->SetArray(*viewPressure,L*W,1);
     imageData->GetPointData()->SetScalars(vtkPressure);
     vtkPressure->Modified();

     // Create a scale lookup table
     double r = 0;
     double g = 0;
     double b = 1.0;
     table->SetRange(-0.02, 0.02); // image intensity range
     table->SetNumberOfTableValues(64);
     for (int x = 0; x < 64 ; x++)
         {
         if (x > 32) {    r = 1.0; g = g - 0.0312; b = b - 0.0312;    }
         if (x == 32){    r = 1.0; g = 1.0; b = 1.0;    }
         if (x < 32)    {    r = r + 0.0312;    g = g + 0.0312;    b = 
1.0;    }
         table->SetTableValue(x,r,g,b,1.0);
     }
     table->SetTableValue(32,1,1,1,1.0); // make sure middle value is 
total black
     table->SetRampToLinear();
     table->Build();

     // Map the image through the lookup table
     color->SetLookupTable(table);
     color->SetInput(imageData);
     imActor->SetInput(color->GetOutput());

     // Setup Interactor Renderer, RenderWindow etc...
     ren->AddActor(imActor);
     ren->SetBackground(0.0,0.0,0.0);
     renWin->SetInteractor(iren);
     renWin->AddRenderer(ren);
     renWin->Render();

     //Initial Condition - Set a pressure point to non-zero
     int xi = 20; int yi = 40; int zi = 36;
     pressure[(zi*L*W)+(yi*W)+xi] = 100.0;

     int T = 50; //Set number of time steps
     int offset1 = numPoints;
     int offset2 = 0;
     int tmp = 0;

     for (int t = 0; t < T; t++)
     {

         cout << t << endl;
         RunMockSimulationStep(L,W,H,offset1,offset2,pressure,numPoints);

         //Update Pressure Field View Here
         // How?
         vtkPressure->Modified();
         renWin->Render();

         // cycle offset
         tmp = offset1;
         offset1 = offset2;
         offset2 = tmp;
     }

     system("pause");
     delete []pressure;
     delete []viewPressure;
     return EXIT_SUCCESS;
}

// --------------- END CPP