Now suppose that in the original image, there is a voxel at point (x,y,z). I now want to calculate/find the coordinates of this point in the new image space.

In the code you provided, the transformation is the identity so the index of a point changes in the resampled_image, but its physical coordinates do not. The various options of moving between physical coordinates and indexes are illustrated below:

original_index = [x, y, z]
original_physical_coords = image.TransformIndexToPhysicalPoint(original_index)
# If you have a continuous value for the index, sub pixel use TransformContinuousIndexToPhysicalPoint
# Obtain the corresponding index in the resampled image, either continuous or discrete
new_index_continous = resampled_image.TransformPhysicalPointToContinuousIndex(original_physical_coords)
new_index = resampled_image.TransformPhysicalPointToIndex(original_physical_coords)

Great thanks! That is exactly what I was looking for. Also thanks for pointing out the difference between physicial space (coordinates) and image space (indices), this cleared it up for me. Thanks again