Pixel manipulations in C#

Beginner Questions
1

Hello,

I am trying to implement the following example from SimpleITK’s notebook.

The Python code is:

img = logo

# Generate random samples inside the image, we will obtain the intensity/color values at these points.
num_samples = 10
physical_points = []
for pnt in zip(*[list(np.random.random(num_samples)*sz) for sz in img.GetSize()]):
    physical_points.append(img.TransformContinuousIndexToPhysicalPoint(pnt))

# Create an image of size [num_samples,1...1], actual size is dependent on the image dimensionality. The pixel
# type is irrelevant, as the image is just defining the interpolation grid (sitkUInt8 has minimal memory footprint).
interp_grid_img = sitk.Image([num_samples] + [1]*(img.GetDimension()-1), sitk.sitkUInt8)

# Define the displacement field transformation, maps the points in the interp_grid_img to the points in the actual
# image.
displacement_img = sitk.Image([num_samples] + [1]*(img.GetDimension()-1), sitk.sitkVectorFloat64, img.GetDimension())
for i, pnt in enumerate(physical_points):
     displacement_img[[i] + [0]*(img.GetDimension()-1)] = np.array(pnt) - np.array(interp_grid_img.TransformIndexToPhysicalPoint([i] + [0]*(img.GetDimension()-1)))

# Actually perform the resampling. The only relevant choice here is the interpolator. The default_output_pixel_value
# is set to 0.0, but the resampling should never use it because we expect all points to be inside the image and this
# value is only used if the point is outside the image extent.
interpolator_enum = sitk.sitkLinear
default_output_pixel_value = 0.0
output_pixel_type = sitk.sitkFloat32 if img.GetNumberOfComponentsPerPixel()==1 else sitk.sitkVectorFloat32
resampled_points = sitk.Resample(img, interp_grid_img, sitk.DisplacementFieldTransform(displacement_img), 
                                 interpolator_enum, default_output_pixel_value, output_pixel_type)

# Print the interpolated values per point
for i in range(resampled_points.GetWidth()):
      print(str(physical_points[i]) + ': ' + str(resampled_points[[i] + [0]*(img.GetDimension()-1)]) + '\n')

The problem is that in C#, neither indexing nor SetPixel work for setting the displacement_img. I have seen examples of buffer access, but only with Int8, not for vector images. How might you do it?

Also, here interp_grid_img is the ReferenceImage and sitk.DisplacementFieldTransform(displacement_img) is the Transform in the object-oriented interface, right?

Thank you!

2

Can you please provide a minimal example in C# to reproduce the problem encountered, and the error message received?

3

Of course, here is the code giving me trouble:

uint[] gridSize = new uint[] { (uint)Math.Floor(insonationWidth * sampleFactor), (uint)Math.Floor(insonationHeight * sampleFactor), 1 };
        Image grid = new Image(new VectorUInt32(gridSize), itk.simple.PixelIDValueEnum.sitkUInt8);
        grid.SetSpacing(new VectorDouble() { pixel_size / sampleFactor, pixel_size / sampleFactor, pixel_size / sampleFactor });
        Image displacement_field = new Image(new VectorUInt32(gridSize), itk.simple.PixelIDValueEnum.sitkVectorFloat64, 3);
        for (int i = 0; i < insonationWidth * sampleFactor; i++)
        {
            for (int j = 0; j < insonationHeight * sampleFactor; j++)
            {
                VectorDouble gridPoint = grid.TransformIndexToPhysicalPoint(new VectorInt64() { i, j, 1 });
                displacement_field[i, j, 1] = physicalSamplePoints[(int)Math.Floor(j + i * insonationHeight * sampleFactor)] - new Vector3((float)gridPoint[0],(float)gridPoint[1],(float)gridPoint[2]);
            }
        }

The error message is “Cannot apply indexing with [] to an expression of type ‘Image’” (when I try to set displacement_field). If I try SetPixel, it tells me the class Image contains no such method.

4

Thank you for the code.

Since C# is a typed language you need to use SetPixelAsVectorUInt32 method.

I’ll need to look into the [] operator further…

1 Like
5

Thank you! That worked.

Also, this code is very slow. Is that to be expected? As I said, I’m trying to resample at a given set of points.

6

Hello @Ale_Ballester,

If the goal is to resample at a set of points, the next version of SimpleITK will have new methods for the Image class EvaluateAtPhysicalPoint and EvaluateAtContinuousIndex (get the compiled release candidate from github). A Python notebook illustrating the use is here relevant section Resampling at a set of locations.

7

I am trying to construct a 2D image from an arbitrary slice of a 3D image. Since I have the 4x4 transformation matrix that can give me arbitrary points in my plane of interest, I thought resampling at a particular set of sample points within the plane and setting a default pixel value was the easiest way to do it. It was easy, but way too slow (I am writing real time software).

8

I think it is possible to do this operation fast enough with numpy arrays. I think you would need to operate in the indices preserving the scalar values of the CT, then you should do a trivial threshold operation. The transform should be something like SliceToRAS on Slicer or it’s inverse

9

Could you use the Resample function and to resample the 3d image onto your 2d target image?

1 Like
10

Yes, but I need to implement it on C#/Unity for a standalone application (no slicer).

11

I am trying to do it using Resample with identity transform, setting the origin, direction and output size, and then use ExtractImageFilter to get my slice. So far no success.

1 Like
12

Keep in mind I am generating 192,000 points per frame. Do you really think this can be made fast enough?

13

Yeah

(If I understood your question well)

14

I made a simple example (in Python) on how to do an oblique 2d cut through a 3d volume using SimpleITK’s Resample function. Here’s how I do it:

import SimpleITK as sitk

# Example volume: https://raw.githubusercontent.com/dave3d/dicom2stl/main/examples/Data/ct_example.nii.gz
img = sitk.ReadImage("ct_example.nii.gz")

cut = sitk.Image([300,300,1], sitk.sitkUInt16)

rotate = sitk.Euler3DTransform()

# 15 degree rotation about X axis
rotate.SetRotation(.2618, 0.0, 0.0)

cut.SetDirection(rotate.GetMatrix())

cut2 = sitk.Resample(img, cut)

I made a Jupyter notebook that uses ITKWidgets to show the result, if you want to try it yourself:
Cut-resample.ipynb · GitHub

3 Likes
15

Thank you!
I ended up solving it using resample with an identity transform. But this is great!

3 Likes