def hemi_merge(self, meth='single', weight=None):
        """

        Parameters
        ----------
        meth : 'single' or 'both'.single, keep roi which appear in a single hemisphere;
        both, only keep roi which appear in both hemisphere
        weight: weight for each roi, n_subj x n_roi np.array

        Returns
        -------

        """
        if self.type is 'scalar':
            self.roi_name = [self.roi_name[i] for i in np.arange(0, len(self.roi_name), 2)]
            odd_f = np.arange(0, len(self.feat_name), 2)
            self.feat_name = [self.feat_name[i] for i in odd_f]

            if weight is None:
                weight = np.ones(self.meas.shape)
                weight[np.isnan(self.meas)] = 0
            else:
                weight = np.repeat(weight, self.meas.shape[1]/weight.shape[1], axis=1)

            if meth is 'single':
                for f in odd_f:
                    meas = self.meas[:, f:f+2]
                    bool_nan = np.isnan(self.meas)
                    index = np.logical_xor(bool_nan[:, 0], bool_nan[:, 1])
                    value = np.where(np.isnan(meas[index, 0]), meas[index, 1], meas[index, 0])
                    meas[index, :] = np.repeat(value[..., np.newaxis], 2, axis=1)
            elif meth is 'both':
                    pass

            odd_meas = self.meas[:, odd_f] * weight[:, odd_f]
            even_meas = self.meas[:, odd_f+1] * weight[:, odd_f+1]
            self.meas = (odd_meas + even_meas)/(weight[:, odd_f] + weight[:, odd_f+1])
        else:
            self.roi_name = [self.roi_name[i] for i in np.arange(0, len(self.roi_name), 2)]
            n_subj, n_feat = self.meas.shape
            meas = np.reshape(self.meas, (n_subj, -1, 3))
            odd_f = np.arange(0, meas.shape[1], 2)
            f_index = []
            for i in np.arange(0, meas.shape[1], 2):
                for j in [0, 1, 2]:
                    f_index.append(i*3+j)
            self.feat_name = [self.feat_name[i] for i in f_index]

            if meth is 'single':
                for f in odd_f:
                    f_meas = meas[:, f:f+2, :]
                    bool_nan = np.isnan(np.prod(f_meas, axis=2))
                    index = np.logical_xor(bool_nan[:, 0], bool_nan[:, 1])
                    value = np.where(np.isnan(f_meas[index, 0, :]), f_meas[index, 1, :], f_meas[index, 0, :])
                    meas[index, f:f+2, :] = np.repeat(value[:, np.newaxis, :], 2, axis=1)
                meas[:, odd_f+1, 0] = -meas[:, odd_f+1, 0]
            elif meth is 'both':
                meas[:, odd_f+1, 0] = -meas[:, odd_f+1, 0]

            self.meas = np.reshape((meas[:, odd_f, :] + meas[:, odd_f+1, :])/2, (n_subj, -1))

def calc_link_dis(base_side, side1, side2):
    # print base_side.shape, side1.shape, side2.shape
    ans = np.zeros_like(base_side, dtype=np.float)
    mask = np.ones_like(base_side, dtype=np.bool)

    #point on the link
    mask_on_line = np.logical_and(base_side == side1+side2, mask)
    mask = np.logical_xor(mask, mask_on_line)
    ans[mask_on_line] = 0

    #the adjaceny points on the link is overlapped
    mask_point = np.logical_and(base_side < 1e-10, mask)
    mask = np.logical_xor(mask, mask_point)
    ans[mask_point] = side1[mask_point]

    side1_sqr = side1 * side1
    side2_sqr = side2 * side2
    base_side_sqr = base_side * base_side

    #obtuse case 1
    mask_obtuse1 = np.logical_and(side1_sqr > base_side_sqr + side2_sqr, mask)
    mask = np.logical_xor(mask, mask_obtuse1)
    ans[mask_obtuse1] = side2[mask_obtuse1]

    #obtuse case 2
    mask_obtuse2 = np.logical_and(side2_sqr > base_side_sqr + side1_sqr, mask)
    mask = np.logical_xor(mask, mask_obtuse2)
    ans[mask_obtuse2] = side1[mask_obtuse2]

    #compute height by Heron's formula
    half_p = (base_side[mask] + side1[mask] + side2[mask]) * 0.5 # half perimeter
    area = np.sqrt(half_p * (half_p - side1[mask]) * (half_p - side2[mask]) * (half_p - base_side[mask]))
    ans[mask] = 2 * area / base_side[mask]
    return ans

