def gen_mesh_from_voxels(voxels, dims, etype='q', mtype='v'):
    """
    Generate FE mesh from voxels (volumetric data).

    Parameters
    ----------
    voxels : array
        Voxel matrix, 1=material.
    dims : array
        Size of one voxel.
    etype : integer, optional
        'q' - quadrilateral or hexahedral elements
        't' - triangular or tetrahedral elements
    mtype : integer, optional
        'v' - volumetric mesh
        's' - surface mesh

    Returns
    -------
    mesh : Mesh instance
        Finite element mesh.
    """

    dims = dims.squeeze()
    dim = len(dims)
    nddims = nm.array(voxels.shape) + 2

    nodemtx = nm.zeros(nddims, dtype=nm.int8)
    vxidxs = nm.where(voxels)
    set_nodemtx(nodemtx, vxidxs, etype)

    ndidx = nm.where(nodemtx)
    del(nodemtx)

    coors = nm.array(ndidx).transpose() * dims
    nnod = coors.shape[0]

    nodeid = -nm.ones(nddims, dtype=nm.int32)
    nodeid[ndidx] = nm.arange(nnod)

    if mtype == 's':
        felems = []
        nn = nm.zeros(nddims, dtype=nm.int8)

    # generate elements
    if dim == 2:
        ix, iy = vxidxs

        if mtype == 'v':
            elems = nm.array([nodeid[ix,iy],
                              nodeid[ix + 1,iy],
                              nodeid[ix + 1,iy + 1],
                              nodeid[ix,iy + 1]]).transpose()
            edim = 2

        else:
            fc = nm.zeros(nddims + (2,), dtype=nm.int32)
            # x
            fc[ix,iy,:] = nm.array([nodeid[ix,iy + 1],
                                    nodeid[ix,iy]]).transpose()
            fc[ix + 1,iy,:] = nm.array([nodeid[ix + 1,iy],
                                        nodeid[ix + 1,iy + 1]]).transpose()
            nn[ix,iy] = 1
            nn[ix + 1,iy] += 1

            idx = nm.where(nn == 1)
            felems.append(fc[idx])
            # y
            fc.fill(0)
            nn.fill(0)
            fc[ix,iy,:] = nm.array([nodeid[ix,iy],
                                    nodeid[ix + 1,iy]]).transpose()
            fc[ix,iy + 1,:] = nm.array([nodeid[ix + 1,iy + 1],
                                        nodeid[ix,iy + 1]]).transpose()
            nn[ix,iy] = 1
            nn[ix,iy + 1] += 1

            idx = nm.where(nn == 1)
            felems.append(fc[idx])

            elems = nm.concatenate(felems)

            edim = 1

    elif dim == 3:
        ix, iy, iz = vxidxs

        if mtype == 'v':
            elems = nm.array([nodeid[ix,iy,iz],
                              nodeid[ix + 1,iy,iz],
                              nodeid[ix + 1,iy + 1,iz],
                              nodeid[ix,iy + 1,iz],
                              nodeid[ix,iy,iz + 1],
                              nodeid[ix + 1,iy,iz + 1],
                              nodeid[ix + 1,iy + 1,iz + 1],
                              nodeid[ix,iy + 1,iz + 1]]).transpose()
            edim = 3

        else:
            fc = nm.zeros(tuple(nddims) + (4,), dtype=nm.int32)

            # x
            fc[ix,iy,iz,:] = nm.array([nodeid[ix,iy,iz],
                                       nodeid[ix,iy,iz + 1],
                                       nodeid[ix,iy + 1,iz + 1],
                                       nodeid[ix,iy + 1,iz]]).transpose()
            fc[ix + 1,iy,iz,:] = nm.array([nodeid[ix + 1,iy,iz],
                                           nodeid[ix + 1,iy + 1,iz],
                                           nodeid[ix + 1,iy + 1,iz + 1],
                                           nodeid[ix + 1,iy,iz + 1]]).transpose()
            nn[ix,iy,iz] = 1
            nn[ix + 1,iy,iz] += 1

            idx = nm.where(nn == 1)
            felems.append(fc[idx])

            # y
            fc.fill(0)
            nn.fill(0)
            fc[ix,iy,iz,:] = nm.array([nodeid[ix,iy,iz],
                                       nodeid[ix + 1,iy,iz],
                                       nodeid[ix + 1,iy,iz + 1],
                                       nodeid[ix,iy,iz + 1]]).transpose()
            fc[ix,iy + 1,iz,:] = nm.array([nodeid[ix,iy + 1,iz],
                                           nodeid[ix,iy + 1,iz + 1],
                                           nodeid[ix + 1,iy + 1,iz + 1],
                                           nodeid[ix + 1,iy + 1,iz]]).transpose()
            nn[ix,iy,iz] = 1
            nn[ix,iy + 1,iz] += 1

            idx = nm.where(nn == 1)
            felems.append(fc[idx])

            # z
            fc.fill(0)
            nn.fill(0)
            fc[ix,iy,iz,:] = nm.array([nodeid[ix,iy,iz],
                                       nodeid[ix,iy + 1,iz],
                                       nodeid[ix + 1,iy + 1,iz],
                                       nodeid[ix + 1,iy,iz]]).transpose()
            fc[ix,iy,iz + 1,:] = nm.array([nodeid[ix,iy,iz + 1],
                                           nodeid[ix + 1,iy,iz + 1],\
                                           nodeid[ix + 1,iy + 1,iz + 1],
                                           nodeid[ix,iy + 1,iz + 1]]).transpose()
            nn[ix,iy,iz] = 1
            nn[ix,iy,iz + 1] += 1

            idx = nm.where(nn == 1)
            felems.append(fc[idx])

            elems = nm.concatenate(felems)

            edim = 2

    # reduce inner nodes
    if mtype == 's':
        aux = nm.zeros((nnod,), dtype=nm.int32)

        for ii in elems.T:
            aux[ii] = 1

        idx = nm.where(aux)

        aux.fill(0)
        nnod = idx[0].shape[0]

        aux[idx] = range(nnod)
        coors = coors[idx]

        for ii in range(elems.shape[1]):
            elems[:,ii] = aux[elems[:,ii]]

    if etype == 't':
        elems = elems_q2t(elems)

    nel = elems.shape[0]
    nelnd = elems.shape[1]

    mesh = Mesh.from_data('voxel_data',
                          coors, nm.ones((nnod,), dtype=nm.int32),
                          {0: nm.ascontiguousarray(elems)},
                          {0: nm.ones((nel,), dtype=nm.int32)},
                          {0: '%d_%d' % (edim, nelnd)})

    return mesh