How to generate vtkImageData from mathematical function

Support
1

Hello.

I’m trying to generate multiple DataSets containing Cayley surfaces in different levels of detail.
The cayley surface is described by following function:

f(x, y, z) = 1 - 16xyz - 4x^2 - 4y^2 - 4z^2

No I don’t have issues generating the values in C++ but im not quite sure on how to load the generated arrays into vtkDataSets. Eventually I want to save them into MHD Files using MetaImageWriter. I couldn’t find a fitting example in the various online tutorial sources so I’m happy if anybody can point me towards a solution since im assuming that this probably was done before with different functions.
Im aiming to create a small command line tool which creates Cayley surfaces in different sizes to output big test data for my GUI application.

I know this isn’t a VTK issue but the Support Topic also is labled with “Questions about VTK usage” so I tought I’d give it a shot.

Also if anybody knows where I can find publically available MetaImage data I’d be happy for pointers.

Best Regards,
David

2

This is what I did so far:

#define X_MAX 100
#define X_MIN -100
#define Y_MAX 100
#define Y_MIN -100
#define Z_MAX 100
#define Z_MIN -100

uint16_t getValue(uint16_t, uint16_t, uint16_t);

int main(int argc, char* argv[]) {
	int *size = new int[3];

	size[0] = abs(X_MIN - X_MAX);
	size[1] = abs(Y_MIN - Y_MAX);
	size[2] = abs(Z_MIN - Z_MAX);
	
	uint16_t***values = new uint16_t **[size[0]];
	for (uint16_t i = 0; i <= size[1]; ++i) {
		values[i] = new uint16_t *[size[2]];

		for (uint16_t j = 0; j <= size[1]; ++j) {
			values[i][j] = new uint16_t[size[2]];
		}
	}

	int i = 0;
	for (uint16_t x = X_MIN; x <= X_MAX; x++) {
		for (uint16_t y = Y_MIN; y <= Y_MAX; y++) {
			for (uint16_t z = Z_MIN; z <= Z_MAX; z++) {
				uint16_t value = getValue(x, y, z);
				std::cout << "(" << x << "|" << y << "|" << z << ") = " << value << std::endl;
				values[x + abs(X_MIN)][y + abs(Y_MIN)][z + abs(Z_MIN)] = value;
			}
		}
	}

	vtkIdType totalSize = (vtkIdType) size[0] * (vtkIdType) size[1] * (vtkIdType) size[2];

	vtkNew<vtkUnsignedShortArray> array;
	array->SetArray(**values, totalSize, 1);

	vtkNew<vtkInformation> info;

	vtkNew<vtkImageData> imageData;
	imageData->GetPointData()->SetScalars(array);
	imageData->SetDimensions(size); 
	imageData->SetScalarType(VTK_UNSIGNED_SHORT, info);
	imageData->SetSpacing(1.0, 1.0, 1.0); 
	imageData->SetOrigin(0.0, 0.0, 0.0);

	vtkNew<vtkMetaImageWriter> writer;
	writer->SetInputData(imageData);
	writer->SetFileName("data.mhd");
	writer->Update();
}

But this throws an exception:

Ausnahme ausgelöst bei 0x00007FFD5C622BE7 (ucrtbased.dll) in minimal.exe: 0xC0000005: Zugriffsverletzung beim Lesen an Position 0x0000000000000001.

Translating to access violation at 0x0…01 in minimal.exe (ucrtbased.dll). Not sure why this is coming up. The execution stack states that the error is thrown in line 66 (writer->Update();).

Update:

PrintSelf output:

Modified Time: 51
Reference Count: 1
Registered Events: (none)
Information: 0000025C2C3574D0
Data Released: False
Global Release Data: Off
UpdateTime: 0
Field Data:
  Debug: Off
  Modified Time: 16
  Reference Count: 1
  Registered Events: (none)
  Number Of Arrays: 0
  Number Of Components: 0
  Number Of Tuples: 0
