Code Monkey home page Code Monkey logo

ndarray-ndimage's People

Contributors

matsjoyce-refeyn avatar nilgoyette avatar

Stargazers

avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar avatar

Watchers

avatar

Forkers

jcmullwh matsjoyce-refeyn

ndarray-ndimage's Issues

Morphology and Kernel

We're currently using the Kernel3d enum (Star and Full) as the only choices for kernels.

The advantage is that we're quite fast because we know exactly where the "weights" are. This makes it possible to hardcode the conditions. For example, we're around 5 times faster than SciPy in morphology operations.

While this is nice, the problem is that the user can't choose it's own special kernel. I have personally never seen a non-start or non-full kernel choice in the wild, but it could happen.

Ideas:

  • A possible idea is to add an enum choice like Other(arr). We would know when it's a special kernel and we would be able to take a special path. It creates no useless allocation. I'm not a fan of templated enum however.
  • Another idea is to always receive an array for the kernel. At the start of the function, we analyze it and determine if it's a Star or Full. If it is, we take the fast paths. It forces us to allocate but we could remove the enum.

In brief, I want to keep being faster than SciPy and add all kernel choices. I tested an algo that offers more choices.

// Indices for the Star kernel
let indices = vec![(0, 1, 1), (1, 0, 1), (1, 1, 0), (1, 1, 1), (1, 1, 2), (1, 2, 1), (2, 1, 1)];
// Erosion
Zip::from(mask.windows((3, 3, 3)))
    .map_assign_into(zone, |w| !indices.iter().any(|&idx| !w[idx]));
// Dilation
Zip::from(mask.windows((3, 3, 3)))
    .map_assign_into(zone, |w| indices.iter().any(|&idx| w[idx]));

It's not too bad. We're still a little faster than SciPy and we're accepting all kernels. I just wonder if there's a way to keep our advantages :)

Is there an equivalent method of `scipy.ndimage.zoom`?

I'm not sure if it's possible to achieve the same effect of resampling a 3D image?

e.g.

In [249]: data = np.array([[1.,2.,3.],[4.,5.,6.],[7.,8.,9.]])

In [250]: from scipy.ndimage import zoom

In [251]: zoom(data, 2.0, mode='nearest')
Out[251]: 
array([[1.0, 1.29710751, 1.74855376, 2.25144624, 2.70289249, 3.0],
       [1.89132253, 2.18843004, 2.63987629, 3.14276878, 3.59421502, 3.89132253],
       [3.24566127, 3.54276878, 3.99421502, 4.49710751, 4.94855376, 5.24566127],
       [4.75433873, 5.05144624, 5.50289249, 6.00578498, 6.45723122, 6.75433873],
       [6.10867747, 6.40578498, 6.85723122, 7.36012371, 7.81156996, 8.10867747],
       [7.0 , 7.29710751, 7.74855376, 8.25144624, 8.70289249, 9.0]])

`convolve` expects both data and weights to have the exact same types

I am trying to convolve over the 2D arrays in an image cube and noticed the following mismatched types error.

pub fn denoise(hsi: &mut Array3<f32>) {
    let corner_kernel = arr2(&[[0.25f32, 0., 0.25], [0., 0., 0.], [0.25, 0., 0.25]]);
    hsi.outer_iter_mut().into_par_iter().for_each(|bnd| {
        let corner_average =
            ndarray_ndimage::convolve(&bnd, &corner_kernel.view(), BorderMode::Mirror, ORIGIN);
    });
// --snip--
}

Error:

error[E0308]: mismatched types
   --> src/lib.rs:93:45
    |
93  |             ndarray_ndimage::convolve(&bnd, &corner_kernel.view(), BorderMode::Mirror, ORIGIN);
    |             -------------------------       ^^^^^^^^^^^^^^^^^^^^^ types differ in mutability
    |             |
    |             arguments to this function are incorrect
    |
    = note: expected reference `&ArrayBase<ViewRepr<&mut f32>, Dim<[usize; 2]>>`
               found reference `&ArrayBase<ViewRepr<&f32>, Dim<[usize; 2]>>`

