调用C函数处理Numpy矩阵时的边界填充问题排查

我希望通过外部C函数对图像执行Numpy核运算，目标是仅在填充范围内遍历图像，避免核越界。调试代码中，Python脚本先生成全1矩阵，将其及维度传入C函数；期望矩阵边界（填充区）设为9，可计算区域设为8。最终目标是对矩阵每个元素进行3×3邻域计算，但当前运行结果与预期不符。

预期结果 #

,0,1,2,3
0,9,9,9,9
1,9,8,8,9
2,9,8,8,9
3,9,9,9,9

实际结果 #

,0,1,2,3
0,9,9,9,9
1,8,8,8,8
2,8,8,8,8
3,8,1,1,1

相关代码 #

Shell脚本（debug_scriptsh.sh） #

a=debug_scriptc && cc -fPIC -shared -o $a.so $a.c && python debug_scriptpy.py && open debug_table.csv

Python脚本（debug_scriptpy.py） #

import ctypes
import pandas as pd
import numpy as np
import sys
import pdb

DLL_NAME = "./debug_scriptc.{:s}".format("dll" if sys.platform[:3].lower() == "win" else "so")

def np_mat_type(rows, cols, element_type=float):
  return np.ctypeslib.ndpointer(dtype=element_type, shape=(rows, cols), flags="C_CONTIGUOUS")

rows0 = 4
cols0 = 4
kernel_size = 3

dll = ctypes.CDLL(DLL_NAME)
matrix_func = dll.matrixFunc4matin
matrix_func.argtypes = (
    ctypes.c_size_t,
    np_mat_type(rows0, cols0), 
    ctypes.c_size_t, 
    ctypes.c_size_t
    )
matrix_func.restype = np_mat_type(rows0, cols0)

dealloc_array = dll.deallocArray
dealloc_array.argtypes = (np_mat_type(rows0, cols0),)
dealloc_array.restype = None
u_rms_array = np.zeros((rows0, cols0))

velocity_image_ones = np.empty((rows0, cols0))
velocity_image_ones[:] = 1

u_rms_array = matrix_func(
    kernel_size, 
    velocity_image_ones,
    rows0, 
    cols0 
    )
debug_DF_u_rms_array_C = pd.DataFrame(u_rms_array.astype(int))
debug_DF_u_rms_array_C.to_csv("debug_table.csv")
dealloc_array(u_rms_array)

C脚本（debug_scriptc.c） #

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#if defined(_WIN32)
#  define DLL00_EXPORT_API __declspec(dllexport)
#else
#  define DLL00_EXPORT_API
#endif

#define ELEMENT_TYPE double

#if defined(__cplusplus)
extern "C" {
#endif

DLL00_EXPORT_API ELEMENT_TYPE* matrixFunc4matin(
  size_t kernel_size,
  const ELEMENT_TYPE *pvelocity_image_ones, 
  size_t rows0, 
  size_t cols0
  );
DLL00_EXPORT_API void deallocArray(ELEMENT_TYPE *pMat);

#if defined(__cplusplus)
}
#endif

ELEMENT_TYPE* matrixFunc4matin(
  size_t kernel_size,
  const ELEMENT_TYPE *pvelocity_image_ones, 
  size_t rows0, 
  size_t cols0 
  )

{
  size_t matSize = sizeof(ELEMENT_TYPE) * cols0 * rows0;
  ELEMENT_TYPE *pRet = (ELEMENT_TYPE*)(malloc(matSize));

  memcpy(pRet, pvelocity_image_ones, matSize);  

  int start_kern_search = kernel_size - 2;
  int end_kern_search_rows = rows0 - (kernel_size-1)/2;
  int end_kern_search_cols = cols0 - (kernel_size-1)/2;

  for (int i_velmap = 0; i_velmap < rows0; i_velmap++){  
    for (int j_velmap = 0; j_velmap < cols0; j_velmap++){  
      if (i_velmap < start_kern_search | i_velmap > end_kern_search_rows | j_velmap < start_kern_search | j_velmap > end_kern_search_cols){
        pRet[i_velmap * cols0 + j_velmap] = 9;
        continue;
      }
      else
      {
        pRet[i_velmap * end_kern_search_cols + j_velmap] = 8;
      }
    }
  }
  return pRet;
}


void deallocArray(ELEMENT_TYPE *pMat)
{
    if (pMat)
        free(pMat);
}

错误分析与修正 #

核心错误点 #

1. 边界判断逻辑错误：

   • 原代码使用|（按位或）替代||（逻辑或），导致边界条件判断失效。按位或会对条件的二进制位进行运算，而非判断逻辑真假，使得本该标记为边界的区域错误进入可计算分支。
   • 边界阈值判断错误：结束索引条件使用>而非>=，对于从0开始的索引，end_kern_search_rows = 4-1=3，索引3属于边界区域，应该被设为9。

2. 内存索引计算错误：

   • 可计算区域赋值时，错误使用end_kern_search_cols作为列数计算内存偏移，导致内存越界和错误赋值，正确的应该使用原始矩阵列数cols0。

修正后的C函数关键代码 #

ELEMENT_TYPE* matrixFunc4matin(
  size_t kernel_size,
  const ELEMENT_TYPE *pvelocity_image_ones, 
  size_t rows0, 
  size_t cols0 
  )

{
  size_t matSize = sizeof(ELEMENT_TYPE) * cols0 * rows0;
  ELEMENT_TYPE *pRet = (ELEMENT_TYPE*)(malloc(matSize));

  memcpy(pRet, pvelocity_image_ones, matSize);  

  int pad = (kernel_size - 1) / 2; // 统一计算填充宽度，适配不同核大小
  int start_idx = pad;
  int end_idx_rows = rows0 - pad;
  int end_idx_cols = cols0 - pad;

  for (int i_velmap = 0; i_velmap < rows0; i_velmap++){  
    for (int j_velmap = 0; j_velmap < cols0; j_velmap++){  
      // 使用逻辑或||，修正边界判断条件
      if (i_velmap < start_idx || i_velmap >= end_idx_rows || j_velmap < start_idx || j_velmap >= end_idx_cols){
        pRet[i_velmap * cols0 + j_velmap] = 9;
      }
      else
      {
        // 使用正确的原始列数计算内存偏移
        pRet[i_velmap * cols0 + j_velmap] = 8;
      }
    }
  }
  return pRet;
}

修正说明 #

• 用pad变量统一计算填充宽度，逻辑更清晰，支持5×5等其他尺寸的核。
• 边界判断改用||逻辑或，确保任一边界条件满足时标记为填充区。
• 结束索引判断改为>=，准确覆盖边界区域。
• 赋值时使用原始矩阵列数cols0计算偏移，避免内存越界和错误赋值。

内容的提问来源于stack exchange，提问作者ecjb