def pre_compute_threshes_uci(features, label, threshes):
    '''

    :param features:
    :param label:
    :param threshes:
    :return:
    '''
    threshes_cheatsheet = []
    n, dim = np.shape(features)
    label_plus_one = np.array(label) + 1
    n_ones = np.ones((1, n))[0]
    for i in range(dim):
        cur_f = np.array([x[i] for x in features])
        # sorted(cur_f, key= lambda x: x[0])
        c_cs = []
        if threshes[i][0]:
            # discrete feature
            for t in threshes[i][1]:
                cur_r = cur_f - t
                cur_r = np.logical_xor(cur_r, n_ones)
                cur_r = np.logical_xor(cur_r, label_plus_one)
                # w_err = np.dot(cur_r, d)
                c_cs.append(cur_r)
        else:
            # continuous feature
            for t in threshes[i][1]:
                cur_r = np.sign(cur_f - t) + 1
                cur_r = np.logical_xor(cur_r, label_plus_one)
                # w_err = np.dot(cur_r, d)
                # n_err = np.dot(cur_r, n_ones)
                c_cs.append(cur_r)
        threshes_cheatsheet.append(c_cs)
    return threshes_cheatsheet

    def export():  # type: () -> None
        node = onnx.helper.make_node(
            'Xor',
            inputs=['x', 'y'],
            outputs=['xor'],
        )

        # 2d
        x = (np.random.randn(3, 4) > 0).astype(np.bool)
        y = (np.random.randn(3, 4) > 0).astype(np.bool)
        z = np.logical_xor(x, y)
        expect(node, inputs=[x, y], outputs=[z],
               name='test_xor2d')

        # 3d
        x = (np.random.randn(3, 4, 5) > 0).astype(np.bool)
        y = (np.random.randn(3, 4, 5) > 0).astype(np.bool)
        z = np.logical_xor(x, y)
        expect(node, inputs=[x, y], outputs=[z],
               name='test_xor3d')

        # 4d
        x = (np.random.randn(3, 4, 5, 6) > 0).astype(np.bool)
        y = (np.random.randn(3, 4, 5, 6) > 0).astype(np.bool)
        z = np.logical_xor(x, y)
        expect(node, inputs=[x, y], outputs=[z],
               name='test_xor4d')

    def get_land_mask(self, points):
        '''Get a landmask respecting lakes, and ponds in island in lake
            as water

        :param points: List of lat, lon pairs
        :type points: :class:`numpy.ndarray` of shape Nx2
        :return: Boolean land mask
        :rtype: :class:`numpy.ndarray` of shape N
        '''

        lats = points[:, 0]
        lons = points[:, 1]
        west, east, south, north = (lons.min(), lons.max(),
                                    lats.min(), lats.max())

        relevant_polygons = self.get_polygons_within(west, east, south, north)
        relevant_polygons.sort()

        mask = num.zeros(points.shape[0], dtype=num.bool)
        for p in relevant_polygons:
            if (p.is_land() or p.is_antarctic_grounding_line() or
               p.is_island_in_lake()):
                mask += p.contains_points(points)
            elif p.is_lake() or p.is_pond_in_island_in_lake():
                water = p.contains_points(points)
                num.logical_xor(mask, water, out=mask)
        return mask