I belive this is because of convolve's signature. Since both data and weights have the same generic parameters, the compiler expects them to both have the exact same types, which, I suspect, is probably not necessary?

pub fn convolve<S, A, D>(
    data: &ArrayBase<S, D>,
    weights: &ArrayBase<S, D>,
    mode: BorderMode<A>,
    origin: isize
) -> Array<A, D>

Functions are not generic

Most functions are using Array when they should use ArrayBase. Thus it's impossible to call them with non-owned arrays (views).

Multithreaded Performance on Convolutions

Looking at the code it looks like it would be pretty simple to implement the par_iter() methods from ndarray into the convolution iterators to allow for multithreading in various processing routines

Morphology

There's a problem with the border management in the morphological erosion. Currently we are doing:

  • Set all borders to 0, then work on all 3x3x3 windows. Thus avoiding to work with outside voxels.
  • This seems logical, but it always gives a black border, whatever the kernel used.

The consequence of our current method can be seen in the following image.
Wrong

  • In red+pink we have ANTs ImageMath result.
  • in pink, we have ndimage result. We can easily see that we "ate" the border.

I read more and I realized that it's not the right behavior: "Pixels beyond the image border are assigned the maximum value afforded by the data type." Thus all the padded voxels should be set to 1.

Our current method (working with windows to avoid padding) can't work. We can keep the window but we must pad with 1 if we want to avoid border effects.

uniform_filter performance

The person who contributed uniform_filter wrote this

It's not going to be as quite as fast as scipy, because scipy's implementation hardcodes the fact that the weights are constant, but I really don't want to duplicate the border modes implementation.

