Pixel manipulations in C#


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()]):

# 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!

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

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.

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

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.

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.

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).

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

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

1 Like

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

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

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


(If I understood your question well)

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)


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


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