Skip to content

add horizontal padding for sinogram Fourier-Wavelets deringer

Pierre Paleo requested to merge munch_padding into master

About

The Fourier-Wavelets deringer ("Munch") often creates "ripple" artefacts.

These artefacts are, to a certain extent, due to the wavelets decomposition with periodic boundary mode. A "aliasing" occurs when reconstructing after filtering the coefficients.

Note that the Fourier transform of the coefficients does not seem to create the artefacts (see figures below), which makes sense as it's performed along the angles (not the horizontal direction).

To do

  • Add padding to MunchDeringer
  • Add padding to CudaMunchDeringer
  • Unit tests
  • Integrate in pipeline
  • End-to-end reconstruction test

Notes

Here are some padding attempts. In the following figures, padding ((v1, v2), (h1, h1)) means that the image is extended by v1 + v2 in the vertical dimension (axis 0), and h1 + h2 in the vertical direction.

a1

a2

Line profiles in the reconstruction indicate that artefacts in the edges occur when using no padding and vertical padding, but disappear when using horizontal padding.

a3

a4

Code:

import numpy as np
from tomoscan.io import HDF5File
from spire.utils import ims
from nabu.reconstruction.fbp import Backprojector
from nabu.thirdparty.pore3d_deringer_munch import munchetal_filter
from silx.image.tomography import get_next_power
from algotom.prep.removal import remove_all_stripe, remove_large_stripe, remove_stripe_based_2d_filtering_sorting, remove_stripe_based_fft, remove_stripe_based_filtering, remove_stripe_based_fitting



def read_slice(fname, h5path, i=0):
    with HDF5File(fname, "r") as f:
        data = f[h5path][i]
    return data

s = read_slice("/data/scisofttmp/paleo/debug/AP/vpan_s08_11_100nm_/rec/sinogram_no_deringer_vpan_s08_11_100nm_rec_.hdf5", "entry/sinogram/results/data", (slice(None, None), slice(None, None), (slice(None, None))))[:, 0, :]


B = Backprojector(s.shape, angles=-np.linspace(0, -np.pi, s.shape[0], False), rot_center=1604, padding_mode="edge", extra_options={"centered_axis": True, "clip_outer_circle": True})
def clip_std(arr, n=3):
    m = np.mean(arr)
    s = np.std(arr)
    return np.clip(arr, m - n*s, m + n*s)


def munch_remove_rings_with_padding(sino, sigma, levels, pad_width):
    n_a, n_x = sino.shape
    n_a2 = get_next_power(n_a * 2)

    sino_p = np.pad(sino, pad_width, mode="edge")
    sino_corr = munchetal_filter(sino_p, levels, sigma)

    ((v1, v2), (h1, h2)) = pad_width
    v1 = v1 or None
    v2 = -v2 if v2 > 0 else None
    h1 = h1 or None
    h2 = -h2 if h2 > 0 else None
    return np.ascontiguousarray(sino_corr[v1:v2, h1:h2], dtype="f")

sinos_deringed = {}
for pw in [
    ((0, 0), (0, 0)), # no padding (was default)
    ((1000, 1000), (0, 0)), # pad only angles
    ((0, 0), (1608, 1608)), # pad only horizontally
    ((1000, 1000), (1608, 1608)), # pad in both directions
]:
    sinos_deringed[pw] = np.ascontiguousarray(munch_remove_rings_with_padding(s, 0.8, 4, pw), dtype="f")
Edited by Pierre Paleo

Merge request reports