如何提高 Python 数组操作性能.

最近正在写一个小系统. 涉及到把 CT 影像转换成 PNG.
CT 影像不是不是标准通用格式, 所以需要从读取 bytes 再转成数组..
这个过程中需要两次处理数组数据, 但是性能很慢..有相关经验的朋友可以给点建议嘛?
以下是详述:

我读出来的 bytes 数据处理成 uint8 后. 是这样的形式:
[208, 4, 208, 4, 208, 4...196, 8]
不懂图像处理的知识, 我的理解就是, 一个灰度, 一个 Alpha 通道值(透明值?).
第一次处理数组是要把上面的数组, 改成如下形式: <- 暂且叫生成数组阶段.
[208, 208, 208, 4, 208, 208, 208, 4.......]
就是把 208 这个灰度值变成 rgb 的形式.
然后再通过一个循环变成 pillow 支持的格式, 如下: <- 暂且叫数组转换阶段吧.
[
[ [r, g, b, a], [r, g, b, a], [r, g, b, a], ],
[ [r, g, b, a], [r, g, b, a], [r, g, b, a], ],
[ [r, g, b, a], [r, g, b, a], [r, g, b, a], ],
]

然后这个数组要通过 numpy 转换成 uint8 类型才能提交给 pillow 或者类似图片处理的库..

所以大致耗时分为三个阶段:
1. 生成数组
2. 转换数组
3. 通过 numpy 转换数组类型为 uint8.

数组大小为: 5022 * 4200 * 4 = 84,369,600

环境 A 介绍:
Win10, 10900K, 64GB 3600, NVME 2TB SSD
Python 3.10.4, pillow 最新版
阶段一耗时: 9 秒左右
阶段二耗时: 45 秒左右
阶段三耗时: 5 秒左右

觉得太慢, 尝试全部改成 numpy 操作.
但是阶段二慢的出奇, 后来查资料发现, numpy 在大量数组下标赋值操作的性能还不如 Python 原生 list.
阶段二, 上厕所, 洗把脸回来还没跑完, 就直接放弃了. 也尝试过 np.insert 操作, 一样慢.

后来尝试环境 B, 基于 pypy 的.
pypy3.8, numpy, pillow, 同样的硬件.
第一阶段耗时: 0.2 秒
第二阶段耗时: 9.7 秒
第三阶段耗时: 77 秒
查阅资料得知, pypy 的 C 扩展接口性能很差, 不如原生 python.
而 Numpy 就是 C 扩展的库..所以导致 Numpy 性能急剧下降.

好了..剩下的办法超出我的知识范畴了..
有朋友能分享一下经验或者给个可能的方向我去看也行.

需要转换的数据量大约有 23 万个 CT 文件. 按照一个文件 1 分钟, 我即便 4 个进程跑(文件转换的服务器 4 核 8GB 内存, 还要跑一些其他服务.), 5 万分钟 / 60 分钟 / 24 小时一点意外没有的情况要跑 42 天...

顺便吐槽以下, 朋友给写了个 nodejs 的版本, 转换一个图 2.3 秒, WTF!!

NumPy

耗时

Python

pillow

70 条回复 • 2022-09-25 15:15:38 +08:00

dlsflh

2022-09-23 16:09:32 +08:00

我不太懂技术，不过那么多深度学习医学影像相关的项目是怎么做的呢？是不是可以参考他们的方法？

hsfzxjy

2022-09-23 16:10:09 +08:00 via Android

建议直接从 bytes 读成一个 numpy array ，这段用 cython 来写，会快得很多，而且不是很麻烦

hsfzxjy

2022-09-23 16:15:42 +08:00 via Android

用 cython 的话，你相当于可以直接用类 python 的语法写个多重循环，把读出来的 byte 放到 array 的对应位置。不用去想什么向量操作，但编译出来会接近 C 的速度。一些越界检查什么的需要关掉以获得最佳性能

MoYi123

2022-09-23 16:17:28 +08:00

