Non isovolumetric voxels in case of rotations may lead to issues

[jakubMitura14]

Becouse rotations.jl seem to have no support of rotating the image for the case when non isovolumetric voxels are present in some cases the resulting image can lose correct representation in real world coordinates. To understand it let's give example let's say we have very "high voxels" so for example voxel has 5 mm in z axis and only 1 mm in x and y axis. We know that the more distant from the center of rotation is the point the bigger the arc it need to travel. Hence in our example the point that will be 10 voxels in z direction should move in far greater arc than the voxel that is 10 voxels from the center in x direction. Hence, we need to work in world coordinate system and take it into account.

On the image on the left it is how rotation would look like if we would ignore spacing in green distance from point to center of rotation; and in yellow arc that a point moves during rotation; on the right when we take the spacing into account. Image is a rough schamtic lenths and proportions are only approximate.

Proposed solution can look somewhat like this:

1. manually using sine and cosines or using rotations.jl obtain rotation matrix
2. get cartesian indices of all points of the voxel array we want to rotate using "get_base_indicies_arr" function from Utils.jl
3. transform cartesian indicies (let's name them points_to_interpolate) into real world coordinates by taking into account spacing and origin sth like below

    points_to_interpolate = get_base_indicies_arr(new_size)
    points_to_interpolate=points_to_interpolate.*new_spacing
     points_to_interpolate = points_to_interpolate.+origin

4. using tensor product apply rotation matrix on coordinates of each point from points_to_interpolate
5. now we have new coordinates of each point we should test it against Rotations.jl (in case of spacing equal (1,1,1) - so isovolumetric results should be identical)
6. we should concatenate the value of each point with its coordinate so we will have vector of length for each voxel in original image
7. now we need to do diffrently depending on weather we want to have the same shape/size of the array as input image (default case) or diffrent. We need to consider a different bigger shape as it can be necessary in order to fit all new points.
   a) in case we want to have the same shape as input image we ignore all points that get outside of physical space of the original image (we can establish it in each axis by mutliplying number of voxels in each axis times thie spacing in this axis and taking into account origin; so for example given origin 10,11,13 and spacing 0.5,0.6,0.8 and image size 90,100,110 the y axis will start at 11 and end at 71).
   b) in case we do not want to have the same shape as original image we need to increase the size of the resulting array and if needed move the origin (so for example as above if the lowest value of y in all points will be 9 we need to set the origin value of y to 9 so 10,9,13 and increase the array size by 4 {11-9=2; ceil(2/0.6)=4 } , analogously for all axes, in case of points getting out down the axis we also need to increase the size of the resulting image but we do not need to change the position of the origin)
8. having now empty result array; points values and their coordinates we need to interpolate them with the grid associated with output image. In default case where origin position do not change we can reuse the points_to_interpolate we had before applying rotation matrix; for case we want to change result array shape we can get output_points in analogous way as we obtained points_to_interpolate
9. We need now to interpolate the values we have in points_to_interpolate onto the output_points via interpolator chosen by the user.
   generally one should adapt my_interpolate from Utils.jl - add necessery arguments in it if any modification of this function is required look into griddeed interpolation Gridded as it will probably show what modifications are needed . In case it would not work We can have 3 strategies how to get this problem solved - 1) check how it is performed in Rotations.jl; 2) use Interpolations.jl via griddeed interpolation Gridded if possible 3) use ScatteredInterpolation.jl
10. test against the simple itk in the test that are present in the test file; visualize the results of simple itk and medimage rotation as nifti files and compare them manually in 3dSlicer

Notes:
We need also to establish is image orientation affecting those calculations, test on diffrent orientations that can be set by change_orientation function from spatial_metadata_change file. in spatial_metadata_change test file there is simple itk function also for testing if needed.

As the inspiration look into resample_to_spacing function implementation

Nearly identical procedure but instead of applying the rotation matrix we would add the translation vector can be used in a variant of translation (other variant is just change the origin; both variants are usefull)