Source code for ClearMap.ImageProcessing.Filter.FilterKernel

# -*- coding: utf-8 -*-
"""
FilterKernel
=============

Implementation of various volumetric filter kernels.


.. _FilterTypes:

Filter Type
-----------

Filter types defined by the ``ftype`` key include: 

=============== =====================================
Type            Descrition
=============== =====================================
``mean``        uniform averaging filter
``gaussian``    Gaussian filter
``log``         Laplacian of Gaussian filter (LoG)
``dog``         Difference of Gaussians filter (DoG)
``sphere``      Sphere filter
``disk``        Disk filter
=============== =====================================
"""
__author__    = 'Christoph Kirst <christoph.kirst.ck@gmail.com>'
__license__   = 'GPLv3 - GNU General Pulic License v3 (see LICENSE.txt)'
__copyright__ = 'Copyright © 2020 by Christoph Kirst'
__webpage__   = 'http://idisco.info'
__download__  = 'http://www.github.com/ChristophKirst/ClearMap2'

import numpy as np
import math

import ClearMap.ImageProcessing.Filter.StructureElement as se


###############################################################################
### Filter kernel
###############################################################################

[docs] def filter_kernel(ftype = 'Gaussian', shape = (5,5), radius = None, sigma = None, sigma2 = None): """Creates a filter kernel of a special type Arguments --------- ftype : str Filter type, see :ref:`FilterTypes`. shape : array or tuple Shape of the filter kernel. radius : tuple or float Radius of the kernel (if applicable). sigma : tuple or float Std for the first gaussian (if applicable). sigma2 : tuple or float Std of a second gaussian (if present). Returns ------- kernel : array Filter kernel. """ ndim = len(shape); if ndim == 2: return filter_kernel_2d(ftype=ftype, shape=shape, sigma=sigma, sigma2=sigma2, radius=radius); elif ndim == 3: return filter_kernel_3d(ftype=ftype, shape=shape, sigma=sigma, sigma2=sigma2, radius=radius); else: raise ValueError('Kernel for %d dimensions not implemented !' % ndim);
[docs] def filter_kernel_2d(ftype = 'Gaussian', shape = (5,5), sigma = None, sigma2 = None, radius = None): """Creates a 2d filter kernel of a special type Arguments --------- ftype : str Filter type, see :ref:`FilterTypes`. shape : array or tuple Shape of the filter kernel. radius : tuple or float Radius of the kernel (if applicable). sigma : tuple or float Std for the first gaussian (if applicable). sigma2 : tuple or float Std of a second gaussian (if present). Returns ------- kernel : array Filter kernel. """ ftype = ftype.lower(); o = se.structure_element_offsets(shape); mo = o.min(axis=1); shape = np.array(shape); if ftype == 'mean': # unifrom return np.ones(shape)/ shape.prod(); elif ftype == 'gaussian': if sigma == None: sigma = shape / 2. / math.sqrt(2 * math.log(2)); sigma = np.array(sigma); if len(sigma) < 3: sigma = np.array((sigma[0], sigma[0])); else: sigma = sigma[0:2]; g = np.mgrid[-o[0,0]:o[0,1], -o[1,0]:o[1,1]]; add = ((shape + 1) % 2) / 2.; x = g[0,:,:,:] + add[0]; y = g[1,:,:,:] + add[1]; ker = np.exp(-(x * x / 2. / (sigma[0] * sigma[0]) + y * y / 2. / (sigma[1] * sigma[1]))); return ker/ker.sum(); elif ftype == 'sphere': if radius == None: radius = mo; radius = np.array(radius); if len(radius) < 3: radius = np.array((radius[0], radius[0])); else: radius = radius[0:2]; g = np.mgrid[-o[0,0]:o[0,1], -o[1,0]:o[1,1]]; add = ((shape + 1) % 2) / 2.; x = g[0,:,:,:] + add[0]; y = g[1,:,:,:] + add[1]; ker = 1 - (x * x / 2. / (radius[0] * radius[0]) + y * y / 2. / (radius[1] * radius[1])); ker[ker < 0] = 0.; return ker / ker.sum(); elif ftype == 'disk': if radius == None: radius = mo; radius = np.array(radius); if len(radius) < 3: radius = np.array((radius[0], radius[0])); else: radius = radius[0:2]; g = np.mgrid[-o[0,0]:o[0,1], -o[1,0]:o[1,1]]; add = ((shape + 1) % 2) / 2.; x = g[0,:,:,:] + add[0]; y = g[1,:,:,:] + add[1]; ker = 1 - (x * x / 2. / (radius[0] * radius[0]) + y * y / 2. / (radius[1] * radius[1])); ker[ker < 0] = 0.; ker[ker > 0] = 1.0; return ker / ker.sum(); elif ftype == 'log': # laplacian of gaussians if sigma == None: sigma = shape / 4. / math.sqrt(2 * math.log(2)); sigma = np.array(sigma); if len(sigma) < 3: sigma = np.array((sigma[0], sigma[0])); else: sigma = sigma[0:2]; g = np.mgrid[-o[0,0]:o[0,1], -o[1,0]:o[1,1]]; add = ((shape + 1) % 2) / 2.; x = g[0,:,:,:] + add[0]; y = g[1,:,:,:] + add[1]; ker = np.exp(-(x * x / 2. / (radius[0] * radius[0]) + y * y / 2. / (radius[1] * radius[1]))); ker /= ker.sum(); arg = x * x / math.pow(sigma[0], 4) + y * y/ math.pow(sigma[1],4) - (1/(sigma[0] * sigma[0]) + 1/(sigma[1] * sigma[1])); ker = ker * arg; return ker - ker.sum()/len(ker); elif ftype == 'dog': if sigma2 == None: sigma2 = shape / 2. / math.sqrt(2 * math.log(2)); sigma2 = np.array(sigma2); if len(sigma2) < 3: sigma2 = np.array((sigma2[0], sigma2[0])); else: sigma2 = sigma2[0:2]; if sigma == None: sigma = sigma2 / 1.5; sigma = np.array(sigma); if len(sigma) < 3: sigma = np.array((sigma[0], sigma[0])); else: sigma = sigma[0:2]; g = np.mgrid[-o[0,0]:o[0,1], -o[1,0]:o[1,1]]; add = ((shape + 1) % 2) / 2.; x = g[0,:,:,:] + add[0]; y = g[1,:,:,:] + add[1]; ker = np.exp(-(x * x / 2. / (sigma[0] * sigma[0]) + y * y / 2. / (sigma[1] * sigma[1]))); ker /= ker.sum(); sub = np.exp(-(x * x / 2. / (sigma2[0] * sigma2[0]) + y * y / 2. / (sigma2[1] * sigma2[1]))); return ker - sub / sub.sum(); else: raise ValueError('Filter type %r not valid!' % ftype);
[docs] def filter_kernel_3d(ftype = 'Gaussian', shape = (5,5,5), sigma = None, sigma2 = None, radius = None): """Creates a 3d filter kernel of a special type Arguments --------- ftype : str Filter type, see :ref:`FilterTypes`. shape : array or tuple Shape of the filter kernel. radius : tuple or float Radius of the kernel (if applicable). sigma : tuple or float Std for the first gaussian (if applicable). sigma2 : tuple or float Std of a second gaussian (if present). Returns ------- kernel : array Filter kernel. """ ftype = ftype.lower(); o = se.structure_element_offsets(shape); mo = o.min(axis=1); shape = np.array(shape); if ftype == 'mean': # differnce of gaussians return np.ones(shape)/ shape.prod(); elif ftype == 'gaussian': if sigma == None: sigma = shape / 2. / math.sqrt(2 * math.log(2)); sigma = np.array(sigma); if len(sigma) < 3: sigma = np.array((sigma[0], sigma[0], sigma[0])); else: sigma = sigma[0:3]; g = np.mgrid[-o[0,0]:o[0,1], -o[1,0]:o[1,1], -o[2,0]:o[2,1]]; add = ((shape + 1) % 2) / 2.; x = g[0,:,:,:] + add[0]; y = g[1,:,:,:] + add[1]; z = g[2,:,:,:] + add[2]; ker = np.exp(-(x * x / 2. / (sigma[0] * sigma[0]) + y * y / 2. / (sigma[1] * sigma[1]) + z * z / 2. / (sigma[2] * sigma[2]))); return ker/ker.sum(); elif ftype == 'sphere': if radius == None: radius = mo; radius = np.array(radius); if len(radius) < 3: radius = np.array((radius[0], radius[0], radius[0])); else: radius = radius[0:3]; g = np.mgrid[-o[0,0]:o[0,1], -o[1,0]:o[1,1], -o[2,0]:o[2,1]]; add = ((shape + 1) % 2) / 2.; x = g[0,:,:,:] + add[0]; y = g[1,:,:,:] + add[1]; z = g[2,:,:,:] + add[2]; ker = 1 - (x * x / 2. / (radius[0] * radius[0]) + y * y / 2. / (radius[1] * radius[1]) + z * z / 2. / (radius[2] * radius[2])); ker[ker < 0] = 0.; return ker / ker.sum(); elif ftype == 'disk': if radius == None: radius = mo; radius = np.array(radius); if len(radius) < 3: radius = np.array((radius[0], radius[0], radius[0])); else: radius = radius[0:3]; g = np.mgrid[-o[0,0]:o[0,1], -o[1,0]:o[1,1], -o[2,0]:o[2,1]]; add = ((shape + 1) % 2) / 2.; x = g[0,:,:,:] + add[0]; y = g[1,:,:,:] + add[1]; z = g[2,:,:,:] + add[2]; ker = 1 - (x * x / 2. / (radius[0] * radius[0]) + y * y / 2. / (radius[1] * radius[1]) + z * z / 2. / (radius[2] * radius[2])); ker[ker < 0] = 0.; ker[ker > 0] = 1.0; return ker / ker.sum(); elif ftype == 'log': # laplacian of gaussians if sigma == None: sigma = shape / 4. / math.sqrt(2 * math.log(2)); sigma = np.array(sigma); if len(sigma) < 3: sigma = np.array((sigma[0], sigma[0], sigma[0])); else: sigma = sigma[0:3]; g = np.mgrid[-o[0,0]:o[0,1], -o[1,0]:o[1,1], -o[2,0]:o[2,1]]; add = ((shape + 1) % 2) / 2.; x = g[0,:,:,:] + add[0]; y = g[1,:,:,:] + add[1]; z = g[2,:,:,:] + add[2]; ker = np.exp(-(x * x / 2. / (radius[0] * radius[0]) + y * y / 2. / (radius[1] * radius[1]) + z * z / 2. / (radius[2] * radius[2]))); ker /= ker.sum(); arg = x * x / math.pow(sigma[0], 4) + y * y/ math.pow(sigma[1],4) + z * z / math.pow(sigma[2],4) - (1/(sigma[0] * sigma[0]) + 1/(sigma[1] * sigma[1]) + 1 / (sigma[2] * sigma[2])); ker = ker * arg; return ker - ker.sum()/len(ker); elif ftype == 'dog': if sigma2 == None: sigma2 = shape / 2. / math.sqrt(2 * math.log(2)); sigma2 = np.array(sigma2); if len(sigma2) < 3: sigma2 = np.array((sigma2[0], sigma2[0], sigma2[0])); else: sigma2 = sigma2[0:3]; if sigma == None: sigma = sigma2 / 1.5; sigma = np.array(sigma); if len(sigma) < 3: sigma = np.array((sigma[0], sigma[0], sigma[0])); else: sigma = sigma[0:3]; g = np.mgrid[-o[0,0]:o[0,1], -o[1,0]:o[1,1], -o[2,0]:o[2,1]]; add = ((shape + 1) % 2) / 2.; x = g[0,:,:,:] + add[0]; y = g[1,:,:,:] + add[1]; z = g[2,:,:,:] + add[2]; ker = np.exp(-(x * x / 2. / (sigma[0] * sigma[0]) + y * y / 2. / (sigma[1] * sigma[1]) + z * z / 2. / (sigma[2] * sigma[2]))); ker /= ker.sum(); sub = np.exp(-(x * x / 2. / (sigma2[0] * sigma2[0]) + y * y / 2. / (sigma2[1] * sigma2[1]) + z * z / 2. / (sigma2[2] * sigma2[2]))); return ker - sub / sub.sum(); else: raise ValueError('Filter type %r not valid!' % ftype);
############################################################################### ### Tests ############################################################################### def _test(): """Tests""" import ClearMap.ImageProcessing.Filter.FilterKernel as fk import ClearMap.Visualization.Plot3d as p3d k = fk.filter_kernel(ftype='dog', shape=(15,15,15), sigma=None, radius=None, sigma2=None); p3d.plot(k)