AllenNeuralDynamics / exa-spim-control

Acquisition control for the exaSPIM microscope system
MIT License
2 stars 1 forks source link

OpenCL flatfield #29

Closed adamkglaser closed 1 year ago

adamkglaser commented 1 year ago

@Poofjunior realized in my old files that I did take a first pass at trying to use OpenCL to do flat fielding on the GPU.

  1. I can't remember how the performance was.
  2. I'm not sure if the arithmetic I am doing with the kernel is optimal. Although it is just pixel wise division so I'm not sure how many options there would be (!)
import numpy as np
from gputools import OCLProgram, OCLArray, get_device

class FlatField():

    def __init__(self):

        # opencl kernel
        self.kernel="""
        __kernel void flatfield(__global short * input1, __global float * input2,
                                   __global short * output1){
          int i = get_global_id(0);
          int j = get_global_id(1);
          int Nx = get_global_size(0);
          output1[Nx*j+i] = input1[Nx*j+i]/input2[Nx*j+i];
        }
        """

        # create a dummy gaussian flatfield map, ff, to test with
        x, y = np.mgrid[-10640//2 + 1:10640//2 + 1, -14192//2 + 1:14192//2 + 1]
        ff = np.exp(-((y**2)/(2.0*5000**2)))
        ff = (ff/ff.max()).astype(np.float32)

        self.x_g2 = OCLArray.from_array(ff)
        self.y_g = OCLArray.empty(tuple((ff.shape[0], ff.shape[1])), np.uint16)

        self.prog = OCLProgram(src_str=self.kernel)

    def compute(self,image):

        x_g1 = OCLArray.from_array(image)
        self.prog.run_kernel(f'flatfield', self.y_g.shape[::-1], None, x_g1.data, self.x_g2.data, self.y_g.data)

        return self.y_g.get()
adamkglaser commented 1 year ago

Decision with Sci Comp to not do this on raw data and always make this part of the post-processing workflow.