Segmentation of nodes:
        we have the existing nodes
        N.F is the ratio of required nodes and existing nodes
        first, add N nodes to each segment
        then, add one more node to the F first segments

    * assume ordered boundaries

def add_stabilizer_nodes(boundaries_raw, electrodes, nr_nodes_between):
    """
    Segmentation of nodes:
        we have the existing nodes
        N.F is the ratio of required nodes and existing nodes
        first, add N nodes to each segment
        then, add one more node to the F first segments

    * assume ordered boundaries
    """
    boundaries = []

    boundaries = boundaries_raw
    # find first electrode in boundary
    for nr in range(electrodes.shape[0] - 1):
        index0 = np.where(
            (boundaries[:, 0] == electrodes[nr, 0]) &
            (boundaries[:, 1] == electrodes[nr, 1])
        )[0]

        index1 = np.where(
            (boundaries[:, 0] == electrodes[nr + 1, 0]) &
            (boundaries[:, 1] == electrodes[nr + 1, 1])
        )[0]

        index0 = index0[0]
        index1 = index1[0]
        if index1 - index0 < 0:
            index0, index1 = index1, index0
        running_index = index0
        nr_nodes = index1 - index0 - 1
        while nr_nodes < nr_nodes_between:
            # determine line equation
            xy0 = boundaries[running_index, 0:2]
            xy1 = boundaries[running_index + 1, 0:2]

            direction = xy1 - xy0
            heading = direction / np.sqrt(np.sum(direction ** 2))

            # new node
            xy_new = xy0 + heading * direction / 2.0
            a = boundaries[running_index, 2][np.newaxis]
            xyb = np.hstack((xy_new, a))
            boundaries = np.insert(boundaries, running_index + 1, xyb, axis=0)

            # 2, because we have to count the new one
            running_index += 2
            index1 += 1
            nr_nodes += 1

            if running_index == index1:
                running_index = index0

    return boundaries