def score(filename):
    """
    Score individual image files for the genetic algorithm.
    The idea is to derive predictive factors for the langmuir performance
    (i.e., max power) based on the connectivity, phase fractions, domain sizes,
    etc. The scoring function should be based on multivariate fits from a database
    of existing simulations. To ensure good results, use robust regression techniques
    and cross-validate the best-fit.

    :param filename: image file name
    :type filename: str

    :return score (ideally as an estimated maximum power in W/(m^2))
    :rtype float
    """
    # this works around a weird bug in scipy.misc.imread with 1-bit images
    # open them with PIL as 8-bit greyscale "L" and then convert to ndimage
    pil_img = Image.open(filename)
    image = misc.fromimage(pil_img.convert("L"))

    width, height = image.shape
    if width != 256 or height != 256:
        print "Size Error: ", filename

    #    isize = analyze.interface_size(image)
    ads1, std1 = analyze.average_domain_size(image)

    # we now need to invert the image to get the second domain size
    inverted = (image < image.mean())
    ads2, std2 = analyze.average_domain_size(inverted)

    #overall average domain size
    ads = (ads1 + ads2) / 2.0

    # transfer distances
    # connectivity
    td1, connect1, td2, connect2 = analyze.transfer_distance(image)

    spots = np.logical_xor(image,
        ndimage.binary_erosion(image, structure=np.ones((2,2))))
    erosion = np.count_nonzero(spots)

    spots = np.logical_xor(image,
        ndimage.binary_dilation(image, structure=np.ones((2,2))))
    dilation = np.count_nonzero(spots)

    # fraction of phase one
    nonzero = np.count_nonzero(image)
    fraction = float(nonzero) / float(image.size)
    # scores zero at 0, 1 and maximum at 0.5
    ps = fraction*(1.0-fraction)

    # from simulations with multivariate nonlinear regression
    return (-1.98566e8) + (-1650.14)/ads + (-680.92)*math.pow(ads,0.25) + \
           1.56236e7*math.tanh(14.5*(connect1 + 0.4)) + 1.82945e8*math.tanh(14.5*(connect2 + 0.4)) \
           + 2231.32*connect1*connect2 \
           + (-4.72813)*td1 + (-4.86025)*td2 \
           + 3.79109e7*ps**8 \
           + 0.0540293*dilation + 0.0700451*erosion

    def render_points_as_image(self,points,bounds,resolution):
        """Inputs:
        Polygon data: a list of coordinates of points that
        define the corners of a polygon
        Bounds: The top right hand corner of the square inside which all
        of the points in the polygon data will fit (x,x) (these are technically
        coordinates, but they should be the same for the sake of squares)
        Resolution: The resolution of the output image; a single number,
        all output images are square
        Output: a black and white image (stored as a matrix of Booleans)."""

        output_image = np.zeros((resolution,resolution), dtype=bool)

        step_size = bounds[1] * 2.0 / resolution

        #Tack the first point onto the end, to make looping through
        #adjacent pairs of points easier
        points = np.append(points,[points[0]],axis = 0)

        #Make sure all points are positive
        points = points + bounds[1]

        #Scale the points so rounding them to whole numbers will place
        #them within the output resolution
        points = points / step_size

        #Round the points to prevent future rounding errors
        points = np.floor(points)

        for i in range(len(points)-1):
            #For each pair of points
            p1 = points[i]
            p2 = points[i+1]
    
            #Calculate the slope
            slope = (p2[1]-p1[1])/(p2[0]-p1[0])

            #Then for each step (of 1) in the y-direction from p1 to p2
            for y_step in range(int(np.abs(p2[1]-p1[1]))):
                if slope:
                    if p2[1] > p1[1]:
                        #Find which x value corresponds to the new y value (using the slope)
                        new_y = int(p1[1] + y_step)
                        new_x = int(p1[0] + y_step/slope)
                    else:
                        new_y = int(p1[1] - y_step)
                        new_x = int(p1[0] - y_step/slope)
                    #Then invert every pixel to the left of the new point.
                    #This very nicely fills in the shape, regardless of concavity/convexity.
                    output_image[-new_y][0:new_x] = np.logical_xor(True,output_image[-new_y][0:new_x])

        for point in points[:-1]:
            #The above algorithm consistently leaves a couple corners with lines not inverted correctly
            #This for loop fixes that with only a small increase in runtime
            if output_image[-point[1]][0]:
                output_image[-point[1]][0:point[0]] = np.logical_xor(True,output_image[-point[1]][0:point[0]])

        return output_image

