An efficient way to rotate a two-dimensional plane and add it to a three-dimensional volume

General Usage
,
1

I am currently working on transforming a two-dimensional array into a three-dimensional array by means of Euler angles (ψ, θ, φ) and adding it to another three-dimensional array. (like this Single particle analysis - Wikipedia)

The immediate method that comes to mind is to add a 2D array to the centre of the Z-axis of a 3D array initialised with zeros and rotate it using ImageTransformations warp. This method is good because you can choose the interpolation method, but it is inefficient when large numbers of zeros are generated and added to another 3D array, and it is slow because you cannot use the GPU’s texture interpolation.(maybe)

Is there any good way to do this?

2

and adding it to another three-dimensional array

Wouldn’t something like this work? (Schematic of a solution)

for I in CartesianIndices(img2d)
    J = map2dcoordto3d(I)    # compute the position of this pixel in img3d
    img3d[J] += img2d[I]
end
3

Thanks for your reply.
Yes, I definitely want to do that, but I am struggling to find an efficient way to do it considering right interpolations.
I will look at the source code of other languages for a while. (Although it’s very painful as I only know Julia!)

4

It should be fairly straightforward:

  • construct your transformation with CoordinateTransformations.jl
  • pad your 2d coord with 0 (for the z coordinate)
  • apply the transformation to the padded coordinate
  • round to the nearest pixel or use interpolation
5

It works,thanks.
However, I want to speed up processing on VRAM as much as possible.
I have simplified the problem.
I want to construct a sparse linear map that rotates a vector of N^2 (2D image) by θ radians to a vector of N^2. I’ll think about it for a while, and if I can do this, I can speed up the process considerably!

6

Final answer.

using CUDA, TestImages, ImageShow
image = testimage("cameraman") 
arr = image .|> Float16

gpu_src = CuTextureArray(CuArray(arr))
gpu_dst = CuArray{Float16}(undef, size(arr))

@inline function calcpoint(blockIdx, blockDim, threadIdx, size)
    i = (blockIdx - 1) * blockDim + threadIdx
    return i, Float32(i)/size
end

function kernel_texture_warp_native(dst::CuDeviceArray{<:Any,2}, texture::CuDeviceTexture{<:Any,2},theta)
    i, u = calcpoint(blockIdx().x, blockDim().x, threadIdx().x, size(dst)[1])
    j, v = calcpoint(blockIdx().y, blockDim().y, threadIdx().y, size(dst)[2])

    # transform coordinates
    u -= 0.5
    v -= 0.5
	tu = u*CUDA.cos(theta) - v*CUDA.sin(theta) + 0.5;
    tv = v*CUDA.cos(theta) + u*CUDA.sin(theta) + 0.5;
	
    dst[i,j] = texture[tu,tv]
    return nothing
end

begin
	kernel = @cuda launch=false kernel_texture_warp_native(gpu_dst, gpu_tex,pi)
	config = launch_configuration(kernel.fun)
	
	dim_x, dim_y= size(gpu_tex, 1), size(gpu_tex, 2)
	threads_x = min(dim_x, config.threads)
	blocks_x = cld(dim_x, threads_x)
	
	kernel(gpu_dst, gpu_tex,pi; threads=threads_x, blocks=(blocks_x, dim_y))
end

I’ll try rotate 2d to 3d tomorrow.