def METHOD_NAME(lam, maxK=500, voxelsize=[1,1,1]):
    #from PYME.DSView import View3D
    lam = lam.astype('d') +  1e-2 #to prevent div/0
    fact = (1./lam)
    #print lam.max()
    #print 'number of NaNs = %d' % isnan(pk).sum()
    dx, dy, dz = np.gradient(lam)
    dx *= 1./voxelsize[0]
    dy *= 1./voxelsize[1]
    dz *= 1./voxelsize[2]
    #View3D(dx*dx*fact)
    Exx = (dx*dx*fact).sum(1).sum(0)
    Eyy = (dy*dy*fact).sum(1).sum(0)
    Ezz = (dz*dz*fact).sum(1).sum(0)
    Exy = (dx*dy*fact).sum(1).sum(0)
    Exz = (dx*dz*fact).sum(1).sum(0)
    Eyz = (dy*dz*fact).sum(1).sum(0)
    FIz = np.array([([[Exx[i], Exy[i], Exz[i]],[Exy[i],Eyy[i], Eyz[i]],[Exz[i],Eyz[i], Ezz[i]]]) for i in range(lam.shape[2])])
    return FIz