def intersect(a, b, c, d) :
    a = np.array(a)
    b = np.array(b)
    c = np.array(c)
    d = np.array(d)
    def ccw(a, b, c) :
        cross = (a[0, :]-b[0, :])*(c[1, :]-b[1, :])
        cross -= (a[1, :]-b[1, :])*(c[0, :]-b[0, :])
        return cross > 0
    c1 = np.logical_xor(ccw(a, b, c), ccw(a, b, d))
    c2 = np.logical_xor(ccw(c, d, a), ccw(c, d, b))
    return np.logical_and(c1, c2)

def chk_c(arr):

        xr=0;
	for sw in range(0,arr.size):
		xr=np.logical_xor(xr,arr[sw])
	
	u_vect=np.zeros(arr.size)
	for it in range(0,arr.size):
		 if np.logical_xor(xr,arr[it])==True:
			 u_vect[it]=1
	
	return(u_vect)

def shape_symmetry(image, center, major_axis, attrs={}, debug=False):
    # pad to make image center coincide with symmetry center
    lesion_mask, _ = pad_for_rotation(image[..., 3], center)

    rotated = rotate(lesion_mask, 90-major_axis.angle)
    flipped = rotated[:,::-1]
    diff = np.logical_xor(rotated, flipped)

    pixels_diff = diff.sum() / 2.
    major_ratio = pixels_diff / rotated.sum()

    if debug:
        print """\
==== Shape Symmetry ====
--- Major Axis ---
num of pixels   : %d
shape sym ratio : %.3f
""" % (pixels_diff, major_ratio)

        plt.subplot(231)
        plt.imshow(rotated)
        plt.subplot(232)
        plt.imshow(flipped)
        plt.subplot(233)
        plt.imshow(diff)

    rotated = rotate(lesion_mask, 180-major_axis.angle)
    flipped = rotated[:,::-1]
    diff = np.logical_xor(rotated, flipped)

    pixels_diff = diff.sum() / 2.
    minor_ratio = pixels_diff / rotated.sum()

    if debug:
        print """\
--- Minor Axis ---
num of pixels   : %d
shape sym ratio : %.3f
""" % (pixels_diff, minor_ratio)

        plt.subplot(234)
        plt.imshow(rotated)
        plt.subplot(235)
        plt.imshow(flipped)
        plt.subplot(236)
        plt.imshow(diff)
        plt.show()

    attrs.update([
        ('Shape Asymmetry Major Ratio', major_ratio),
        ('Shape Asymmetry Minor Ratio', minor_ratio),
        ('--Shape Asymmetry Score', (major_ratio > 0.13)*1 + (minor_ratio > 0.15)*1),
    ])

