0910---scipy.signal.argrelextrema 计算多维数组局部极值问题

scipy signal argrelextrema 计算多维数组局部极值问题

在科学计算中，求一组数据的极值非常普遍，利用python scipy 工具包可以非常高效的计算局部极值，但在处理多维数组数据时，因选取的axis轴向不同，结果具有不同的意义，在理解上也比较困难。本文结合实例，进行透彻分析。

运行环境： Anaconda3

一 1维数组

1 引入工具包

import numpy as np
from scipy.signal import argrelextrema

2 案例

# 例子1
x = np.array([2, 1, 2, 3, 2, 0, 1, 0])
trema = argrelextrema(x, np.greater)
trema

# 输出结果
out (array([3, 6], dtype=int64),)

print('return type:',type(e_1D))
print('tuple lenth:',len(e_1D))
print('tuple element type:',type(e_1D[0]))

return type: <class 'tuple'>
tuple lenth: 1
tuple element type: <class 'numpy.ndarray'>
first data of first tuple element: 3

解释：由于上述例子是求一维数组的局部极大值，结果是返回一个元组，返回元组的元素类型是numpy.ndarray，里面存储局部极值的位置标记，即索引编号为3和6的位置标记信息。因一维数组只有一个X轴(0轴，如果给定 axis=1，会产生错误。

二 2维数组

对于2维数组求极值问题，官方文档给出了一个非常容易引起歧义的例子，下面就该例子进行分析。

# 例子2
y = np.array([ [1, 2, 1, 2],
               [2, 2, 0, 0],
               [5, 3, 4, 4]])
argrelextrema(y, np.less, axis=1)

# out 输出
(array([0, 2], dtype=int64), array([2, 1], dtype=int64))

说明：当 axis=1时，按y轴的方向进行局部极值判断，即以行为单元进行数据比较。

由此可见，二维数组的局部极小值有2个，一个是位于第0行第2列的1，另一个是第2行第一列的3，以蓝色圆圈标出，与输出结果 (array([0, 2], dtype=int64), array([2, 1], dtype=int64)) 相符。

当 axis =0 时，按x轴的方向进行局部极小值判断，即以列为单元进行数值比较。

# 例子3
argrelextrema(y, np.less, axis=0)
# out 输出 
(array([1, 1], dtype=int64), array([2, 3], dtype=int64))

结果应该是元素(1，2)和元素(1,3), 但输出的结果是 (array([1, 1], dtype=int64), array([2, 3], dtype=int64))，怎么都对应不上。

查找官方文档，signal.argrelextrema 的返回值有如下解释

返回值extrema : tuple of ndarrays
Indices of the maxima in arrays of integers. extrema[k] is
the array of indices of axis k of data. Note that the
return value is a tuple even when data is 1-D.

返回值extrema：由 ndarrays数组作为元素组成的元组，极值的位置由ndarrays数组中的整数值给出。extrema[k] 代表数据在k轴上的位置标记数组。注意，在数据是1维数组的情况下，返回结果仍然是元组。

对于例子3，返回结果 (array([1, 1], dtype=int64), array([2, 3], dtype=int64))是含有2个元素的元组，

extrema[0]代表的是 x 轴极值位置标记(行标记)，extrema[1]代表的是y轴极值位置标记(列标记)。转换维人们习惯的行列坐标形式应该是：(extrema[0] [0], extrema[1] [0]) , (extrema[0] [1], extrema[1] [1]) , 即 (1,2) , (1,3).

回过头来解释例子2，返回结果是 (array([0, 2], dtype=int64), array([2, 1], dtype=int64))，对应的行列坐标是(0,2),(2,1),正好和返回的元组元素一样，这就是官方例子的不恰当之处。

三 2维数组的扩展案例

1 数据及运行结果(axis=1)

y1 = np.array([
    [ 1 , 2 , 1 , 2 , 1 ],
    [ 2 , 2 , 0 , 5 , 0 ],
    [ 5 , 3 , 4 , 4 , 2 ],
    [ 1 , 2 , 1 , 2 , 3 ],
    [ 4 , 3 , 0 , 3 , 2 ],
    [ 6 , 2 , 7 , 1 , 0 ]
    ])
    
extrema =  argrelextrema(y1, np.less, axis=1)

(array([0, 1, 2, 3, 4, 5], dtype=int64),
 array([2, 2, 1, 2, 2, 1], dtype=int64))
 对应的行列坐标是：(0,2),(1,2),(2,1),(3,2),(4,2),(5,1)

2 图示(axis=1)

3 数据及运行结果(axis=0)

y1 = np.array([
    [ 1 , 2 , 1 , 2 , 1 ],
    [ 2 , 2 , 0 , 5 , 0 ],
    [ 5 , 3 , 4 , 4 , 2 ],
    [ 1 , 2 , 1 , 2 , 3 ],
    [ 4 , 3 , 0 , 3 , 2 ],
    [ 6 , 2 , 7 , 1 , 0 ]
    ])
    
extrema = argrelextrema(y1, np.less, axis=0)

(array([1, 1, 3, 3, 3, 4], dtype=int64),
 array([2, 4, 0, 1, 3, 2], dtype=int64))
 
 对应的行列坐标是：(1,2),(1,4),(3,0),(3,1),(3,3),(4,2)

4 图示(axis=0)

                              2022-09-12