difference between resample result in c++ version and python version

Hello all,
I am using ITK to resample CT 3D volume. I have done two version, one is c++, one is python. The difference between their result is not small. The bigger difference is existed in boundary part in the image.BTW, both of them used linear interpolator.
As I know, the python itk is wrapped on the c++ library. Am I right? So they should not have difference from in theory. I have also tried the SimpleITK in python. There is no difference between simpleITK and ITK results in python.
Here is code segments :

def resize_itk(data, new_shape, order=3):
    # Resample images to 2mm spacing with ITK
    import itk
    itk_image.SetSpacing( [0.791,0.791,1] )
    original_spacing = itk.spacing(itk_image)
    original_size =itk.size(itk_image)
    dimension = itk_image.GetImageDimension()
    out_spacing = [2.0,2.0,2.0]
    # out_size = [203,203,396]
    out_size = [
        int(np.round(original_size[dim] * (original_spacing[dim] / out_spacing[dim])+0.5)) for dim in range(dimension)

    transform = itk.IdentityTransform[itk.D, dimension].New()
    ImageType = type(itk_image)
    interpolator = itk.LinearInterpolateImageFunction[ImageType, itk.D].New()
    out_image = itk.resample_image_filter(
        # interpolator = interpolator,

    itk.imwrite(out_image, "/mnt/store/data/BR_DATASET/resampled-839316-py-itk-default-interpolator-itkF.nii.gz")
    return itk.GetArrayFromImage(out_image)
def resize_sitk(data, new_shape, order=3):
        # Resample images to 2mm spacing with SimpleITK
    import SimpleITK as sitk
    itk_image = sitk.GetImageFromArray(data)
    itk_image.SetSpacing( [0.791,0.791,1] )
    original_spacing = itk_image.GetSpacing()
    original_size = itk_image.GetSize()
    out_spacing = [2.0,2.0,2.0]
    # out_size = [203,203,396]
    out_size = [
        int(np.round(original_size[0] * (original_spacing[0] / out_spacing[0])+0.5)),
        int(np.round(original_size[1] * (original_spacing[1] / out_spacing[1])+0.5)),
        int(np.round(original_size[2] * (original_spacing[2] / out_spacing[2])+0.5))]

    resample = sitk.ResampleImageFilter()

    out_img = resample.Execute(itk_image)

    return sitk.GetArrayFromImage(out_img)

template< typename InputPixelType, typename OutputPixelType, typename InterpolatorPrecisionType = double,
		typename TransformPrecisionType = InterpolatorPrecisionType, uint16_t dim = 3 >
	static int Resampling(itk::SmartPointer<itk::Image<InputPixelType, dim>>& InputImg,
		double vfOutputSpacing[], const double vfOutputOrigin[],itk::Size<dim>& vnOutputSize, 
		itk::SmartPointer<itk::Image<OutputPixelType, dim>>& OutputImg,uint8_t nInterpolator=0)
		using InputImageType = itk::Image<InputPixelType, dim>;
		using OutputImageType = itk::Image<OutputPixelType, dim>;
		//auto inputRegion = InputImg->GetLargestPossibleRegion();
		//auto vnInputSize = inputRegion.GetSize();
		//auto oriSpacing = InputImg->GetSpacing();

		// Identity transform
		using T_Transform = itk::IdentityTransform<TransformPrecisionType, dim>;

		// The resampler type itself.
		using T_ResampleFilter = itk::ResampleImageFilter<InputImageType, OutputImageType,
			InterpolatorPrecisionType, TransformPrecisionType>;

		// Prepare the resampler.

		// Instantiate the transform and specify it should be the id transform.
		T_Transform::Pointer pTransform = T_Transform::New();

		// Instantiate the resampler. Wire in the transform and the interpolator.
		T_ResampleFilter::Pointer pResizeFilter = T_ResampleFilter::New();

		// If ITK resampler determines there is something to interpolate which is
		// usually the case when upscaling (!) then we must specify the interpolation algorithm.
		if (nInterpolator == LINEAR)
			using T_Interpolator = itk::LinearInterpolateImageFunction<InputImageType, InterpolatorPrecisionType>;
			// Instantiate the linear interpolator.
			T_Interpolator::Pointer pInterpolator = T_Interpolator::New();

		else if (nInterpolator == BSPLINE)
			//default order is 3
			using T_Interpolator = itk::BSplineInterpolateImageFunction<InputImageType, InterpolatorPrecisionType,
			T_Interpolator::Pointer pInterpolator = T_Interpolator::New();
				L"The interplator is unexpected.", std::to_wstring(UnexpInterplatorErr));
			return UnexpInterplatorErr;
		//if (nInterpolator == NEAREST)
		//	using T_Interpolator = itk::NearestNeighborInterpolateImageFunction<InputImageType, double>;
		//	T_Interpolator::Pointer pInterpolator = T_Interpolator::New();
		//	pResizeFilter->SetInterpolator(pInterpolator);

		// Set the output origin. You may shift the original image "inside" the
		// new image size by specifying something else than 0.0, 0.0 here.


		// Set the output spacing.

		// Set the output size.


		OutputImg = pResizeFilter->GetOutput();
		return Success;

Here is one slice in the difference result:

Here is the histogram of the difference result:

Is there any problems in the code? Can I get the same result in python with c++?
Looking forward to get your reply. Thanks a lot.

What is the meta-data of the output images?

Thanks @blowekamp . You gave me important hint.
I have checked the metadata from itksnap. and I found the pixdim is different in dim 1,2,3, the difference is 0.0014. I updated the target spacing in C++,but not in python. After adjusted them to the same value, currently the diff range is [-1.5368,1.68865]. This maybe caused by the different sform_code?
BTW, I encountered the other trouble. There is not metadata in the image which has been got from function GetImageFromArray. But if I write this image to nifti file, and read it to image, there is metadata in the image. It means the metadata is added in write function.
Here is the code segments:

import SimpleITK as sitk
import matplotlib.pyplot as plt
nii = r"C:\Users\xxx\Downloads\resampled-839316-py-sitk-2.0014.nii.gz"
img_itk = sitk.ReadImage(nii)
arr = sitk.GetArrayFromImage(img_itk)
print("read image:",img_itk.GetMetaDataKeys())
img_out = sitk.GetImageFromArray(arr)
print("get from array:",img_out.GetMetaDataKeys())
print("get from array:",img_out.GetSpacing())

output result:
read image: ('ITK_FileNotes', 'ITK_original_direction', 'ITK_original_spacing', 'aux_file', 'bitpix', 'cal_max', 'cal_min', 'datatype', 'descrip', 'dim[0]', 'dim[1]', 'dim[2]', 'dim[3]', 'dim[4]', 'dim[5]', 'dim[6]', 'dim[7]', 'dim_info', 'intent_code', 'intent_name', 'intent_p1', 'intent_p2', 'intent_p3', 'nifti_type', 'pixdim[0]', 'pixdim[1]', 'pixdim[2]', 'pixdim[3]', 'pixdim[4]', 'pixdim[5]', 'pixdim[6]', 'pixdim[7]', 'qform_code', 'qform_code_name', 'qoffset_x', 'qoffset_y', 'qoffset_z', 'quatern_b', 'quatern_c', 'quatern_d', 'scl_inter', 'scl_slope', 'sform_code', 'sform_code_name', 'slice_code', 'slice_duration', 'slice_end', 'slice_start', 'srow_x', 'srow_y', 'srow_z', 'toffset', 'vox_offset', 'xyzt_units')
get from array: ()
get from array: (1.0, 1.0, 1.0)
Backend TkAgg is interactive backend. Turning interactive mode on.

Here are my new questions:
How to set metadata before write it to file? Or I make some mistakes in the test above?
Any other hint?:slight_smile: Appreciated.

A new question:
I modify the sform_code in image with SetMetadata function, and write to file, and read to image, the sform_code is still 1.

import SimpleITK as sitk
import matplotlib.pyplot as plt
nii = r"C:\Users\xxx\Downloads\extremities\resampled-839316-py-sitk-2.0014.nii.gz"
img_itk = sitk.ReadImage(nii)
print("sform_code :",img_itk.GetMetaData('sform_code'))
arr = sitk.GetArrayFromImage(img_itk)
nii_out = r"C:\Users\xxx\Downloads\extremities\resampled-839316-py-sitk-2.0014-0.nii.gz"
img_sform_0 = sitk.ReadImage(nii_out)
print("sform_code :",img_sform_0.GetMetaData('sform_code'))

the output:
sform_code : 1
sform_code : 1
1 Like