def black_tophat(image, selem=None, out=None):
    """Return black top hat of an image.

    The black top hat of an image is defined as its morphological closing minus
    the original image. This operation returns the dark spots of the image that
    are smaller than the structuring element. Note that dark spots in the
    original image are bright spots after the black top hat.

    Parameters
    ----------
    image : ndarray
        Image array.
    selem : ndarray, optional
        The neighborhood expressed as a 2-D array of 1's and 0's.
        If None, use cross-shaped structuring element (connectivity=1).
    out : ndarray, optional
        The array to store the result of the morphology. If None
        is passed, a new array will be allocated.

    Returns
    -------
    out : array, same shape and type as `image`
        The result of the morphological black top hat.

    Examples
    --------
    >>> # Change dark peak to bright peak and subtract background
    >>> import numpy as np
    >>> from skimage.morphology import square
    >>> dark_on_grey = np.array([[7, 6, 6, 6, 7],
    ...                          [6, 5, 4, 5, 6],
    ...                          [6, 4, 0, 4, 6],
    ...                          [6, 5, 4, 5, 6],
    ...                          [7, 6, 6, 6, 7]], dtype=np.uint8)
    >>> black_tophat(dark_on_grey, square(3))
    array([[0, 0, 0, 0, 0],
           [0, 0, 1, 0, 0],
           [0, 1, 5, 1, 0],
           [0, 0, 1, 0, 0],
           [0, 0, 0, 0, 0]], dtype=uint8)

    """
    if out is image:
        original = image.copy()
    else:
        original = image
    out = closing(image, selem, out=out)
    if np.issubdtype(out.dtype, np.bool_):
        np.logical_xor(out, original, out=out)
    else:
        out -= original
    return out

    def _select_coords(self, coords, prm):
        """
        Refinement of QuadricGM._select_coords; we want to choose the correct
        intersection point for a set of rays and our *truncated* quadric cone
        surface.

        Arguments:
        coords - a 2 by 3 by n array whose each column is the global coordinates
            of one intersection point of a ray with the surface.
        prm - a 2 by n array (CHECK THIS) giving the parametric location on the
            ray where the intersection occurs.

        Returns:
        The index of the selected intersection, or None if neither will do.
        """
        select = N.empty(prm.shape[1])
        select.fill(N.nan)
        # Projects the hit coordinates in a local frame on the z axis.
        height = N.sum(N.linalg.inv(self._working_frame)[None,2,:,None]*N.concatenate((coords, N.ones((2,1,coords.shape[-1]))), axis=1), axis=1)
        # Checks if the local_z-projected hit coords are in the actual height of the furstum
        # and if the parameter is positive so that the ray goes ahead.
        inside = (self.z1 <= height) & (height <= self.z2)
        positive = prm > 1e-10
        hitting = inside & positive
        # Choses between the two intersections offered by the surface.
        select[N.logical_and(*hitting)] = 1
        one_hitting = N.logical_xor(*hitting)
        select[one_hitting] = N.nonzero(hitting.T[one_hitting,:])[1]

        return select

 def encodeMask(M):
     """
     Encode binary mask M using run-length encoding.
     :param   M (bool 2D array)  : binary mask to encode
     :return: R (object RLE)     : run-length encoding of binary mask
     """
     [h, w] = M.shape
     M = M.flatten(order='F')
     N = len(M)
     counts_list = []
     pos = 0
     # counts
     counts_list.append(1)
     diffs = np.logical_xor(M[0:N-1], M[1:N])
     for diff in diffs:
         if diff:
             pos +=1
             counts_list.append(1)
         else:
             counts_list[pos] += 1
     # if array starts from 1. start with 0 counts for 0
     if M[0] == 1:
         counts_list = [0] + counts_list
     return {'size':      [h, w],
            'counts':    counts_list ,
            }

    def _shift2boolean(self,
                       q_mesh_shift,
                       is_gamma_center=False,
                       tolerance=1e-5):
        """
        Tolerance is used to judge zero/half gird shift.
        This value is not necessary to be changed usually.
        """
        if q_mesh_shift is None:
            shift = np.zeros(3, dtype='double')
        else:
            shift = np.array(q_mesh_shift, dtype='double')
    
        diffby2 = np.abs(shift * 2 - np.rint(shift * 2))
        if (diffby2 < 0.01).all(): # zero/half shift
            if is_gamma_center:
                is_shift = [0, 0, 0]
            else: # Monkhorst-pack
                diff = np.abs(shift - np.rint(shift))
                is_shift = list(np.logical_xor((diff > 0.1),
                                               (self._mesh % 2 == 0)) * 1)
        else:
            is_shift = None

        return is_shift

    def _support_recovery_norm(self, X_test, relative=False):
        """ accuracy related error pseudo-norm

        Parameters
        ----------
        X_test : positive-definite, symmetric numpy.ndarray of shape (p, p)
            the target precision matrix

        relative: boolean
            whether the error is given as a percentage or as an absolute
            number of counts


        Returns
        -------
        ell0 pseudo-norm between X_test and the estimator

        """
        if relative:
            p = X_test.shape[0]
            c = p * (p - 1)
        else:
            c = 2.
        return np.sum(np.logical_xor(
            np.abs(self.auxiliary_prec_) > machine_eps(0),
            np.abs(X_test) > machine_eps(0))) / c

	def _select_coords(self, coords, prm):
		"""
		Choose between two intersection points on a quadric surface.
		This implementation extends QuadricGM's behaviour by not choosing
		intersections outside the rectangular aperture.
		
		Arguments:
		coords - a 2 by 3 by n array whose each column is the global coordinates
			of one intersection point of a ray with the sphere.
		prm - the corresponding parametric location on the ray where the
			intersection occurs.

		Returns:
		The index of the selected intersection, or None if neither will do.
		"""
		select = QuadricGM._select_coords(self, coords, prm)

		coords = N.concatenate((coords, N.ones((2,1,coords.shape[2]))), axis = 1)
		# assumed no additional parameters to coords, axis = 1
		local = N.sum(N.linalg.inv(self._working_frame)[None,:2,:,None] * \
			coords[:,None,:,:], axis=2)

		abs_x = abs(local[:,0,:])
		abs_y = abs(local[:,1,:])
		outside = abs_x > self._w
		outside |= abs_y > self._h
		inside = (~outside) & (prm > 1e-6)

		select[~N.logical_or(*inside)] = N.nan
		one_hit = N.logical_xor(*inside)
		select[one_hit] = N.nonzero(inside.T[one_hit,:])[1]

		return select

