python调用c之3D三线性插值 重采样

python调用c之3D三线性插值 重采样

实际的参考是根据cityscapsscripts修改的cityscapesscripts/cityscapesscripts/evaluation at master · mcordts/cityscapesscripts · GitHub

1、首先编写的C语言脚本

volume_resample_impl.c

我挺喜欢这个地方加上_impl的，这样方便和打包调用的名称混淆，其他涉及到的细节我就不写了。

void volume_resample_nearast(
			const unsigned short*  f_oriImg_p ,
			unsigned short *       f_resImg_p ,
			const unsigned int*    f_size_p ,
			const float*           f_origin_p ,
			const float*           f_spacing_p ,
			const unsigned int*    f_tsize_p,
			const float*           f_torigin_p,
			const float*           f_tspacing_p
	)
	{
	const unsigned int out_nvox = f_tsize_p[0] * f_tsize_p[1] * f_tsize_p[2];
    for (unsigned int out_idx = 0; out_idx < out_nvox; out_idx = out_idx +1)
    {
        unsigned int*  out_ijk[3];
        COORDS_FROM_INDEX(out_ijk, out_idx, f_tsize_p);
        float*  xyz[3];
        volume2patient(out_ijk, f_torigin_p, f_tspacing_p, xyz);
        unsigned int*  in_ijk[3];
        patient2continusvolume(xyz, f_origin_p, f_spacing_p, in_ijk );

        if(is_inside(in_ijk,f_size_p))
        {
            unsigned int in_idx = INDEX_FROM_COORDS(in_ijk, f_size_p);
            f_resImg_p[out_idx] = f_oriImg_p[in_idx] ;

        }
      }
	}

2、编写pxy，这个类似于dll中export文件，之后这个pxy文件会被变异成pyd，

# cython methods to speed-up resample
from enum import Enum
import numpy as np
cimport cython
cimport numpy as np
import ctypes

np.import_array()
# ctypedef   signed short int16_t

cdef extern from "volume_resample_impl.c":
	void volume_resample_linear(
			const signed short*  f_oriImg_p ,
			signed short *       f_resImg_p ,
			const unsigned int*    f_size_p ,
			const float*           f_origin_p ,
			const float*           f_spacing_p ,
			const unsigned int*    f_tsize_p,
			const float*           f_torigin_p,
			const float*           f_tspacing_p,
			const signed short   f_defValue_us
	)
	void mask_resample_linear(
			const unsigned  char*  f_oriImg_p ,
			unsigned  char *       f_resImg_p ,
			const unsigned int*    f_size_p ,
			const float*           f_origin_p ,
			const float*           f_spacing_p ,
			const unsigned int*    f_tsize_p,
			const float*           f_torigin_p,
			const float*           f_tspacing_p,
			const unsigned  char  f_defValue_us
	)
cdef tonumpymask(unsigned char * data, unsigned int* size):
	if not (data and size ): raise ValueError
	return np.PyArray_SimpleNewFromData(3, [size[0], size[1], size[2]], np.NPY_UINT8, data)


cdef tonumpyarray(signed short * data, unsigned int* size):
	if not (data and size ): raise ValueError
	return np.PyArray_SimpleNewFromData(3, [size[0], size[1], size[2]], np.NPY_INT16, data)

@cython.boundscheck(False)
def cvolume_resample_linear( np.ndarray[np.int16_t , ndim=3] oriArr ,
							 np.ndarray[np.int16_t , ndim=3] resArr ,
							 np.ndarray[np.uint32_t,   ndim=1] sizeArr ,
							 np.ndarray[np.float32_t,   ndim=1] originArr ,
							 np.ndarray[np.float32_t,   ndim=1] spacingArr ,
							 np.ndarray[np.uint32_t,   ndim=1] tsizeArr ,
							 np.ndarray[np.float32_t,   ndim=1] toriginArr ,
							 np.ndarray[np.float32_t,   ndim=1] tspacingArr,
							 np.int16_t defValue
							 ):
	cdef np.ndarray[np.int16_t    , ndim=3, mode="c"] oriArr_c
	cdef np.ndarray[np.int16_t    , ndim=3, mode="c"] resArr_c
	cdef np.ndarray[np.uint32_t     , ndim=1, mode="c"] sizeArr_c
	cdef np.ndarray[np.float32_t     , ndim=1, mode="c"] originArr_c
	cdef np.ndarray[np.float32_t     , ndim=1, mode="c"] spacingArr_c
	cdef np.ndarray[np.uint32_t     , ndim=1, mode="c"] tsizeArr_c
	cdef np.ndarray[np.float32_t     , ndim=1, mode="c"] toriginArr_c
	cdef np.ndarray[np.float32_t     , ndim=1, mode="c"] tspacingArr_c
	cdef np.int16_t defValue_c

	oriArr_c       = np.ascontiguousarray(oriArr ,      dtype=np.int16)
	resArr_c       = np.ascontiguousarray(resArr ,      dtype=np.int16)
	sizeArr_c      = np.ascontiguousarray(sizeArr ,     dtype=np.uint32)
	originArr_c    = np.ascontiguousarray(originArr ,   dtype=np.float32)
	spacingArr_c   = np.ascontiguousarray(spacingArr ,  dtype=np.float32)
	tsizeArr_c     = np.ascontiguousarray(tsizeArr ,    dtype=np.uint32)
	toriginArr_c   = np.ascontiguousarray(toriginArr ,  dtype=np.float32)
	tspacingArr_c  = np.ascontiguousarray(tspacingArr , dtype=np.float32)
	defValue_c = np.int16(defValue)

	volume_resample_linear( &oriArr_c[0,0,0] , &resArr_c[0,0,0] , &sizeArr_c[0] , &originArr_c[0] , &spacingArr_c[0] , &tsizeArr_c[0], &toriginArr_c[0], &tspacingArr_c[0], defValue_c )
	resArr = np.ascontiguousarray(tonumpyarray(&resArr_c[0,0,0], &tsizeArr_c[0]))
	return np.copy(resArr)