你的问题我估计肯定是能直接调 numpy 的 api 或者其他的一些图像处理库解决的, 但是大家都不一定了解你的 CT 影像具体的细节, 所以你要先说明一下慢的那一步你具体是怎么做的. 而不是一个简单的“数组转换阶段”就概括了.

jeeyong

2022-09-23 16:20:22 +08:00

@dlsflh 感谢回复..
我猜..他们只是读...并不涉及到数组大范围的修改, 或者大量的下标赋值操作吧?所以并没有那么慢..
哦对, 另外可能还直接调用了 GPU 做计算...而我目前不能这么做..服务器的显卡很垃圾.

@hsfzxjy cython....是不是不如直接 nodejs 搞一个了...研究以下打包发布, python 直接拉起一个 node 的进程后台转换..这个东西有开发周期的约束, 快到了..

lizytalk

2022-09-23 16:22:24 +08:00

可以试一下 numba?

jeeyong

2022-09-23 16:22:50 +08:00

@MoYi123
我读出来的 bytes 数据处理成 uint8 后. 是这样的形式:
[208, 4, 208, 4, 208, 4...196, 8]
不懂图像处理的知识, 我的理解就是, 一个灰度, 一个 Alpha 通道值(透明值?).
第一次处理数组是要把上面的数组, 改成如下形式: <- 暂且叫生成数组阶段.
就是再赋值两次灰度..写入数组.构建成如下形式:

[208, 208, 208, 4, 208, 208, 208, 4.......]

再下一步就是把 208 这个灰度值变成 rgb 的形式.
然后再通过一个循环变成 pillow 支持的格式, 如下: <- 暂且叫数组转换阶段吧.
[
[ [r, g, b, a], [r, g, b, a], [r, g, b, a], ],
[ [r, g, b, a], [r, g, b, a], [r, g, b, a], ],
[ [r, g, b, a], [r, g, b, a], [r, g, b, a], ],
]

我补充了点文字..你看一下能理解吗?

yoohwzy

2022-09-23 16:23:09 +08:00 via iPhone

不要用 for 循环来处理

```Python
w = 5022
h = 4200

img = np.asarray(img).reshape((w * h, -1))

rgba_img = np.stack((img[:, 0],

lookStupiToForce

2022-09-23 16:23:19 +08:00

看你的意思，是转换 5022*4200 像素的灰度图片，到 rgb 色彩空间？

往这方面搜搜，看有没有结果
www.google.com/search?q=python+图像+色彩空间+转换

上面第一条搜索结果 https://blog.csdn.net/hallobike/article/details/120385129 里，有“读入灰度图片”

中文搜不到也可以去搜英文
www.google.com/search?q=python+color+space+convert

jeeyong

2022-09-23 16:24:42 +08:00

@lizytalk 试过.. 仅限于 @jit 装饰. 一堆警告和报错.
主要是不支持一些方法. 包括不限于 numpy 的...
网上翻的教程, 所有的数组转 uint8 类型, 都是 numpy 直接干的...
其他我就不会了...有其他方法我也想试试.... 不用 numpy 转换, 或许可以直接上 pypy

hsfzxjy

2022-09-23 16:24:58 +08:00

@jeeyong Cython 不难的，花一个小时看看文档就能上手，写出来肯定比 nodejs 快，因为只涉及简单的内存存储。你用 nodejs 写还要考虑进程间通信，以及两个语言中格式的切换，这部分开销也不小。

threebr

2022-09-23 16:25:41 +08:00

不懂 numpy 为什么会慢。用 numpy 的话要用向量操作代替循环，也就是说你处理一张图片的时候（你的三个阶段）是不写 for 循环的，然后底层 numpy 可能是用 c 编译的库，可能直接用 Intel 或者 AMD 专门处理向量的二进制库，不应该比你自己拿 c 写更慢

jeeyong

2022-09-23 16:26:53 +08:00

@yoohwzy
@lookStupiToForce 感谢回复..先尝试 np.stack 的用法...
再看 google 结果...
对..如果不用转换, 直接读取生成, 是正确的思路..只是当时没找到方法, 就先硬上了..
先解决,再优化~

jeeyong

2022-09-23 16:28:58 +08:00

@hsfzxjy 是吗? 一直很怕上 C 相关的东西...哈哈
nodejs 可以直接读数据库获取路径, 去转换完了写回对应字段...
解决问题来说, 应该这个方法是最快的..
但是不服气啊...60 秒和 2.3 秒的差距...

stein42

2022-09-23 16:29:49 +08:00

直接用 numpy 就可以了，python 的循环太慢。
```
import numpy as np

