Python：带约束的3D数组转2D结果的中值滤波技术问询

Alright, let's dive into solving this problem—your 3D array is massive, so memory efficiency and speed are going to be non-negotiable here. Let's recap your requirement clearly: for every (x,y) position in your 12000×10000 spatial grid, you need the median of all non-zero values from the 100×100 spatial neighborhood centered at (x,y), across all 750 z-slices. Here's how to approach this effectively:

A 750×12000×10000 array stored as float32 takes up ~360GB of memory—way too big to load all at once. So we'll use memory-mapped arrays to work with slices of the data without loading everything into RAM.

• Use numpy.memmap to read your input 3D array and write the output 2D array directly to disk, avoiding memory overload.

We'll cover three approaches, from basic to high-performance, depending on your hardware and speed needs.

Option 1: Basic NumPy (For Smaller Test Runs) #

This is straightforward but slow for full-scale data—great for testing logic on a smaller subset first.

import numpy as np

# Configure paths and shapes
input_path = "your_3d_array.npy"
output_path = "output_median.npy"
input_shape = (750, 12000, 10000)
output_shape = (12000, 10000)
win_size = 100
half_win = win_size // 2

# Open memory-mapped arrays
arr_3d = np.memmap(input_path, dtype=np.float32, mode='r', shape=input_shape)
output_2d = np.memmap(output_path, dtype=np.float32, mode='w+', shape=output_shape)

# Process internal grid (avoid boundary edges first)
for x in range(half_win, input_shape[1] - half_win):
    for y in range(half_win, input_shape[2] - half_win):
        # Extract the 100x100 spatial window across all z-slices
        window = arr_3d[:, x-half_win:x+half_win, y-half_win:y+half_win]
        # Filter out zero values
        non_zero_vals = window[window != 0]
        # Calculate median (fallback to 0 if no non-zero values exist)
        output_2d[x, y] = np.median(non_zero_vals) if len(non_zero_vals) > 0 else 0

# Fill boundary edges with nearest valid values (adjust this logic to your needs)
output_2d[:half_win, :] = output_2d[half_win, :]
output_2d[input_shape[1]-half_win:, :] = output_2d[input_shape[1]-half_win-1, :]
output_2d[:, :half_win] = output_2d[:, half_win]
output_2d[:, input_shape[2]-half_win:] = output_2d[:, input_shape[2]-half_win-1]

# Ensure all data is written to disk
del output_2d

Option 2: Parallel Processing (Speed Up CPU Execution) #

Since each window's calculation is independent, we can split the grid into blocks and process them in parallel with multiple CPU cores.

import numpy as np
from multiprocessing import Pool

def process_block(block_args):
    arr_3d, output_2d, x_start, x_end, y_start, y_end, win_size, half_win = block_args
    input_h, input_w = arr_3d.shape[1], arr_3d.shape[2]
    
    for x in range(x_start, x_end):
        for y in range(y_start, y_end):
            # Skip boundary positions to avoid index errors
            if x < half_win or x >= input_h - half_win or y < half_win or y >= input_w - half_win:
                continue
            window = arr_3d[:, x-half_win:x+half_win, y-half_win:y+half_win]
            non_zero_vals = window[window != 0]
            output_2d[x, y] = np.median(non_zero_vals) if len(non_zero_vals) > 0 else 0
    return

# Initialize memory-mapped arrays
input_path = "your_3d_array.npy"
output_path = "output_median.npy"
input_shape = (750, 12000, 10000)
output_shape = (12000, 10000)
win_size = 100
half_win = win_size // 2

arr_3d = np.memmap(input_path, dtype=np.float32, mode='r', shape=input_shape)
output_2d = np.memmap(output_path, dtype=np.float32, mode='w+', shape=output_shape)

# Split grid into manageable blocks (adjust block size based on your RAM)
x_blocks = np.linspace(0, input_shape[1], 13, dtype=int)  # 12 blocks of 1000 each
y_blocks = np.linspace(0, input_shape[2], 11, dtype=int)  # 10 blocks of 1000 each

# Build task list for parallel processing
tasks = []
for i in range(len(x_blocks)-1):
    for j in range(len(y_blocks)-1):
        tasks.append((
            arr_3d, output_2d,
            x_blocks[i], x_blocks[i+1],
            y_blocks[j], y_blocks[j+1],
            win_size, half_win
        ))

# Run parallel processing (adjust processes to match your CPU core count)
with Pool(processes=8) as pool:
    pool.map(process_block, tasks)

# Fill boundary edges
output_2d[:half_win, :] = output_2d[half_win, :]
output_2d[input_shape[1]-half_win:, :] = output_2d[input_shape[1]-half_win-1, :]
output_2d[:, :half_win] = output_2d[:, half_win]
output_2d[:, input_shape[2]-half_win:] = output_2d[:, input_shape[2]-half_win-1]

del output_2d

Option 3: GPU Acceleration (For Maximum Speed) #

If you have access to an NVIDIA GPU, use CuPy (a GPU-accelerated NumPy drop-in replacement) to process the data orders of magnitude faster.

import cupy as cp

# Load data to GPU (use memmap if even GPU memory is insufficient)
arr_3d_gpu = cp.load("your_3d_array.npy")
win_size = 100
half_win = win_size // 2

# Create sliding window view across spatial dimensions
# Output shape: (12000-99, 10000-99, 750, 100, 100)
window_view = cp.lib.stride_tricks.sliding_window_view(
    arr_3d_gpu, 
    window_shape=(750, win_size, win_size),
    axis=(0, 1, 2)
)

# Flatten each window into a 1D array of values
window_flat = window_view.reshape(window_view.shape[0], window_view.shape[1], -1)

# Define function to compute median of non-zero values in a window
def compute_non_zero_median(arr):
    non_zero = arr[arr != 0]
    return cp.median(non_zero) if len(non_zero) > 0 else 0

# Apply median calculation across all windows
output_gpu = cp.apply_along_axis(compute_non_zero_median, axis=2, arr=window_flat)

# Pad edges to match original 12000×10000 shape (using edge values)
output_2d = cp.pad(output_gpu, pad_width=half_win, mode='edge').get()

# Save result to disk
cp.save("output_median.npy", output_2d)

• Stick to smaller data types: Use float32 instead of float64 if your precision allows—cuts memory usage in half.
• Adjust block size: For parallel processing, tweak block sizes to balance memory usage and CPU utilization.
• Custom boundary logic: Replace edge padding with zero/NaN or other values if that's more appropriate for your use case.

内容的提问来源于stack exchange，提问作者Alisson do Carmo