Based on the code (https://github.com/scipy/scipy/blob/main/scipy/ndimage/src/ni_filters.c#L359-L362), scipy is using a single accumulator as it moves along lines, and just adds the new value and subtracts the oldest value each time the filter moves to a new sample. Then it just divides by the filter size, which should be a lot more efficient than a correlation. Looking again at the inner_collerate1d code, it should be possible to do something similar here, but it probably makes sense to do that in a subsequent PR for simplicity. I might have a go at that, it's made me curious now 😁.

No benchmark have been done/provided, but we can guess that uniform_filter can be optimized.

Padding performance

Compared to NumPy, pad is slow. correlate1d is less slow but we still need to make it faster. Those functions are particularly important because they are the basis for several other functions:

  • correlate1d: convolve1d, gaussian_filter/1d, prewitt, sobel
  • pad: correlate1d, min/max filters

I'm not sure what NumPy does to be so fast. Of course, Cython is used, so it's compiled, optimized and there's no bound checks. but still, we're around 2-6x slower, depending in the mode. IMO, we should be able to be as fast, if not faster. It's a mystery to me how NumPy manages to be so fast (around 10x) on the reflect, symmetric and wrap modes.

To replicate the results, please use

use std::time::Instant;

use ndarray::{Array1, Array3, ShapeBuilder};

use ndarray_ndimage::{pad, PadMode};

fn main() {
    println!("C order");
    let data_c = (0..200 * 200 * 200)
        .map(|v| v as u32)
        .collect::<Array1<_>>()
        .into_shape((200, 200, 200))
        .unwrap();
    test(&data_c);

    println!("\nF order");
    let mut data_f = Array3::zeros((200, 200, 200).f());
    data_f.assign(&data_c);
    test(&data_f);
}

fn test(data: &Array3<u32>) {
    let now = Instant::now();
    let padded = pad(&data, &[[1, 1]], PadMode::Constant(0));
    let elapsed = now.elapsed();
    println!(
        "Constant 0 {}s {}ms  {}",
        elapsed.as_secs(),
        elapsed.subsec_millis(),
        padded[(0, 0, 0)]
    );

    let now = Instant::now();
    let padded = pad(&data, &[[1, 1]], PadMode::Constant(42));
    let elapsed = now.elapsed();
    println!(
        "Constant 42 {}s {}ms  {}",
        elapsed.as_secs(),
        elapsed.subsec_millis(),
        padded[(0, 0, 0)]
    );

    let now = Instant::now();
    let padded = pad(&data, &[[1, 1]], PadMode::Edge);
    let elapsed = now.elapsed();
    println!("Edge      {}s {}ms  {:?}", elapsed.as_secs(), elapsed.subsec_millis(), padded.dim());

    let now = Instant::now();
    let padded = pad(&data, &[[1, 1]], PadMode::Maximum);
    let elapsed = now.elapsed();
    println!("Maximum   {}s {}ms  {:?}", elapsed.as_secs(), elapsed.subsec_millis(), padded.dim());

    let now = Instant::now();
    let padded = pad(&data, &[[1, 1]], PadMode::Mean);
    let elapsed = now.elapsed();
    println!("Mean      {}s {}ms  {:?}", elapsed.as_secs(), elapsed.subsec_millis(), padded.dim());

    let now = Instant::now();
    let padded = pad(&data, &[[1, 1]], PadMode::Median);
    let elapsed = now.elapsed();
    println!("Median    {}s {}ms  {:?}", elapsed.as_secs(), elapsed.subsec_millis(), padded.dim());

    let now = Instant::now();
    let padded = pad(&data, &[[1, 1]], PadMode::Minimum);
    let elapsed = now.elapsed();
    println!("Minimum   {}s {}ms  {:?}", elapsed.as_secs(), elapsed.subsec_millis(), padded.dim());

    let now = Instant::now();
    let padded = pad(&data, &[[1, 1]], PadMode::Reflect);
    let elapsed = now.elapsed();
    println!("Reflect   {}s {}ms  {:?}", elapsed.as_secs(), elapsed.subsec_millis(), padded.dim());

    let now = Instant::now();
    let padded = pad(&data, &[[1, 1]], PadMode::Symmetric);
    let elapsed = now.elapsed();
    println!("Symmetric {}s {}ms  {:?}", elapsed.as_secs(), elapsed.subsec_millis(), padded.dim());

    let now = Instant::now();
    let padded = pad(&data, &[[1, 1]], PadMode::Wrap);
    let elapsed = now.elapsed();
    println!("Wrap      {}s {}ms  {:?}", elapsed.as_secs(), elapsed.subsec_millis(), padded.dim());
}

And for Python

from datetime import datetime
import time

import numpy as np


def test(a):
    start = datetime.now()
    np.pad(a, 1, mode='constant', constant_values=42)
    print(' constant', datetime.now() - start)

    start = datetime.now()
    np.pad(a, 1, mode='edge')
    print('     edge', datetime.now() - start)

    start = datetime.now()
    np.pad(a, 1, mode='maximum')
    print('  maximum', datetime.now() - start)

    start = datetime.now()
    np.pad(a, 1, mode='mean')
    print('     mean', datetime.now() - start)

    start = datetime.now()
    np.pad(a, 1, mode='median')
    print('   median', datetime.now() - start)

    start = datetime.now()
    np.pad(a, 1, mode='minimum')
    print('  minimum', datetime.now() - start)

    start = datetime.now()
    np.pad(a, 1, mode='reflect')
    print('  reflect', datetime.now() - start)

    start = datetime.now()
    np.pad(a, 1, mode='symmetric')
    print('symmetric', datetime.now() - start)

    start = datetime.now()
    np.pad(a, 1, mode='wrap')
    print('     wrap', datetime.now() - start)
    

print('C order')
a = np.arange(200*200*200).reshape((200, 200, 200)).astype(np.uint32)
test(a)

print('\nF order')
a = np.asfortranarray(a)
test(a)

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.