Trying to load custom loaded images into vtkImageImport

Hello everyone,
I’m trying to load custom loaded images (because at my workplace i must use this library) to vtkImageImport (and at the end to vtkImageData, to render a Volume).

I have this object dicomStack that already contain all the images loaded from a dicom folder, already ordered. So I am pretty sure we can assume the images and their pixels are good and ready to be used.

I have already seen topic as: this one but that doesn’t helped me much.

This is the current “pseudocode” that I am using, and that doesn’t work:

//dicomStack is the object with all the images from the dicom folder loaded
//for every image append the raw data (for now is single frame)
QVector<void*> imageRawData;
for(int i=0; i < dicomStack.size(); i++){
    DcmImage currentImage = dicomStack.image(i);
    if(currentImage.numFrames()>0){
        imageRawData.append(currentImage.getRawPixelData(0));
    }
}
//Create imageImport
vtkSmartPointer<vtkImageImport> imageImport =
    vtkSmartPointer<vtkImageImport>::New();
imageImport->SetDataSpacing(spacing.x(), spacing.y(), spacing.z());
imageImport->SetDataOrigin(origin.x(), origin.y(), origin.z());
imageImport->SetWholeExtent(0, imageSize.width(), 0, imageSize.height(), 0, dicomStack.size() - 1);
imageImport->SetDataExtentToWholeExtent();
imageImport->SetDataScalarTypeToUnsignedChar();
imageImport->SetNumberOfScalarComponents(1);
imageImport->SetImportVoidPointer(imageRawData.data());
imageImport->Update();
[...]
m_imageData->ShallowCopy(imageImport->GetOutput());

The problem is… I can’t understand why it doesn’t work. The import doesn’t cause any error, when i call renderWindow()->Render() the execution get stuck here: (Common\Core\vtkDataArrayPrivate.txx)


And i can’t see why. I tried to understand where that error comes from but i had not results.

Any suggestion? What I could be missing?

PS: That dicom folder is valid, i can succesfully render a volume loading it with vtkDICOMImageReader

As someone may notice, having QVector<void*> might be wrong because it’s a list of pointers and SetImportVoidPointer might not like it.
But i tried using QVector<short> / short* = new short[size] too and always the same error appear

(I converted the data using *(reinterpret_cast<short*>(currentImage.getRawPixelData(0)));)

Hi Bandark, I definitely sympathize with being forced to use a tool that might not be the best tool for the job.

So, no, a vector of pointers will definitely not work with vtkImageImport. The SetImportVoidPointer() method requires a single pointer to a contiguous array that contains all of the images. Therefore you will need to create vector<unsigned char>, vector<short>, or vector<unsigned short>. I’ll assume that you know which.

To this array, you will need to append the entire contents of each of your images (i.e. all the pixels). This could be done by calling insert() repeatedly, but probably the most efficient way would be to resize() the vector to make it large enough to hold all the pixels of all the images, and then use memcpy() calls to copy the pixels from the images into the array.

I have to ask, what type does getRawPixelData(x) return? And what does its argument do? And how do you get the size of the data that it returns?

First of all, thank you for your time and your reply ^^

  • The folder I am using to test this import contain 16bit images, so either short or unsigned short. I can check with image.bitsStored() that in this case is 16. Also image.pixelRepresentation() returns int16 (so probably short because the flag UInt16 exist for unsigned).
  • The argument of getRawPixelData(x) is the frame, this because there might be a case where all the frames of the volume are in a single “image”. In this case every image has a single frame so it’s always 0.
  • getRawPixelData(x) returns a buffer void* that points to the raw pixel data of the frame
  • About the size, I am not 100% sure, the best way seems doing image.pixelData(0)->size(). This returns 524’288 (512x512x2), and I am being told that is *2 because is 16bit.
    Why? How should i use this size?

From what you’ve told me, I’m going to use these assumptions:

  • image.pixelData(0)->size() returns the buffer size (in bytes) of the first frame
  • image.getRawPixelData(0) returns the pointer to the buffer itself

What VTK needs is a buffer for the entire image volume. To keep things simple, let’s assume that each image has just one frame. So the buffer size that VTK needs is the buffer size per image frame, multiplied by the total number of images.

So you’re going to need a huge buffer, large enough to store all the images. At this point, let’s forget about vtkImageImport. The most efficient way to approach this is to directly access the vtkImageData’s buffer, rather than creating a temporary buffer. Then vtkImageImport becomes completely unnecessary.

// VTK will allocate a buffer of the desired size and type
int dataType = VTK_SHORT; // make sure this is correct!
int dataTypeSize = 2; // each 'short' is 2 bytes
int numberOfComponents = 1;
m_imageData->SetDimensions(imageSize.width(), imageSize.height(), dicomStack.size());
m_imageData->AllocateScalars(dataType, numberOfComponents); // allocate buffer
m_imageData->SetSpacing(spacing.x(), spacing.y(), spacing.z());

// compute the number of bytes required to store each slice of the volume
size_t bytesPerSlice = dataTypeSize * numberOfComponents;
bytesPerSlice *= imageSize.width();
bytesPerSlice *= imageSize.height();

for (int i=0; i < dicomStack.size(); i++) {
    // are you sure that you don't want to use a reference here?
    DcmImage currentImage = dicomStack.image(i);
    if (currentImage.numFrames() > 0) {
        size_t bytesPerFrame = image.pixelData(0)->size();
        // do a sanity check on the size
        if (bytesPerFrame == bytesPerSlice) {
            // copy the image into the correct slice of the vtkImageData buffer 
            memcpy(m_imageData->GetScalarPointer(0, 0, i), currentImage.getRawPixelData(0), bytesPerSlice);
        }
    }
}

This works because m_imageData->GetScalarPointer(k, j, i) returns a void* to the correct position within the vtkImageData’s buffer (where k=column, j=row, i=slice).

Please note that I just typed this code and haven’t compiled it, so it probably contains typos.

1 Like

Thank you so much for your suggestions, I will try this version in the next days and I will see what i can get.

In the worst case i’ll read the pixels one by one to create the correct buffer.

I will update you as soon as i will be able to make something work.

First of all, thank you so much again.
I couldn’t resist but try your solution as soon as possible.
The result is interesting, it actually work, for the most part. The result is this (on the right):

As you can see, most of the data kinda seems to be there but looks like there’s something wrong.

Any idea of what could cause this problem, or is something about spacing/order i should investigate on my end?

The spacing looks correct, but the order looks wrong. What guarantees do you have that the dicomStack is in the correct order? Take a look at the InstanceNumber of each image, or even better look at the ImagePositionPatient of each image.

The other thing to look at is the RescaleIntercept and RescaleSlope of each image in the stack. The vtkDICOMImageReader will use these to rescale the raw data (multiply by slope, and add intercept). So if you directly copy the raw DICOM data into a vtkImageData, without rescaling the data values, then the values you get will not be the same as vtkDICOMImageReader (i.e. you will get the raw, uncalibrated data values).

1 Like

Can confirm the order was wrong, reordering the images result in a more “correct” result:
image

Will investigate more stuff stuch as RescaleIntercept and RescaleSlope in the next days to see if i can understand why the result is still different. (Every other suggestion is welcome if you have any)

Thank you again and have a nice weekend

Edit: Changing the transfer function to visualize more inside in that model, all the bones and the data is correctly there, but why it need a different preset? I’d like to utilize the presets i already have, so maybe i have to edit the data as you said(?)

Let’s assume I get RescaleIntercept and RescaleSlope data.
There’s any VTK filter/way to apply them or should I just find a way to apply them to raw data?
Edit1: (Like using vtkImageShiftScale?)
Edit2: Made it work using vtkImageShiftScale, that’s all for now, thank you so much

Yes, you can use vtkImageShiftScale, but be careful because vtkImageShiftScale applies the equation y = (x + Shift)*Scale, but DICOM uses the equation y = x*RescaleSlope + RescaleIntercept. Usually you won’t see a difference because the RescaleSlope will almost always be 1 for Hounsfield units, but it’s still important to get the math right.

Also, be sure to use vtkImageShiftScale::ClampOverflowOn() to ensure that the rescaling doesn’t overflow the data type.


I’ll add a few caveats here.

  • DICOM defines the CT RescaleSlope and RescaleIntercept per image, not per stack, and in some rare cases the values won’t be the same for all images in the stack. So if you want your code to be 100%, you would need to do the rescaling ‘per image’ when you create the vtkImageData, instead of doing it afterwards with vtkImageShiftScale.
  • In general the RescaleSlope and RescaleIntercept should not be used for MR images.
  • To be absolutely correct, the DICOM standard requires that you mask the raw data values so that all bits except the BitsStored are zero. Then, if PixelRepresentation is 1, it is necessary to sign extend the value. This should be done before applying any rescaling.
1 Like

Thank you so much for this explanation

How would you implement point 3?
Because i have a problem in importing 12bit images and i think it could be related to that somehow

The bit-masking code that I use in my own DICOM reader is located here and here on github.

That looks very helpful but apparently is not the problem (or i am calling it wrong).

My problem is that a 12bit TC with some transfer functions looks like this:

Instead of this:

Any idea where the error comes from?

After you apply the RescaleSlope and RescaleIntercept, make sure that the output data type is signed, because CT values are signed. Usually they range from -1000 to around +3000.

If you do not do this, but instead use unsigned values and clamp all negative values to zero, then the result will look like your first image above.

If you are using vtkImageShiftScale, use SetOutputScalarTypeToShort() or even SetOutputScalarTypeToFloat().

1 Like

This solved it, thank you so much again