Having troubles saving some tags

[linako111]

I'm able to save images and most of the tags. But position and orientation are behaving not as I would expect. Btw, I made them work yesterday and then modified something and now cannot figure out the reason the code below does not work. When the code executes it prints correct values for updated tags "0020|0032" and "0020|0037". But when I open the metadata info in fiji it shows something different (perhaps some default values)

      for image in images:
            image = sitk.Cast(image, sitk.sitkInt16)
            direction = image.GetDirection()
            new_direction = '\\'.join(map(str, (direction[0], direction[2], direction[1], direction[3])))
            point = image.TransformIndexToPhysicalPoint((0,0))
            new_position = '\\'.join(map(str, (point[0], point[1])))

            for tag, value in series_tag_values:
                image.SetMetaData(tag, value)

            image.SetMetaData("0008|0012", time.strftime("%Y%m%d")) # Instance Creation Date
            image.SetMetaData("0008|0013", time.strftime("%H%M%S")) # Instance Creation Time
            image.SetMetaData("0008|0021", modification_date) # Series Date,
            image.SetMetaData("0008|0031", modification_time) # Series Time
            image.SetMetaData("0020|000e", seriesID) # Series Instance UID
            image.SetMetaData("0020|0013", f'{i+1:04d}') # Instance Number
            image.SetMetaData("0020|0032", new_position) # Image Position (Patient)                          
            image.SetMetaData("0020|0037", new_direction) # Image Orientation (Patient) 

            print("new position ", image.GetMetaData("0020|0032")) # this line prints correct info !!!
            print("new direction ", image.GetMetaData("0020|0037")) # this line prints correct info !!!

            writer.SetFileName(os.path.join(out_dir,'reco.' + f'{t:03d}' + '.' + f'{i*angle_step:03d}' + '.dcm'))
            writer.Execute(image)
            i += 1

[zivy]

Hello @linako111,

I am assuming that you are writing a 3D image as a set of 2D slices, new_direction should be a 6 element vector and new_position a 3 element one. While the slice is 2D it is actually [x,y,1] in 3D. If you look up the 0020|0032 and 0020|0037 you'll see that this is what they are expected to contain.

Two things that you will likely need to change:

1. image should be 3D.
2. Add the missing elements from direction when constructing new_direction and use (0,0,0) when computing point.

[linako111]

Hello @zivy,

I construct the image as MIP, using this code:

for angle in rotation_angles:
            rad = angle * np.pi / 180.0
            new_dir = (np.cos(rad), 0, np.sin(rad), 0, 1, 0, -np.sin(rad), 0, np.cos(rad))
            rotation_transform.SetRotation(rotation_axis, angle) 
            resampled_image = sitk.Resample(image1=image, 
                                            size=new_sz, 
                                            transform=rotation_transform, 
                                            interpolator=sitk.sitkLinear, 
                                            outputOrigin=min_bounds, 
                                            outputSpacing=new_spc, 
 #                                           outputDirection = [1,0,0,0,1,0,0,0,1])
                                           outputDirection=new_dir)
            proj_image = projection[ptype](resampled_image, paxis)
            extract_size = list(proj_image.GetSize())
            extract_size[paxis]=0
            extracted_image = sitk.Extract(proj_image, extract_size)   # this is now a 2D image!!!
            proj_images.append(extracted_image)  

The images become 2D. I want to save correct image orientation for each angle.

extracted_image.GetDirection() it gives me: (np.cos(rad), np.sin(rad), -np.sin(rad), np.cos(rad)) - 4 values only.

What is the correct way to fix it? Add a fake axis to the 2D image?

[zivy]

Hello @linako111,

Thanks for the clarification. Not really sure what is the correct thing to do with a MIP as it is not a volumetric image. I am not sure how DICOM encodes orientation for 2D images or if it is even part of the standard (my experience with cone-beam CT was that the image index in the cine loop matched the index in the set of calibration matrices which encoded the orientation).

Is there any reason to call the Extract method? Leaving proj_image as is seems to be the correct representation for the MIP, a 2D image(size_x, size_y, 1) in 3D space.

[linako111]

Thank you @zivy

I actually used your answer (as a template) in this thread: https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226/10

[zivy]

Hi @linako111,

The reason the images were extracted to 2D in that code was that I created a 3D faux volume there for visualization purposes, so didn't really care about the 3D aspect of things.

[dave3d]

SimpleITK 2d images have 2d orientation and direction. But when we read a DICOM slice, it creates a 3d image, with zdim=1.

Also I guess that you need to use the Image's SetOrigin and SetDirection methods, rather than trying to edit the meta data dictionary. I think when it writes out those particular entries into the DICOM, it pulls the information from the SimpleITK Image geometry parameters, not from the meta data dictionary.

[linako111]

Hi @zivy and @dave3d Thanks for your answers. I modified the code:

 for angle in rotation_angles:
            rad = angle * np.pi / 180.0
            new_dir = (np.cos(rad), 0, np.sin(rad), 0, 1, 0, -np.sin(rad), 0, np.cos(rad))
            rotation_transform.SetRotation(rotation_axis, angle) 
            resampled_image = sitk.Resample(image1=image, 
                                            size=new_sz, 
                                            transform=rotation_transform, 
                                            interpolator=sitk.sitkLinear, 
                                            outputOrigin=min_bounds, 
                                            outputSpacing=new_spc, 
#                                             outputDirection = [1,0,0,0,1,0,0,0,1])
                                            outputDirection=new_dir)
            proj_image = projection[ptype](resampled_image, paxis)
            extract_size = list(proj_image.GetSize())
            extract_size[paxis]=0
            extracted_image = sitk.Extract(proj_image, extract_size)
            slice_volume = sitk.JoinSeries(extracted_image)
            proj_images.append(slice_volume)
        processed_images.append(proj_images)

It works now. When I removed "Extract" part, because my paxis=1, the image's size became 290x1x277 which was wrong, that's why first I extracted 2D image and then added z direction. I hope it makes sense.

I also removed setting MetaData for these tags and you are right - the information still comes from the image.

Thank you very much for the help.