Number Of Points: 8000000
Number Of Cells: 7880599
Cell Data:
  Debug: Off
  Modified Time: 19
  Reference Count: 1
  Registered Events:
    Registered Observers:
      vtkObserver (0000025C2C372960)
        Event: 33
        EventName: ModifiedEvent
        Command: 0000025C2C357EF0
        Priority: 0
        Tag: 1
  Number Of Arrays: 0
  Number Of Components: 0
  Number Of Tuples: 0
  Copy Tuple Flags: ( 1 1 1 1 1 0 1 1 1 1 1 )
  Interpolate Flags: ( 1 1 1 1 1 0 0 1 1 1 1 )
  Pass Through Flags: ( 1 1 1 1 1 1 1 1 1 1 1 )
  Scalars: (none)
  Vectors: (none)
  Normals: (none)
  TCoords: (none)
  Tensors: (none)
  GlobalIds: (none)
  PedigreeIds: (none)
  EdgeFlag: (none)
  Tangents: (none)
  RationalWeights: (none)
  HigherOrderDegrees: (none)
Point Data:
  Debug: Off
  Modified Time: 36
  Reference Count: 1
  Registered Events:
    Registered Observers:
      vtkObserver (0000025C2C372780)
        Event: 33
        EventName: ModifiedEvent
        Command: 0000025C2C357EF0
        Priority: 0
        Tag: 1
  Number Of Arrays: 1
  Array 0 name = nullptr
  Number Of Components: 1
  Number Of Tuples: 8000000
  Copy Tuple Flags: ( 1 1 1 1 1 0 1 1 1 1 1 )
  Interpolate Flags: ( 1 1 1 1 1 0 0 1 1 1 1 )
  Pass Through Flags: ( 1 1 1 1 1 1 1 1 1 1 1 )
  Scalars:
    Debug: Off
    Modified Time: 11
    Reference Count: 2
    Registered Events: (none)
    Name: (none)
    Data type: unsigned short
    Size: 8000000
    MaxId: 7999999
    NumberOfComponents: 1
    Information: 0000000000000000
    Name: (none)
    Number Of Components: 1
    Number Of Tuples: 8000000
    Size: 8000000
    MaxId: 7999999
    LookupTable: (none)
  Vectors: (none)
  Normals: (none)
  TCoords: (none)
  Tensors: (none)
  GlobalIds: (none)
  PedigreeIds: (none)
  EdgeFlag: (none)
  Tangents: (none)
  RationalWeights: (none)
  HigherOrderDegrees: (none)
Bounds:
  Xmin,Xmax: (0, 199)
  Ymin,Ymax: (0, 199)
  Zmin,Zmax: (0, 199)
Compute Time: 59
Spacing: (1, 1, 1)
Origin: (0, 0, 0)
Direction: (1, 0, 0, 0, 1, 0, 0, 0, 1)
Dimensions: (200, 200, 200)
Increments: (0, 0, 0)
Extent: (0, 199, 0, 199, 0, 199)
3

On way to simplify the creation of an image is with image iterators, such as vtkImagePointIterator. This iterator will iterate over each pixel/voxel and give you the (x,y,z) coords that you can use to evaluate your function.

I haven’t seen a good example (mea culpa, since I wrote the iterator class), but the usage is roughly as follows:

#include <vtkImageData.h>
#include <vtkImagePointIterator.h>

// set image properties for an image with a domain of 
//  [ -10.0, +10.0] in x, y and z
double spacing[3] = { 0.1, 0.1, 0.1 };
int size[3] = { 201, 201, 201 };
double origin[3] = { -10.0, -10.0, -10.0 };

// create an image and allocate the pixel data
auto imageData = vtkSmartPointer<vtkImageData>::New();
imageData->SetDimensions(size);
imageData->SetSpacing(spacing);
imageData->SetOrigin(origin);
// replace '1' with '2' if you want to create a complex image
imageData->AllocateScalars(VTK_FLOAT, 1);

// get the scalars from the image (they're allocated but not initialized).
vtkFloatArray *scalars = vtkFloatArray::SafeDownCast(imageData->GetPointData()->GetScalars());
// alternatively, use vtkDataArray to avoid the cast
//vtkDataArray *scalars = imageData->GetPointData()->GetScalars();

vtkImagePointIterator iter(imageData);
while (!iter.IsAtEnd())
{
    double point[3];
    iter.GetPosition(point);
    // as an example, use this simple function
    double value = point[0]*point[0] + point[1]*point[1] + point[2]*point[2];
    // insert the value into the image
    scalars->SetValue(iter.GetId(), value);
    iter.Next();
}

// as an optional setp, vtkImageShiftScale can be used with an
// appropriate shift and scale to convert the image from float to
// unsigned short for rendering

Edit: Regarding the crash that you saw, it was probably due to the fact the vtkImageData requires a contiguous array, but your code allocated non-contiguous segments. A vtkImageData would need a single array of length size[0]*size[1]*size[2].