def confuse(soll, ist):
    assert(soll.shape==ist.shape)
    df = pd.DataFrame(zip(soll, ist, np.logical_not(np.logical_xor(soll,ist))), columns=['soll','ist', 'TF'])
    df['PN'] = df['soll'].apply(lambda x: 'Pos' if x==1 else 'Neg')
    df['confuse'] = df[['TF', 'PN']].apply(lambda (x,y): "{}{}".format(x,y), axis=1)

    confusion_matrix = np.zeros([3,3])
    confusion_matrix[0,0] = df.confuse.value_counts().TrueNeg
    confusion_matrix[1,1] = df.confuse.value_counts().TruePos
    try:
        confusion_matrix[0,1] = df.confuse.value_counts().FalsePos
    except Exception:
        print "have no False Positives"
    try:
        confusion_matrix[1,0] = df.confuse.value_counts().FalseNeg
    except Exception:
        print "have no False Negatives"
    confusion_matrix[0,2] = confusion_matrix[0,0] + confusion_matrix[0,1]
    confusion_matrix[1,2] = confusion_matrix[1,0] + confusion_matrix[1,1]
    confusion_matrix[2,0] = confusion_matrix[0,0] + confusion_matrix[1,0]
    confusion_matrix[2,1] = confusion_matrix[0,1] + confusion_matrix[1,1]
    confusion_matrix[2,2] = confusion_matrix[0,0] + confusion_matrix[0,1] + \
                            confusion_matrix[1,0] + confusion_matrix[1,1]

    accuracy = float(confusion_matrix[0,0] + confusion_matrix[1,1])/confusion_matrix[2,2]
    misclass_rate = float(confusion_matrix[1,0] + confusion_matrix[0,1])/confusion_matrix[2,2]
    assert((accuracy+misclass_rate)==1)

    recall = float(confusion_matrix[1,1]) / confusion_matrix[1,2]
    specificity = float(confusion_matrix[0,0]) / confusion_matrix[0,2]

    print "accuracy={}, misclassifaction={}, recall={}, specificity={}".format(accuracy, misclass_rate, recall, specificity)

    return accuracy, misclass_rate, recall, specificity

