def _get_fibre_image(network, cpores, vox_len, fibre_rad):
    r"""
    Produce image by filling in voxels along throat edges using Bresenham line
    Then performing distance transform on fibre voxels to erode the pore space
    """

    cthroats = network.find_neighbor_throats(pores=cpores)

    # Below method copied from geometry model throat.vertices
    # Needed now as network may not have all throats assigned to geometry
    # i.e network['throat.vertices'] could return garbage
    verts = _sp.ndarray(network.num_throats(), dtype=object)
    for i in range(len(verts)):
        verts[i] = _sp.asarray(list(network["throat.vert_index"][i].values()))
    cverts = verts[cthroats]
    [vxmin, vxmax, vymin, vymax, vzmin, vzmax] = _get_vertex_range(cverts)
    # Translate vertices so that minimum occurs at the origin
    for index in range(len(cverts)):
        cverts[index] -= np.array([vxmin, vymin, vzmin])
    # Find new size of image array
    cdomain = np.around(np.array([(vxmax-vxmin),
                                  (vymax-vymin),
                                  (vzmax-vzmin)]), 6)
    logger.info("Creating fibre domain range: " + str(np.around(cdomain, 5)))
    lx = np.int(np.around(cdomain[0]/vox_len)+1)
    ly = np.int(np.around(cdomain[1]/vox_len)+1)
    lz = np.int(np.around(cdomain[2]/vox_len)+1)
    # Try to create all the arrays we will need at total domain size
    try:
        pore_space = np.ones([lx, ly, lz], dtype=np.uint8)
        fibre_space = np.zeros(shape=[lx, ly, lz], dtype=np.uint8)
        dt = np.zeros([lx, ly, lz], dtype=float)
        # Only need one chunk
        cx = cy = cz = 1
        chunk_len = np.max(np.shape(pore_space))
    except:
        logger.info("Domain too large to fit into memory so chunking domain"
                    "to process image, this may take some time")
        # Do chunking
        chunk_len = 100
        if (lx > chunk_len):
            cx = np.ceil(lx/chunk_len).astype(int)
        else:
            cx = 1
        if (ly > chunk_len):
            cy = np.ceil(ly/chunk_len).astype(int)
        else:
            cy = 1
        if (lz > chunk_len):
            cz = np.ceil(lz/chunk_len).astype(int)
        else:
            cz = 1

    # Get image of the fibres
    line_points = bresenham(cverts, vox_len/2)
    line_ints = (np.around((line_points/vox_len), 0)).astype(int)
    for x, y, z in line_ints:
        try:
            pore_space[x][y][z] = 0
        except __HOLE__:
            logger.warning("Some elements in image processing are out" +
                           "of bounds")

    num_chunks = np.int(cx*cy*cz)
    cnum = 1
    for ci in range(cx):
        for cj in range(cy):
            for ck in range(cz):
                # Work out chunk range
                logger.info("Processing Fibre Chunk: "+str(cnum)+" of " +
                            str(num_chunks))
                cxmin = ci*chunk_len
                cxmax = np.int(np.ceil((ci+1)*chunk_len + 5*fibre_rad))
                cymin = cj*chunk_len
                cymax = np.int(np.ceil((cj+1)*chunk_len + 5*fibre_rad))
                czmin = ck*chunk_len
                czmax = np.int(np.ceil((ck+1)*chunk_len + 5*fibre_rad))
                # Don't overshoot
                if cxmax > lx:
                    cxmax = lx
                if cymax > ly:
                    cymax = ly
                if czmax > lz:
                    czmax = lz
                dt = ndimage.distance_transform_edt(pore_space[cxmin:cxmax,
                                                               cymin:cymax,
                                                               czmin:czmax])
                fibre_space[cxmin:cxmax,
                            cymin:cymax,
                            czmin:czmax][dt <= fibre_rad] = 0
                fibre_space[cxmin:cxmax,
                            cymin:cymax,
                            czmin:czmax][dt > fibre_rad] = 1
                del dt
                cnum += 1
    del pore_space
    return fibre_space
dataset/ETHPy150Open PMEAL/OpenPNM/OpenPNM/Geometry/models/pore_volume.py/_get_fibre_image