width = 5022
height = 4200

# buffer 直接从文件读取
buffer = bytes([208, 4]) * (width * height)
assert len(buffer) == width * height * 2

b = np.frombuffer(buffer, dtype=np.uint8)
grey = b[::2][..., None]
alpha = b[1::2][..., None]
image = np.hstack((grey, grey, grey, alpha)).reshape((width, height, 4))
```

jeeyong

2022-09-23 16:30:08 +08:00

@threebr 我的理解是, python 要通过对应的接口和 numpy 通信, 而 numpy 本身对于下标赋值这种操作的过程是很复杂的.
这基本就是昨天看到的一篇文章的原话. 作者还列出了对应的处理过程, 因为是 C, 直接放弃阅读了..

yoohwzy

2022-09-23 16:30:32 +08:00 via iPhone

不要用 for 循环来处理

```Python
w = 5022
h = 4200
c = 4

img = np.asarray(img).reshape((w * h, -1))

# 转换 GA 图片到 RGBA
rgba_img = np.stack((img[:, 0], img[:, 0], img[:, 0], img[:, 1]), axis=1)
# 转换数据格式为 np.uint8 并生成 HWC 排列的图片
rgba_img = rgba_img.astype(np.uint8).reshape((h, w, c))

# 如果需要 CHW 排列的
rgba_img = rgba_img.transpose((2, 0, 1))
```

jeeyong

2022-09-23 16:32:30 +08:00

@stein42 卧槽卧槽, 果然大神多....
我需要消化一下代码...
这之前都没用过 numpy...感谢感谢..

另外请教一下,
buffer = bytes([208, 4]) * (width * height)
文件内容不只是 [208, 4] 还有其他的灰度值和通道值, 这一步不理解..

hsfzxjy

2022-09-23 16:32:35 +08:00

看你新的描述应该用不到 cython ，我盲写了一段，你可以根据文件格式调一下

```python
import numpy as np

with open("file", "rb") as f:
<TAB>img = f.read()

img = np.frombuffer(img, dtype=np.uint8) # 2HW
img = img.reshape(H, W, 2) # H x W x 2
img = np.stack([img[..., 0], img[..., 0], img[..., 0], img[..., 1]]) # H x W x 4
```

lmshl

2022-09-23 16:32:58 +08:00

我是直接写个 tensor 扔给 torch 跑，我知道它会自动用 SIMD 或者 CUDA (如果启动 gpu 的话)

hsfzxjy

2022-09-23 16:34:02 +08:00

@hsfzxjy #19 修正 np.stack([...], axis=-1)

jeeyong

2022-09-23 16:34:15 +08:00

@dlsflh
我是直接写个 tensor 扔给 torch 跑，我知道它会自动用 SIMD 或者 CUDA (如果启动 gpu 的话)

@lmshl 这大哥的做法解释了你的疑问..

jeeyong

2022-09-23 16:35:27 +08:00

各位的代码, 够我消化一夜了....
太感谢了...太感谢了!
@hsfzxjy
@yoohwzy
@stein42

wcsjtu

2022-09-23 16:38:33 +08:00

应该是代码里出现了大量的 for 循环，以及大量的__getitem__/__setitem__操作才慢的。numpy.ndarray 的随机读取性能确实不如 builtins.list 。因为`ndarray[i]` 需要 new 一个 PyLongObject 出来，而`builtins.list[i]`只需要 refcnt++。

楼主这个问题，如果用 numpy 的话，就得摆脱面向过程的思想，用函数式来做。numpy 的 broadcast 机制应该能实现楼主想要的功能。需要稍微学习一下。

如果不想用 numpy 的话，只能用预编译或者 jit 方案来加速了。既然楼主已经试过 numba 了，我推荐另一个工具 pythran. 性能与 numba 差不多, 但是比 numba 好用

jeeyong

2022-09-23 16:40:34 +08:00

@wcsjtu 我先尝试 numpy 的方案....再看预编译或者 jit 类的...
谢谢~

lmshl

2022-09-23 16:43:13 +08:00

前段时间处理过 DICOM 格式，不过我是转换为 Voxel ，量也不大不需要注重性能。

但我不精通 numpy ，理论上 numpy 应该不太能把常用 operator 解释为向量指令集，如果是简单求和之类的操作还好。

stein42

2022-09-23 16:47:40 +08:00

@jeeyong
这里是举例，实际应该从文件读取。
buffer = open(name, 'rb').read()

或者直接用 np.fromfile

思路很直接：
先得到 numpy 的数组。
再提取灰度和 alpha 的数组。
再拼接成二维数组，灰度重复 3 次。
最后 reshape 成需要的形状。

vicalloy

2022-09-23 17:10:21 +08:00

用 numpy 就要全程使用 numpy ，不然数据处理过程中会大量构造和销毁 python 对象，速度会很慢。

newmlp

2022-09-23 17:20:34 +08:00

直接 cpp 一把梭。。。用什么 python

faterazer

2022-09-23 17:23:08 +08:00

听你的描述，用 numpy 去搞，应该不会有什么性能问题呀，最后的数据维度应该是（ 4 通道数 * 高 * 宽），能否给一两个具体的 case 呢？就像 LeetCode 示例那样（数据不敏感的话可以放出一部分）。不过有几点要注意一下：

- 使用 numpy ，就多用向量化编程，避免 for 循环之类的（ numpy 做了向量操作优化，底层是经过优化的计算库，不用担心性能问题）。
- 如果真的要频繁的修改，建议先在 Python 的 list 对象上修改好，再转成 numpy

faterazer

2022-09-23 17:26:10 +08:00

@newmlp numpy 可不代表 python 哦

Cy86

2022-09-23 18:26:31 +08:00

楼主方便给一张 CT 影像么, 这样大伙可以本地跑测下

CamD

2022-09-23 18:37:43 +08:00 via iPhone

numpy ，或者 torch 移到 gpu 里操作。

jdhao

2022-09-23 19:45:35 +08:00 via Android

直接给个样例文件，让网友帮你看看，你这么说不好定位问题

chashao

2022-09-23 19:47:29 +08:00

求楼主发个 CT 影像，我试试 c++

iloveios

2022-09-23 19:53:19 +08:00 via iPhone

php 才是王道

Cy86

2022-09-23 19:56:27 +08:00

from numpy import maximum,uint8
from pydicom import dcmread
from PIL.Image import fromarray
from time import perf_counter

def convert_dcm_jpg():
im = dcmread('test.DCM')
rescaled_image = im.pixel_array.astype(float)
# rescaled_image = (maximum(rescaled_image,0)/rescaled_image.max())*255 # float pixels
final_image = uint8(rescaled_image) # integers pixels
final_image = fromarray(final_image)
return final_image
if __name__ == '__main__':
t = perf_counter()
image = convert_dcm_jpg()
image.save('test1.jpg')
print(F'coast:{perf_counter() - t:.8f}s')

Cy86

2022-09-23 19:57:40 +08:00

Cy86

2022-09-23 19:58:21 +08:00

我用的 http://www.rubomedical.com/dicom_files/dicom_viewer_Mrbrain.zip

paopjian

2022-09-23 20:13:38 +08:00

有考虑提供一张图片和原始代码?这样可能做优化更快

LnTrx

2022-09-23 20:21:13 +08:00

所有运算尽量用矩阵实现，然后直接用 pytorch 就行了

shinsekai

2022-09-23 20:22:59 +08:00

同样结构的代码，matlab 比 numpy 快 20 倍，楼主可以试试

Muniesa

2022-09-23 20:30:34 +08:00

我应该没理解错吧

inframe

2022-09-23 21:44:39 +08:00

简单，就是向量化思想方法
resharp 2 列 N 行，[v0,v1],v 是个列向量
然后拼装列向量为[v0,v0,v0,v1]

raw=np.random.randint(0,10, 10**8).reshape((-1,2))
v1=np.vstack([raw[:,0],raw[:,0],raw[:,1]])
output=v1.T

mizuBai

2022-09-23 21:45:04 +08:00 via iPhone

楼主可以试试写 cython 扩展，或者可以试试 f2py ，不过感觉 Fortran 这么上古的语言楼主会的概率不大（

Purelove

2022-09-23 23:29:53 +08:00

rust pyo3 了解一下

jeeyong

2022-09-24 00:29:51 +08:00

@Cy86
@jdhao
@chashao
@paopjian
老哥们, 辛苦啊...
我现在是直接用真实的患者 CT 片子写..
不能外传....

另外, 我看有人已经上 pydicom 了...
另外我请教一下, 我现在处理的 dicom 格式不是通用的 dicom 格式...
是排版好之后的 dcm 文件.
很多信息不能直接获取..最重要的 pydicom.pixel_array 和 pixelData 是没有的..
我是通过 pydicom
dcm = pydicom.read_file(filePath)
imageInfo = dcm.ReferencedImageBoxSequence[0].BasicGrayscaleImageSequence[0]
这种方式获取的 pixelData 的 bytes
有人能科普一下这个问题吗?

jeeyong

2022-09-24 00:57:20 +08:00

@Muniesa 老哥, 图完全看不清楚

lucays

2022-09-24 01:40:11 +08:00 via Android

dicom 是可以写入 pixel array 生成新的 dcm 文件的
只是模板没有的话需要摸索下
画图直接上 numpy+opencv

lucays

2022-09-24 01:43:26 +08:00 via Android

真实的 ct 片子也可以用 pydicom 读取然后抹掉 patient 信息生成新的不记名的不就完了

你这是完全从零开始啊，没点时间搞不定的，听起来一点 dicom 的基础建设都没

Nugine0

2022-09-24 01:49:09 +08:00 via Android

cython ，rust pyo3 ，c++，nodejs ，c 语言动态库，都可以试试。
如果只需要一次性转换，还可以开好点的临时云服务器，大力出奇迹。

Mayye

2022-09-24 01:55:44 +08:00

@Nugine0 感觉现在不是机器配置的问题而是代码效率的问题了云服务器应该帮助不大

Xs0ul

2022-09-24 02:01:40 +08:00

test_image = [208, 4, 208, 4, 208, 4, 100, 4, 100, 4, 100, 4]
test_image = np.array(test_image).reshape(2,3,2).astype(np.uint8)
image = PIL.Image.fromarray(test_image, mode='LA')
# if need RGBA
image = image.convert('RGBA')

pillow 本身就支持灰度 + alpha 的格式: https://pillow.readthedocs.io/en/stable/handbook/concepts.html

jeeyong

2022-09-24 02:57:36 +08:00

@lucays 第三天....T_T

jeeyong

2022-09-24 03:57:58 +08:00

@Xs0ul Hey~`
这方法效率很高..
但是有个问题....
他显示出来的图片...如果是 LA 模式, 很浅. 黑色的几乎看不清楚
我尝试用循环把 Alpha 通道值修改为 255 - 原值.
结果图像感觉噪点又很多
LA 模式:
https://imgur.com/QmmmcX9

转为 RGBA 模式
https://imgur.com/GPtekhA

期望得到的结果:
https://imgur.com/undefined

Xs0ul

2022-09-24 04:25:02 +08:00

@jeeyong #55 你最后一个链接是 undefined,我看不了期望的是怎么样的。但你开头给的那些数字，alpha 看起来非常的小？那显示出来就几乎是透明的。或许 208 是 alpha ，4 是灰度？具体 ct 影像的格式我不了解

milkpuff

2022-09-24 10:52:44 +08:00

上面说的向量化操作都是可以的。不要用循环索引的方式就行了。

lambdaq

2022-09-24 11:02:47 +08:00

python 要性能第一个原则就是只要你写 for 就已经输了。

ferstar

2022-09-24 11:38:04 +08:00 via Android

向量化加太极😏
https://docs.taichi-lang.org/blog/accelerate-python-code-100x

jeeyong

2022-09-24 13:17:51 +08:00

@Xs0ul 太着急了, 没粘好.
https://imgur.com/WWKXy9A

jeeyong

2022-09-24 13:24:10 +08:00

@Xs0ul 元素对调有什么好方法吗?
就是我想尝试把 208 当成 Alpha 通道, 4 当成灰度, 解析一次图像

Xs0ul

2022-09-24 13:31:11 +08:00