def contiguous_regions(condition):
    """Find contiguous ``True`` regions of the boolean array ``condition``.

    Return a 2D array where the first column is the start index of the region
    and the second column is the end index, found on [so-contiguous]_.

    Parameters
    ----------
    condition : bool array

    Returns
    -------
    idx : ``[[i0_0, i0_1], [i1_0, i1_1], ...]``
        A list of integer couples, with the start and end of each ``True`` blocks
        in the original array

    Notes
    -----
    .. [so-contiguous] http://stackoverflow.com/questions/4494404/find-large-number-of-consecutive-values-fulfilling-condition-in-a-numpy-array
    """
    # Find the indices of changes in "condition"
    diff = np.logical_xor(condition[1:], condition[:-1])
    idx, = diff.nonzero()
    # We need to start things after the change in "condition". Therefore,
    # we'll shift the index by 1 to the right.
    idx += 1
    if condition[0]:
        # If the start of condition is True prepend a 0
        idx = np.r_[0, idx]
    if condition[-1]:
        # If the end of condition is True, append the length of the array
        idx = np.r_[idx, condition.size]
    # Reshape the result into two columns
    idx.shape = (-1, 2)
    return idx

 def intervalls(self,t,d, condition, frames = False):
     selection = condition(d)
     #print 'selection',selection
     #selection = numpy.append(selection, [False])
     #~ if the condition evaluates to true for the last frame (selection[-1] == True and selection[0] == False) the following roll-xor combination will lead switch_frame[0] == True
     switch_frames = numpy.logical_xor(numpy.roll(selection, 1), selection)
     switch_frames[0] = selection[0] 
     switch_frames[-1] = False # always drop unfinished intervalls
     
     #print 'switchframes',switch_frames
     # detect where the the condition changes from true to false, the roll will directly mark the first and last frame where the condition is true
     start_end = switch_frames.nonzero()[0] # make the returned 0-dimensional tuple an array
     # we condition is true up to the end, we need drop the last transition to condition = true, since we cannot return a closed interval
     if start_end.shape[0] % 2 == 1:
         start_end = numpy.reshape(start_end[:-1],[start_end.size/2,2])                       # reshape the array to contain start-end pairs
     else:
         start_end = numpy.reshape(start_end,[start_end.size/2,2])                       # reshape the array to contain start-end pairs
     
     # always drop intervalls already started at t=0
     if selection[0]:
         start_end = start_end[1:]
         
     if frames:
         return start_end
     else:
         # return the intervalls where the condition is true
         return numpy.array([t[start_end[:,0]],t[start_end[:,1]]]).transpose()

    def _select_coords(self, coords, prm):
        """
        Choose between two intersection points on a quadric surface.
        This implementation extends QuadricGM's behaviour by not choosing
        intersections higher or lower than half the cylinder height.

        Arguments:
        coords - a 2 by 3 by n array whose each column is the global coordinates
            of one intersection point of a ray with the sphere.
        prm - the corresponding parametric location on the ray where the
            intersection occurs.

        Returns:
        The index of the selected intersection, or None if neither will do.
        """
        select = N.empty(prm.shape[1])
        select.fill(N.nan)
        
        height = N.sum(N.linalg.inv(self._working_frame)[None,2,:,None] * \
            N.concatenate((coords, N.ones((2,1,coords.shape[-1]))), axis=1), axis=1)
        inside = (abs(height) <= self._half_h)
        positive = prm > 1e-10
        
        hitting = inside & positive
        select[N.logical_and(*hitting)] = 1
        one_hitting = N.logical_xor(*hitting)
        select[one_hitting] = N.nonzero(hitting.T[one_hitting,:])[1]

        return select