Numba最近邻插值（CPU+ GPU + Z轴切块 + XYZ轴切块 + 多线程）

为什么使用 GPU 反而更慢了：

• （1）数据传输开销：GPU 的计算速度快，但数据在 CPU 和 GPU 之间传输时会有较大的开销。
• （2）任务并行性不高：GPU 适合大规模并行计算，如果任务的并行性不够高，比如 Z 轴切块后的任务处理，GPU 的潜力可能无法得到充分发挥。相比之下，CPU 在处理较小规模任务时可能表现得更有效率。
• （3）专用 GPU 内存不足，自动转共享 GPU 内存（时耗增加）。

最近邻插值（加速方法）

（1）scipy.ndimage.zoom

from scipy.ndimage import zoom
import time
import numpy as npif __name__ == "__main__":# （1）创建数组input_data = np.zeros((1024, 1024, int(1024 * 0.5))).astype(bool)input_data[50:600, 200:1000, 5:30] = 1# （2）插值计算start_time = time.time()#############################################################zoom_factor = [2, 2, 2]  # 指定插值因子   output_data = zoom(input_data, zoom_factor, order=0)#############################################################print("计算时间:", time.time() - start_time)print("获取非零元素的数量（输入）:", np.count_nonzero(input_data))print("获取非零元素的数量（输出）:", np.count_nonzero(output_data))print("原始数据:", input_data.shape)print("插值结果:", output_data.shape)"""
##### 插值两倍 #####
计算时间: 21.160449504852295
获取非零元素的数量（输入）: 11000000
获取非零元素的数量（输出）: 88000000
原始数据: (1024, 1024, 512)
插值结果: (2048, 2048, 1024)##### 插值四倍 #####
计算时间: 172.16581082344055
获取非零元素的数量（输入）: 11000000
获取非零元素的数量（输出）: 704760200
原始数据: (1024, 1024, 512)
插值结果: (4096, 4096, 2048)
"""

（2）Numba-CPU加速

import numba
import numpy as np@numba.jit()
def nearest_neighbor_interpolation(input_data, scale_factors, output_data):"""最近邻插值算法输入参数：3D图像 + 插值因子 + 预分配的输出数据"""# 对目标数组中的每个点进行插值for z in range(output_data.shape[0]):zz = min(round(z / scale_factors[0]), input_data.shape[0] - 1)  # round四舍五入可能会超出数据边界的值for y in range(output_data.shape[1]):yy = min(round(y / scale_factors[1]), input_data.shape[1] - 1)  # round四舍五入可能会超出数据边界的值for x in range(output_data.shape[2]):# 计算在原始数据中的对应位置，并限制在原始数据范围内xx = min(round(x / scale_factors[2]), input_data.shape[2] - 1)  # round四舍五入可能会超出数据边界的值output_data[z, y, x] = input_data[zz, yy, xx]  # 最近邻插值return output_dataif __name__ == "__main__":# （1）创建数组input_data = np.zeros((1024, 1024, int(1024*0.5)), dtype=bool)  # 创建3D数组input_data[50:600, 200:1000, 5:30] = 1# （2）插值计算import timestart_time = time.time()############################################################## 计算目标形状并创建目标数组output_shape = np.round(np.array(input_data.shape) * np.array(scale_factors)).astype(int)output_data = np.zeros(output_shape, dtype=input_data.dtype)scale_factors = [4, 4, 4]  # 指定插值因子nearest_neighbor_interpolation(input_data, scale_factors, output_data)  # 执行3D最近邻插值#############################################################print("计算时间:", time.time() - start_time)print("获取非零元素的数量（输入）:", np.count_nonzero(input_data))print("获取非零元素的数量（输出）:", np.count_nonzero(output_data))print("原始数据:", input_data.shape)print("插值结果:", output_data.shape)"""
##### 插值两倍 #####
计算时间: 7.694857835769653
获取非零元素的数量（输入）: 11000000
获取非零元素的数量（输出）: 86240000
原始数据: (1024, 1024, 512)
插值结果: (2048, 2048, 1024)##### 插值四倍 #####
计算时间: 56.58441758155823
获取非零元素的数量（输入）: 11000000
获取非零元素的数量（输出）: 696960000
原始数据: (1024, 1024, 512)
插值结果: (4096, 4096, 2048)
"""

（3）Numba-GPU加速

import numpy as np
from numba import cuda@cuda.jit
def nearest_neighbor_interpolation_gpu(input_data, output_data, scale_factors):"""最近邻插值的CUDA版本"""z, y, x = cuda.grid(3)  # 获取3D网格中线程的位置 (z, y, x)if z < output_data.shape[0] and y < output_data.shape[1] and x < output_data.shape[2]:zz = min(round(z / scale_factors[0]), input_data.shape[0] - 1)yy = min(round(y / scale_factors[1]), input_data.shape[1] - 1)xx = min(round(x / scale_factors[2]), input_data.shape[2] - 1)output_data[z, y, x] = input_data[zz, yy, xx]def gpu_nearest_neighbor_interpolation(input_data, scale_factors):# 创建目标数组的形状output_shape = np.round(input_data.shape * scale_factors).astype(int)output_data = np.zeros(output_shape, dtype=input_data.dtype)# 将输入数据从CPU上传到GPUinput_data_gpu = cuda.to_device(input_data)output_data_gpu = cuda.to_device(output_data)######################################################################## 定义线程和块的大小threads_per_block = (16, 16, 4)blocks_per_grid = ((output_shape[0] + threads_per_block[0] - 1) // threads_per_block[0],(output_shape[1] + threads_per_block[1] - 1) // threads_per_block[1],(output_shape[2] + threads_per_block[2] - 1) // threads_per_block[2])# 执行CUDA核函数nearest_neighbor_interpolation_gpu[blocks_per_grid, threads_per_block](input_data_gpu, output_data_gpu, scale_factors)######################################################################## 将结果从GPU下载回CPUoutput_data = output_data_gpu.copy_to_host()return output_dataif __name__ == "__main__":# （1）创建3D数组input_data = np.zeros((1024, 1024, int(1024 * 0.5))).astype(bool)input_data[50:600, 200:1000, 5:30] = 1# （2）执行插值import timestart_time = time.time()#######################################################################scale_factors = np.array([4, 4, 4])  # 指定插值因子output_data = gpu_nearest_neighbor_interpolation(input_data, scale_factors)#######################################################################print("计算时间:", time.time() - start_time)print("获取非零元素的数量（输入）:", np.count_nonzero(input_data))print("获取非零元素的数量（输出）:", np.count_nonzero(output_data))print("原始数据:", input_data.shape)print("插值结果:", output_data.shape)"""
##### 插值两倍 #####
计算时间: 3.6601688861846924
获取非零元素的数量（输入）: 11000000
获取非零元素的数量（输出）: 86240000
原始数据: (1024, 1024, 512)
插值结果: (2048, 2048, 1024)##### 插值四倍 #####
计算时间: 122.51327633857727
获取非零元素的数量（输入）: 11000000
获取非零元素的数量（输出）: 696960000
原始数据: (1024, 1024, 512)
插值结果: (4096, 4096, 2048)
"""

（4）Numba-CPU加速（Z轴切块）

import numba
import numpy as np@numba.jit()
def nearest_neighbor_interpolation(input_data, scale_factors, output_block):# 对目标数组中的每个点进行插值for z in range(output_block.shape[0]):zz = min(round(z / scale_factors[0]), input_data.shape[0] - 1)for y in range(output_block.shape[1]):yy = min(round(y / scale_factors[1]), input_data.shape[1] - 1)for x in range(output_block.shape[2]):# 计算在原始数据中的对应位置，并限制在原始数据范围内xx = min(round(x / scale_factors[2]), input_data.shape[2] - 1)output_block[z, y, x] = input_data[zz, yy, xx]  # 最近邻插值return output_blockdef segment_blocks(input_data, block_size, scale_factors):num_blocks = int(np.ceil(input_data.shape[0] / block_size))print("切块数量 =", num_blocks)zoomed_blocks = []for i in range(num_blocks):start_idx = i * block_sizeend_idx = min((i + 1) * block_size, input_data.shape[0])block = input_data[start_idx:end_idx, :, :]##########################################################################output_shape = np.round(np.array(block.shape) * scale_factors).astype(int)  # 计算目标形状output_block = np.zeros(output_shape, dtype=block.dtype)  # 创建目标形状的数组zoomed_block = nearest_neighbor_interpolation(block, scale_factors, output_block)# zoomed_block = scipy.ndimage.zoom(block, zoom_factor, order=1)##########################################################################zoomed_blocks.append(zoomed_block)output_data = np.concatenate(zoomed_blocks, axis=0)  # 沿着第一个轴进行拼接得到合并结果return output_dataif __name__ == "__main__":input_data = np.zeros((1024, 1024, int(1024 * 0.5))).astype(bool)  # 创建一个示例的3D数组input_data[50:600, 200:1000, 5:30] = 1import timestart_time = time.time()###############################################################scale_factors = np.array([4, 4, 4])  # 指定插值因子block_size = 100  # 切割小块的Z轴尺度，会自动分配并处理越界问题（切块数量 = input_data[0] / block_size）output_data = segment_blocks(input_data, block_size, scale_factors)  # 分块插值##########################################################################print("计算时间:", time.time() - start_time)print("获取非零元素的数量（输入）:", np.count_nonzero(input_data))print("获取非零元素的数量（输出）:", np.count_nonzero(output_data))print("原始数据:", input_data.shape)print("插值结果:", output_data.shape)"""
##### 插值两倍 #####
切块数量 = 11
总共耗时: 6.803957223892212
Non-zero Count: 11000000
Non-zero Count: 86318400
原始数据: (1024, 1024, 512)
插值结果: (2048, 2048, 1024)##### 插值四倍 #####
切块数量 = 11
计算时间: 52.44191575050354
获取非零元素的数量（输入）: 11000000
获取非零元素的数量（输出）: 697593600
原始数据: (1024, 1024, 512)
插值结果: (4096, 4096, 2048)
"""

（5）Numba-CPU加速（XYZ轴切块）

import numba
import numpy as np@numba.jit()
def nearest_neighbor_interpolation(input_data, scale_factors, output_block):for z in range(output_block.shape[0]):zz = min(round(z / scale_factors[0]), input_data.shape[0] - 1)for y in range(output_block.shape[1]):yy = min(round(y / scale_factors[1]), input_data.shape[1] - 1)for x in range(output_block.shape[2]):xx = min(round(x / scale_factors[2]), input_data.shape[2] - 1)output_block[z, y, x] = input_data[zz, yy, xx]return output_blockdef segment_blocks(input_data, block_sizes, scale_factors):num_blocks_z = int(np.ceil(input_data.shape[0] / block_sizes[0]))num_blocks_y = int(np.ceil(input_data.shape[1] / block_sizes[1]))num_blocks_x = int(np.ceil(input_data.shape[2] / block_sizes[2]))zoomed_blocks_zz = []for z in range(num_blocks_z):start_z = z * block_sizes[0]end_z = min((z + 1) * block_sizes[0], input_data.shape[0])zoomed_blocks_yy = []for y in range(num_blocks_y):start_y = y * block_sizes[1]end_y = min((y + 1) * block_sizes[1], input_data.shape[1])zoomed_blocks_xx = []for x in range(num_blocks_x):start_x = x * block_sizes[2]end_x = min((x + 1) * block_sizes[2], input_data.shape[2])block = input_data[start_z:end_z, start_y:end_y, start_x:end_x]print(f"{num_blocks_z, num_blocks_y, num_blocks_x}/{z + 1, y + 1, x + 1}", "block.shape =", block.shape)##########################################################################output_shape = np.round(block.shape * scale_factors).astype(int)output_block = np.zeros(output_shape, dtype=block.dtype)zoomed_block_x = nearest_neighbor_interpolation(block, scale_factors, output_block)##########################################################################zoomed_blocks_xx.append(zoomed_block_x)zoomed_block_y = np.concatenate(zoomed_blocks_xx, axis=2)  # 沿着X轴拼接zoomed_blocks_yy.append(zoomed_block_y)zoomed_block_z = np.concatenate(zoomed_blocks_yy, axis=1)  # 沿着Y轴拼接zoomed_blocks_zz.append(zoomed_block_z)zoomed_block = np.concatenate(zoomed_blocks_zz, axis=0)  # 沿着Z轴拼接print("切块数量 =", num_blocks_z * num_blocks_y * num_blocks_x)return zoomed_blockif __name__ == "__main__":input_data = np.zeros((1024, 1024, int(1024 * 0.5))).astype(bool)input_data[50:600, 200:1000, 5:30] = 1import timestart_time = time.time()##########################################################################scale_factors = np.array([4, 4, 4])block_sizes = (100, 100, 100)  # 小块的尺度，沿着z轴、y轴、x轴zoomed_block = segment_blocks(input_data, block_sizes, scale_factors)##########################################################################print("计算时间:", time.time() - start_time)print("获取非零元素的数量（输入）:", np.count_nonzero(input_data))print("获取非零元素的数量（输出）:", np.count_nonzero(zoomed_block))print("原始数据:", input_data.shape)print("插值结果:", zoomed_block.shape)"""
##### 插值两倍 #####
切块数量 = 726
计算时间: 9.834398746490479
获取非零元素的数量（输入）: 11000000
获取非零元素的数量（输出）: 86318400
原始数据: (1024, 1024, 512)
插值结果: (2048, 2048, 1024)##### 插值四倍 #####
切块数量 = 726
计算时间: 72.6902551651001
获取非零元素的数量（输入）: 11000000
获取非零元素的数量（输出）: 697593600
原始数据: (1024, 1024, 512)
插值结果: (4096, 4096, 2048)
"""

（6）Numba-CPU加速（XYZ轴切块）+ 多线程

import numba
import numpy as np
import concurrent.futures@numba.jit()
def nearest_neighbor_interpolation(input_data, scale_factors, output_block):for z in range(output_block.shape[0]):zz = min(round(z / scale_factors[0]), input_data.shape[0] - 1)for y in range(output_block.shape[1]):yy = min(round(y / scale_factors[1]), input_data.shape[1] - 1)for x in range(output_block.shape[2]):xx = min(round(x / scale_factors[2]), input_data.shape[2] - 1)output_block[z, y, x] = input_data[zz, yy, xx]return output_blockdef process_block(block, scale_factors):output_shape = np.round(block.shape * scale_factors).astype(int)output_block = np.zeros(output_shape, dtype=block.dtype)zoomed_block = nearest_neighbor_interpolation(block, scale_factors, output_block)return zoomed_blockdef segment_blocks(input_data, block_sizes, scale_factors):num_blocks_z = int(np.ceil(input_data.shape[0] / block_sizes[0]))num_blocks_y = int(np.ceil(input_data.shape[1] / block_sizes[1]))num_blocks_x = int(np.ceil(input_data.shape[2] / block_sizes[2]))zoomed_blocks = []with concurrent.futures.ThreadPoolExecutor() as executor:futures = []for z in range(num_blocks_z):start_z = z * block_sizes[0]end_z = min((z + 1) * block_sizes[0], input_data.shape[0])for y in range(num_blocks_y):start_y = y * block_sizes[1]end_y = min((y + 1) * block_sizes[1], input_data.shape[1])for x in range(num_blocks_x):start_x = x * block_sizes[2]end_x = min((x + 1) * block_sizes[2], input_data.shape[2])##########################################################################block = input_data[start_z:end_z, start_y:end_y, start_x:end_x]futures.append(executor.submit(process_block, block, scale_factors))########################################################################### for future in concurrent.futures.as_completed(futures):#     zoomed_blocks.append(future.result())print("切块数量 =", num_blocks_z * num_blocks_y * num_blocks_x)return zoomed_blocksif __name__ == "__main__":input_data = np.zeros((1024, 1024, int(1024 * 0.5))).astype(bool)input_data[50:600, 200:1000, 5:30] = 1import timestart_time = time.time()##########################################################################scale_factors = np.array([4, 4, 4])block_sizes = (100, 100, 100)  # 小块的尺度，沿着z轴、y轴、x轴output_data = segment_blocks(input_data, block_sizes, scale_factors)##########################################################################print("计算时间:", time.time() - start_time)print("获取非零元素的数量（输入）:", np.count_nonzero(input_data))# print("获取非零元素的数量（输出）:", np.count_nonzero(output_data))print("原始数据:", input_data.shape)# print("插值结果:", output_data.shape)"""
##### 插值两倍 #####
切块数量 = 726
计算时间: 6.778625011444092
获取非零元素的数量（输入）: 11000000
原始数据: (1024, 1024, 512)##### 插值四倍 #####
切块数量 = 726
计算时间: 50.20111918449402
获取非零元素的数量（输入）: 11000000
原始数据: (1024, 1024, 512)
"""