Skip to content

Commit

Permalink
Added the gaussian_filter source code for inspiration
Browse files Browse the repository at this point in the history
  • Loading branch information
@aaschwanden
aaschwanden committed Jun 5, 2024
1 parent a55f981 commit 6bc28af
Showing 1 changed file with 201 additions and 0 deletions.
201 changes: 201 additions & 0 deletions notebooks/dem_smoothing.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -311,6 +311,183 @@
"Maybe we can use https://github.com/scipy/scipy/blob/v1.13.1/scipy/ndimage/_filters.py#L186 for inspiration. We need to compute $\\sigma$ scaled by ice thickness in a sensible way. "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "01e19cb4-64bd-4024-a05a-45439ada7c6f",
"metadata": {},
"outputs": [],
"source": [
"def _gaussian_weighted_kernel1d(l, H, x):\n",
" \"\"\"\n",
" Computes a 1-D ice-thickness weighted Gaussian convolution kernel.\n",
" \"\"\"\n",
"\n",
" L = l * (H @ H.T)**(1/2) # I am not sure this is correct\n",
" sigma = (1 + D**2 / (2 * L**2))**(-1)\n",
" sigma2 = sigma * sigma\n",
" phi_x = numpy.exp(-0.5 / sigma2 * x ** 2)\n",
" phi_x = phi_x / phi_x.sum()\n",
"\n",
" return phi_x\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "707bfde1-f8b0-497e-b984-9f192d3b3b6c",
"metadata": {},
"outputs": [],
"source": [
"def gaussian_weighted_filter1d(input, l, H, x, axis=-1, order=0, output=None,\n",
" mode=\"reflect\", cval=0.0):\n",
" \"\"\"1-D Gaussian filter.\n",
"\n",
" Parameters\n",
" ----------\n",
" %(input)s\n",
" l : number of ice thicknesses\n",
" ice thickness multiplier\n",
" %(axis)s\n",
" order : int, optional\n",
" An order of 0 corresponds to convolution with a Gaussian\n",
" kernel. A positive order corresponds to convolution with\n",
" that derivative of a Gaussian.\n",
" %(output)s\n",
" %(mode_reflect)s\n",
" %(cval)s\n",
"\n",
" Returns\n",
" -------\n",
" gaussian_weighted_filter1d : ndarray\n",
"\n",
"\n",
" \"\"\"\n",
" sd = float(sigma)\n",
" # make the radius of the filter equal to truncate standard deviations\n",
" lw = int(truncate * sd + 0.5)\n",
" if radius is not None:\n",
" lw = radius\n",
" if not isinstance(lw, numbers.Integral) or lw < 0:\n",
" raise ValueError('Radius must be a nonnegative integer.')\n",
" # Since we are calling correlate, not convolve, revert the kernel\n",
" weights = _gaussian_weighted_kernel1d(l, H, x)[::-1]\n",
" return correlate1d(input, weights, axis, output, mode, cval, 0)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1ee43d50-da19-41aa-a2b6-cdc447223e84",
"metadata": {},
"outputs": [],
"source": [
"def gaussian_weighted_filter(input, l, H, x, output=None,\n",
" mode=\"reflect\", cval=0.0,\n",
" axes=None):\n",
" \"\"\"Multidimensional Gaussian filter.\n",
"\n",
" Parameters\n",
" ----------\n",
" %(input)s\n",
" sigma : scalar or sequence of scalars\n",
" Standard deviation for Gaussian kernel. The standard\n",
" deviations of the Gaussian filter are given for each axis as a\n",
" sequence, or as a single number, in which case it is equal for\n",
" all axes.\n",
" order : int or sequence of ints, optional\n",
" The order of the filter along each axis is given as a sequence\n",
" of integers, or as a single number. An order of 0 corresponds\n",
" to convolution with a Gaussian kernel. A positive order\n",
" corresponds to convolution with that derivative of a Gaussian.\n",
" %(output)s\n",
" %(mode_multiple)s\n",
" %(cval)s\n",
" truncate : float, optional\n",
" Truncate the filter at this many standard deviations.\n",
" Default is 4.0.\n",
" radius : None or int or sequence of ints, optional\n",
" Radius of the Gaussian kernel. The radius are given for each axis\n",
" as a sequence, or as a single number, in which case it is equal\n",
" for all axes. If specified, the size of the kernel along each axis\n",
" will be ``2*radius + 1``, and `truncate` is ignored.\n",
" Default is None.\n",
" axes : tuple of int or None, optional\n",
" If None, `input` is filtered along all axes. Otherwise,\n",
" `input` is filtered along the specified axes. When `axes` is\n",
" specified, any tuples used for `sigma`, `order`, `mode` and/or `radius`\n",
" must match the length of `axes`. The ith entry in any of these tuples\n",
" corresponds to the ith entry in `axes`.\n",
"\n",
" Returns\n",
" -------\n",
" gaussian_filter : ndarray\n",
" Returned array of same shape as `input`.\n",
"\n",
" Notes\n",
" -----\n",
" The multidimensional filter is implemented as a sequence of\n",
" 1-D convolution filters. The intermediate arrays are\n",
" stored in the same data type as the output. Therefore, for output\n",
" types with a limited precision, the results may be imprecise\n",
" because intermediate results may be stored with insufficient\n",
" precision.\n",
"\n",
" The Gaussian kernel will have size ``2*radius + 1`` along each axis. If\n",
" `radius` is None, the default ``radius = round(truncate * sigma)`` will be\n",
" used.\n",
"\n",
" Examples\n",
" --------\n",
" >>> from scipy.ndimage import gaussian_filter\n",
" >>> import numpy as np\n",
" >>> a = np.arange(50, step=2).reshape((5,5))\n",
" >>> a\n",
" array([[ 0, 2, 4, 6, 8],\n",
" [10, 12, 14, 16, 18],\n",
" [20, 22, 24, 26, 28],\n",
" [30, 32, 34, 36, 38],\n",
" [40, 42, 44, 46, 48]])\n",
" >>> gaussian_filter(a, sigma=1)\n",
" array([[ 4, 6, 8, 9, 11],\n",
" [10, 12, 14, 15, 17],\n",
" [20, 22, 24, 25, 27],\n",
" [29, 31, 33, 34, 36],\n",
" [35, 37, 39, 40, 42]])\n",
"\n",
" >>> from scipy import datasets\n",
" >>> import matplotlib.pyplot as plt\n",
" >>> fig = plt.figure()\n",
" >>> plt.gray() # show the filtered result in grayscale\n",
" >>> ax1 = fig.add_subplot(121) # left side\n",
" >>> ax2 = fig.add_subplot(122) # right side\n",
" >>> ascent = datasets.ascent()\n",
" >>> result = gaussian_filter(ascent, sigma=5)\n",
" >>> ax1.imshow(ascent)\n",
" >>> ax2.imshow(result)\n",
" >>> plt.show()\n",
" \"\"\"\n",
" input = numpy.asarray(input)\n",
" output = _ni_support._get_output(output, input)\n",
"\n",
" axes = _ni_support._check_axes(axes, input.ndim)\n",
" num_axes = len(axes)\n",
" orders = _ni_support._normalize_sequence(order, num_axes)\n",
" sigmas = _ni_support._normalize_sequence(sigma, num_axes)\n",
" modes = _ni_support._normalize_sequence(mode, num_axes)\n",
" radiuses = _ni_support._normalize_sequence(radius, num_axes)\n",
" axes = [(axes[ii], sigmas[ii], orders[ii], modes[ii], radiuses[ii])\n",
" for ii in range(num_axes) if sigmas[ii] > 1e-15]\n",
" if len(axes) > 0:\n",
" for axis, sigma, order, mode, radius in axes:\n",
" gaussian_weighted_filter1d(input, l, H, x, axis, output,\n",
" mode, cval)\n",
" input = output\n",
" else:\n",
" output[...] = input[...]\n",
" return output"
]
},
{
"cell_type": "markdown",
"id": "2dca7e41-2e62-4576-8996-13b25bafa88d",
Expand All @@ -328,6 +505,30 @@
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "624e6cb3-5752-4513-92c4-649cb195f4c2",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "a0d0802b-22f4-40db-bb09-f5d443ccd645",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "dc007bbe-44d6-4d4d-aef7-3cacf88405ba",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand Down

0 comments on commit 6bc28af

Please sign in to comment.