Commit 4586367f authored by Pierre Paleo's avatar Pierre Paleo
Browse files

Start Cuda -log

parent 1b346a61
......@@ -60,3 +60,29 @@ __global__ void inplace_complexreal_mul_2Dby2D(complex* arr2D_out, float* arr2D_
arr2D_out[i]._M_re *= b;
arr2D_out[i]._M_im *= b;
}
#ifndef DO_CLIP_MIN
#define DO_CLIP_MIN 0
#endif
#ifndef DO_CLIP_MAX
#define DO_CLIP_MAX 0
#endif
// arr = -log(arr)
__global__ void nlog(float* array, float* output, int Nx, int Ny, int Nz, float clip_min, float clip_max) {
int x = blockDim.x * blockIdx.x + threadIdx.x;
int y = blockDim.y * blockIdx.y + threadIdx.y;
int z = blockDim.z * blockIdx.z + threadIdx.z;
if ((x >= Nx) || (y >= Ny) || (z >= Nz)) return;
int pos = (z*Ny + y)*Nx + x;
float val = array[pos];
#if DO_CLIP_MIN
val = fmaxf(val, clip_min);
#endif
#if DO_CLIP_MAX
val = fminf(val, clip_max);
#endif
output[pos] = -logf(val);
}
......@@ -49,6 +49,38 @@ class CCDProcessing(object):
self.shape = (n_z, n_x)
def take_logarithm(self, radios=None, output=None, clip_min=None, clip_max=None):
"""
Take the negative logarithm of a radios chunk.
Parameters
-----------
radios: numpy.ndarray, optional
If not provided, the radios chunk provided at class instantiation is
taken.
output: numpy.ndarray, optional
Output of -log(radios).
If not provided, a new copy is created.
clip_min: float, optional
Before taking the logarithm, the values are clipped to this minimum.
clip_max: float, optional
Before taking the logarithm, the values are clipped to this maximum.
"""
if radios is None:
radios = self.radios
if output is None:
output = np.zeros_like(radios)
if (clip_min is not None) or (clip_max is not None):
np.clip(radios, clip_min, clip_max, out=output)
np.log(output, out=output)
else:
np.log(radios, out=output)
output[:] *= -1;
return output
# Possible inputs scenarios:
# (1) radios: array, flats/darks: dict of DataUrl
# In this case, we should keep the mapping between range(radios.shape[0])
......
......@@ -13,6 +13,9 @@ except ImportError:
class CudaFlatField(FlatField):
def __init__(
self,
......@@ -113,8 +116,6 @@ class CudaFlatField(FlatField):
return radios
class CudaCCDCorrection(CCDCorrection):
def __init__(
self,
......@@ -182,3 +183,70 @@ class CudaCCDCorrection(CCDCorrection):
else:
radios = self.radios
return self.cuda_median_filter.medfilt2(radios, output=output)
class CudaLog(CCDProcessing):
"""
Helper class to take -log(radios)
"""
def __init__(self, radios, clip_min=None, clip_max=None):
super().__init__(radios)
self.clip_min = clip_min
self.clip_max = clip_max
self._init_kernels()
def _init_kernels(self):
self._do_clip_min = int(self.clip_min is not None)
self._do_clip_max = int(self.clip_max is not None)
self._nlog_srcfile = get_cuda_srcfile("ElementOp.cu")
nz, ny, nx = radios.shape
self._nx = np.int32(nx)
self._ny = np.int32(ny)
self._nz = np.int32(nz)
self.nlog_kernel = CudaKernel(
"nlog",
filename=self._nlog_srcfile,
signature="PPiiiff",
options=[
"-DDO_CLIP_MIN=%d" % self._do_clip_min,
"-DDO_CLIP_MAX=%d" % self._do_clip_max,
]
)
def take_logarithm(self, radios=None, output=None, clip_min=None, clip_max=None):
"""
Take the negative logarithm of a radios chunk.
Parameters
-----------
radios: numpy.ndarray, optional
If not provided, the radios chunk provided at class instantiation is
taken.
output: pycuda.gpuarray.GPUArray, optional
Output of -log(radios).
If not provided, a new GPU array is created.
clip_min: float, optional
Before taking the logarithm, the values are clipped to this minimum.
clip_max: float, optional
Before taking the logarithm, the values are clipped to this maximum.
"""
if radios is None:
radios = self.radios
if output is None:
output = garray.zeros_like(radios)
self.nlog_kernel(
radios,
output,
self._nx,
self._ny,
self._nz,
self.clip_min,
self.clip_max
)
return output
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment