<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def insert_cylinder(im, xyz0, xyz1, r):
 r""" Inserts a cylinder of given radius onto a given image Parameters im : array_like Original voxelated image xyz0, xyz1 : 3-by-1 array_like Voxel coordinates of the two end points of the cylinder r : int Radius of the cylinder Returns ------- im : ND-array Original voxelated image overlayed with the cylinder Notes ----- This function is only implemented for 3D images """

    if im.ndim != 3:
        raise Exception('This function is only implemented for 3D images')
    # Converting coordinates to numpy array
    xyz0, xyz1 = [sp.array(xyz).astype(int) for xyz in (xyz0, xyz1)]
    r = int(r)
    L = sp.absolute(xyz0 - xyz1).max() + 1
    xyz_line = [sp.linspace(xyz0[i], xyz1[i], L).astype(int) for i in range(3)]

    xyz_min = sp.amin(xyz_line, axis=1) - r
    xyz_max = sp.amax(xyz_line, axis=1) + r
    shape_template = xyz_max - xyz_min + 1
    template = sp.zeros(shape=shape_template)

    # Shortcut for orthogonal cylinders
    if (xyz0 == xyz1).sum() == 2:
        unique_dim = [xyz0[i] != xyz1[i] for i in range(3)].index(True)
        shape_template[unique_dim] = 1
        template_2D = disk(radius=r).reshape(shape_template)
        template = sp.repeat(template_2D, repeats=L, axis=unique_dim)
        xyz_min[unique_dim] += r
        xyz_max[unique_dim] += -r
    else:
        xyz_line_in_template_coords = [xyz_line[i] - xyz_min[i] for i in range(3)]
        template[tuple(xyz_line_in_template_coords)] = 1
        template = spim.distance_transform_edt(template == 0) <= r

    im[xyz_min[0]:xyz_max[0]+1,
       xyz_min[1]:xyz_max[1]+1,
       xyz_min[2]:xyz_max[2]+1] += template

    return im