4

Thank you for your example. But using your code I stillget the same crash:

double spacing[3] = { 0.1, 0.1, 0.1 };
int size[3] = { 201, 201, 201 };
double origin[3] = { -10.0, -10.0, -10.0 };

auto imageData = vtkSmartPointer<vtkImageData>::New();
imageData->SetDimensions(size);
imageData->SetSpacing(spacing);
imageData->SetOrigin(origin);
imageData->AllocateScalars(VTK_FLOAT, 1);

vtkFloatArray* scalars = vtkFloatArray::SafeDownCast(imageData->GetPointData()->GetScalars());
vtkImagePointIterator iter(imageData);

while (!iter.IsAtEnd())	{
	double point[3];
	iter.GetPosition(point);
	double value = point[0] * point[0] + point[1] * point[1] + point[2] * point[2];
	scalars->SetValue(iter.GetId(), value);
	iter.Next();
}
	
vtkNew<vtkMetaImageWriter> writer;
writer->SetInputData(imageData);
writer->SetFileName("data.mhd");
writer->Update();

Error stays exactly the same, even with your simplified function.

5

If you’re not tied to VTK, it’s pretty straightforward to do in SimpleITK, which allows mathematical operations on images.

Here’s how I compute your function using SimpleITK in Python:

# compute f(x, y, z) = 1 - 16xyz - 4x^2 - 4y^2 - 4z^2

import SimpleITK as sitk

diam = 200
radius = diam/2

# A 3d vector image where each voxel has its (x,y,z) location as its value
loc_img = sitk.PhysicalPointSource(sitk.sitkVectorFloat32, [diam,diam,diam])

# Split the vector image into X, Y and Z images
x_img = sitk.VectorIndexSelectionCast(loc_img, 0)
y_img = sitk.VectorIndexSelectionCast(loc_img, 1)
z_img = sitk.VectorIndexSelectionCast(loc_img, 2)

# subtract out the radius to make the range [-radius, radius]
x_img = x_img - radius
y_img = y_img - radius
z_img = z_img - radius

x2_img = x_img*x_img*(-4.0)
y2_img = y_img*y_img*(-4.0)
z2_img = z_img*z_img*(-4.0)

xyz_img = x_img*y_img*z_img*(-16.0)

result = xyz_img + x2_img + y2_img + z2_img + 1.0

sitk.WriteImage(result, 'cayley.mha')

I prefer this approach versus iterating over the pixels and computing the formula at each pixel.

SimpleITK also has C++ bindings, so the same code can be done in C++. I just find Python quicker and easier to develop in.

6

But how would I then write this with MetaImageWriter as an mhd-File? I wan’t to create the datasets so they can be later imported to a VTK program which runs them throug a VTK Pipeline.

7

I write out a MHA file, which is a MetaImage file. An MHA file is a MetaImage file with the pixel data in the file. MHD is a just a MetaImage header that points to a raw data file.

If you prefer MHD, just change the output file name to “Cayley.mhd”, and it will write that out.

8

I honestly didn’t notice the code was scrollable. I will try this in a second, thank you very much.

9

The problem might be writer->Update(). The VTK writers are designed to write the file with a call to writer->Write(). Some writers misbehave if you call Update() instead of Write().

10

The fixed the issue. Thank you!
I’ll still update later with the python solution.

11

Just in case you’re curious, the vtkImagePointIterator also works in Python, though with Python I’d almost recommend just doing the math in numpy, since numpy is faster than looping in Python. It depends on what your speed requirements are.

import vtk

# set image properties for an image with a domain of 
#  [ -10.0, +10.0] in x, y and z
spacing = [ 0.1, 0.1, 0.1 ]
size = [ 201, 201, 201 ]
origin = [ -10.0, -10.0, -10.0 ]

# create an image and allocate the pixel data
imageData = vtk.vtkImageData()
imageData.SetDimensions(size)
imageData.SetSpacing(spacing)
imageData.SetOrigin(origin)
# replace '1' with '2' if you want to create a complex image
imageData.AllocateScalars(vtk.VTK_FLOAT, 1)

# get the scalars from the image (they're allocated but not initialized).
scalars = imageData.GetPointData().GetScalars())

iter = vtk.vtkImagePointIterator(imageData)
while  not iter.IsAtEnd():
    x, y, z = iter.GetPosition()
    scalars.SetValue(iter.GetId(), x**2 + y**2 + z**2)
    iter.Next();
2 Likes