@cython.boundscheck(False)
def cmask_resample_linear( np.ndarray[np.uint8_t , ndim=3] oriArr ,
							 np.ndarray[np.uint8_t , ndim=3] resArr ,
							 np.ndarray[np.uint32_t,   ndim=1] sizeArr ,
							 np.ndarray[np.float32_t,   ndim=1] originArr ,
							 np.ndarray[np.float32_t,   ndim=1] spacingArr ,
							 np.ndarray[np.uint32_t,   ndim=1] tsizeArr ,
							 np.ndarray[np.float32_t,   ndim=1] toriginArr ,
							 np.ndarray[np.float32_t,   ndim=1] tspacingArr,
							 np.uint8_t defValue
							 ):
	cdef np.ndarray[np.uint8_t    , ndim=3, mode="c"] oriArr_c
	cdef np.ndarray[np.uint8_t    , ndim=3, mode="c"] resArr_c
	cdef np.ndarray[np.uint32_t     , ndim=1, mode="c"] sizeArr_c
	cdef np.ndarray[np.float32_t     , ndim=1, mode="c"] originArr_c
	cdef np.ndarray[np.float32_t     , ndim=1, mode="c"] spacingArr_c
	cdef np.ndarray[np.uint32_t     , ndim=1, mode="c"] tsizeArr_c
	cdef np.ndarray[np.float32_t     , ndim=1, mode="c"] toriginArr_c
	cdef np.ndarray[np.float32_t     , ndim=1, mode="c"] tspacingArr_c
	cdef np.uint8_t defValue_c

	oriArr_c       = np.ascontiguousarray(oriArr ,      dtype=np.uint8)
	resArr_c       = np.ascontiguousarray(resArr ,      dtype=np.uint8)
	sizeArr_c      = np.ascontiguousarray(sizeArr ,     dtype=np.uint32)
	originArr_c    = np.ascontiguousarray(originArr ,   dtype=np.float32)
	spacingArr_c   = np.ascontiguousarray(spacingArr ,  dtype=np.float32)
	tsizeArr_c     = np.ascontiguousarray(tsizeArr ,    dtype=np.uint32)
	toriginArr_c   = np.ascontiguousarray(toriginArr ,  dtype=np.float32)
	tspacingArr_c  = np.ascontiguousarray(tspacingArr , dtype=np.float32)
	defValue_c = np.uint8(defValue)
	#
	mask_resample_linear( &oriArr_c[0,0,0] , &resArr_c[0,0,0] , &sizeArr_c[0] , &originArr_c[0] , &spacingArr_c[0] , &tsizeArr_c[0], &toriginArr_c[0], &tspacingArr_c[0], defValue_c )
	resArr = np.ascontiguousarray(tonumpymask(&resArr_c[0,0,0], &tsizeArr_c[0]))
	return np.copy(resArr)

c和Python调用最坑的地方就是字符串类型和numpy数据类型的存储形式。

3、编译

实际上的编译方式有很多，这里采用setup.py，下面的程序大部分没有，但是留着便于以后好改。

# Enable cython support for slightly faster eval scripts:
# python -m pip install cython numpy
# CYTHONIZE_eval= python setup.py build_ext --inplace

import os
from setuptools import setup, find_packages

include_dirs = []
ext_modules = []

from Cython.Build import cythonize
import numpy as np
include_dirs = [np.get_include()]

os.environ["CC"] = "g++"
os.environ["CXX"] = "g++"

pyxFile = os.path.join('volume_resample',"volume_resample.pyx")
ext_modules = cythonize(pyxFile)



config = {
    'name': '',
    'description': '',
    'long_description': '',
    'long_description_content_type': "text/markdown",
    'author': '',
    'url': '',
    'author_email': '',
    'license': '',
    'version': '0.0.0',
    'install_requires': [],
    'setup_requires': [],
    'extras_require': {
    },
    'packages': find_packages(),
    'scripts': [],
    'entry_points': {},
    'package_data': {},
    'ext_modules': ext_modules,
    'include_dirs': include_dirs
}

setup(**config)

python -m pip install cython numpy这就是编译语句，之前也做过相应的Python打包成pyd大差不差把。

4、调用

pyd调用和以前Python调用是一样的，直接import就好了

import numpy as np


CSUPPORT = True
originArr =np.array([1.0,1.0,1.0],np.float32)
toriginArr =np.array([1.0,1.0,1.0],np.float32)

sizeArr=np.array([4,3,2],np.uint32)
tsizeArr=np.array([8,6,4],np.uint32)

spacingArr =np.array([1.0,1.0,1.0],np.float32)
tspacingArr =np.array([0.5,0.5,0.5],np.float32)

# Check if C-Support is available for better performance
oriArr = np.random.randint(0,10,sizeArr,dtype=np.int16)
# resArr = np.zeros_like(oriArr)
resArr = np.random.randint(0,10,tsizeArr,dtype=np.int16)

defValue=np.int16(0)
if CSUPPORT:
    try:
        from volume_resample import volume_resample
    except:
        CSUPPORT = False
if (CSUPPORT):
        # using cython
        volume_resample.cvolume_resample_linear(oriArr, resArr, sizeArr, originArr,spacingArr,tsizeArr,toriginArr,tspacingArr,defValue)