python - scipy.ndimage.laplace 中的边界模式

Question

函数scipy.ndimage.laplace可用于计算应用于 N 维数组的拉普拉斯算子。如果要使用此功能，例如，用于物理应用，关键部分通常是如何处理边界。

查看函数文档，支持的选项有很多:

• “反射(reflect)”(d c b a | a b c d | d c b a)
• ‘常量’ (k k k k | a b c d | k k k k)
• “最近”(a a a a | a b c d | d d d d)
• ‘镜像’ (d c b | a b c d | c b a)
• ‘包裹’(a b c d | a b c d | a b c d)

我个人缺少的一件事是“平稳延续” 选项。在上面的表示法中，它将对应于镜像选项的一个稍微修改的版本:

• '平稳连续'((a-d)(a-c)(a-b)| a b c d |(d-c)(d-b)(d-a))

这样做的动机是不在边界处引入扭结，所有支持的选项似乎都是针对任意输入做的。

问题:在scipy/numpy中有什么办法可以做这个拉普拉斯运算吗？

Answer 1

我深入研究了 scipy 源代码。我的发现让我认为仅使用 python scypi/numpy 添加您的自定义模式并不容易:您最终将在 python 中重新编码 scipy 的大部分主干代码。

您所指的模式深入到底层 scipy 的 C 代码中进行管理。您将需要更改 C 枚举值，根据您的自定义模式逻辑更新源文件，并生成一个新包。

这可能是要在 github 上打开的问题/演变:要么允许自定义构建数组方法，要么实现新模式。或者每人一张票。

参见 here .引用:

/* Extend a line in memory to implement boundary conditions: */
int NI_ExtendLine(double *buffer, npy_intp line_length,
                  npy_intp size_before, npy_intp size_after,
                  NI_ExtendMode extend_mode, double extend_value)
{
    double *first = buffer + size_before;
    double *last = first + line_length;
    double *src, *dst, val;

    switch (extend_mode) {
        /* aaaaaaaa|abcd|dddddddd */
        case NI_EXTEND_NEAREST:
            src = first;
            dst = buffer;
            val = *src;
            while (size_before--) {
                *dst++ = val;
            }
            src = last - 1;
            dst = last;
            val = *src;
            while (size_after--) {
                *dst++ = val;
            }
            break;
        /* abcdabcd|abcd|abcdabcd */
        case NI_EXTEND_WRAP:
            src = last - 1;
            dst = first - 1;
            while (size_before--) {
                *dst-- = *src--;
            }
            src = first;
            dst = last;
            while (size_after--) {
                *dst++ = *src++;
            }
            break;
        /* abcddcba|abcd|dcbaabcd */
        case NI_EXTEND_REFLECT:
            src = first;
            dst = first - 1;
            while (size_before && src < last) {
                *dst-- = *src++;
                --size_before;
            }
            src = last - 1;
            while (size_before--) {
                *dst-- = *src--;
            }
            src = last - 1;
            dst = last;
            while (size_after && src >= first) {
                *dst++ = *src--;
                --size_after;
            }
            src = first;
            while (size_after--) {
                *dst++ = *src++;
            }
            break;
        /* cbabcdcb|abcd|cbabcdcb */
        case NI_EXTEND_MIRROR:
            src = first + 1;
            dst = first - 1;
            while (size_before && src < last) {
                *dst-- = *src++;
                --size_before;
            }
            src = last - 2;
            while (size_before--) {
                *dst-- = *src--;
            }
            src = last - 2;
            dst = last;
            while (size_after && src >= first) {
                *dst++ = *src--;
                --size_after;
            }
            src = first + 1;
            while (size_after--) {
                *dst++ = *src++;
            }
            break;
        /* kkkkkkkk|abcd]kkkkkkkk */
        case NI_EXTEND_CONSTANT:
            val = extend_value;
            dst = buffer;
            while (size_before--) {
                *dst++ = val;
            }
            dst = last;
            while (size_after--) {
                *dst++ = val;
            }
            break;
        default:
            PyErr_Format(PyExc_RuntimeError,
                         "mode %d not supported", extend_mode);
            return 0;
    }
    return 1;
}

请注意，我绝对不是 scipy 方面的专家，所以也许我遗漏了一些重要的东西，或者有非常有据可查的原因来解释建议的模式。

---

编辑 2019-Ap-24 一些 documentation here我没有得到所有的细节，但据我了解:

• 它指的是 Philippe Thevenez 博士的学术著作，作为算法某些部分的潜在来源。
• 这张ticket似乎表明 future 可能会开发一些其他的模式。

This ticket还处理新模式。

鉴于门票的开放日期，它可能需要一些时间才能到来。