@jeeyong #61 test_image = np.array(test_image).reshape(2,3,2).astype(np.uint8)[:, :, ::-1]

encro

2022-09-24 14:32:39 +08:00

用 opencv 直接转应该效率会高一些？

encro

2022-09-24 14:34:27 +08:00

如果我没记错，转成 PNG 最慢的是保存！！！

jeeyong

2022-09-24 14:34:34 +08:00

@Xs0ul 嗯...208 是灰度, 4 是 alpha 我提高了 alpha 值, 255 - 原 alpha 值也还是不如我期望的那张图那么清晰.
奇怪了...

jeeyong

2022-09-24 14:35:21 +08:00

@encro 存 jpg 会报错, 就没去管它了, 现存成 png, 解决问题再优化

MaybeRichard

2022-09-24 14:47:56 +08:00

同样是做医学图像的，我们老师带我们做的项目只用 c++，CT 格式给的是 DICOM 和 MHD ，读取用的是 VtkMetaImageReader

jeeyong

2022-09-24 17:53:23 +08:00

结帖! 感谢中间提供帮助和建议的大佬们.

中间经历的几个需要反思的问题点.
1. 没用弄懂图片的格式相关知识的情况下, 盲目上手, 导致做了很多很多无效的计算.
2. 没有用好 numpy

附上源码
from numpy import maximum,uint8
from pydicom import dcmread
from PIL.Image import fromarray
from time import perf_counter
import numpy as np
import PIL
import time

start = time.time()
dcm = dcmread('dcm/aaa.dcm')
imageSources = dcm.ReferencedImageBoxSequence[0].BasicGrayscaleImageSequence[0]
imageData = imageSources.PixelData
# *** 下面这句是最重要和相对来说最耗时的
image = np.frombuffer(imageData, dtype=np.uint8).reshape(imageSources.Rows, imageSources.Columns, 2)

img = PIL.Image.fromarray(image, mode='LA')
img = img.convert('RGB')
img.save('test1.png')
print(time.time() - start)

总计用时 1.6S, numpy 果然强大.

@Xs0ul 谢谢! 另外, 我最后问的关于生成图片的问题, 是我弄错了, 我用两个完全不同部位的图片做的对比..
头部的 CT 结果轮廓很清晰, 腹部的, 看起来就不清楚了..纯粹是我粗心导致的错误.

xgdgsc

2022-09-25 10:05:50 +08:00 via Android

用 julia 写扩展 https://github.com/Suzhou-Tongyuan/jnumpy

laqow

2022-09-25 15:15:38 +08:00

转个 png 直接用 imagej 的宏不就完了，bioformat 把什么事情都做完了，就算 python 后面处理也是先把格式转统一了比如 tiff 之类的再做