Python / Cython optimization (3D kernel filter)

I'm tryin to optimize this loop in python:

def apply_kernel(X,Z,kernel):
    d=kernel.shape
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            for k in range(X.shape[2]):

                X[i,j,k]=np.sum(kernel*Z[i:(i+d[0]),j:(j+d[1]),k:(k+d[2])])
    return X

I tried with cython:

import cython
import numpy as np
cimport numpy as np


ctypedef np.float64_t cpl_t
cpl = np.float64


@cython.boundscheck(False) # compiler directive
@cython.wraparound(False) # compiler directive
def appKernel( np.ndarray[ cpl_t, ndim=3] X,np.ndarray[ cpl_t, ndim=3] Z,np.ndarray[ cpl_t, ndim=3] kernel):
    cdef Py_ssize_t i, j, k, ik, jk, kk
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            for k in range(X.shape[2]):

                for ik in range(kernel.shape[0]):
                    for jk in range(kernel.shape[1]):
                        for kk in range(kernel.shape[2]):
                            X[i,j,k]=X[i,j,k]+Z[i+ik,j+jk,k+kk]*kernel[ik,jk,kk]


    return X

There is an improvement (about 10 times faster); but I think that it could be better. Do you have any ideas on how to improve this?

Next step would be paralelize this using mpi4py.