def remove_noise_with_hsv(self, img):
     # Use number of occurrences to find the standard h, s, v
     # Convert to int so we can sort the colors
     # noinspection PyTypeChecker
     img_int = np.dot(np.rint(img * 255), np.power(256, np.arange(3)))
     color_array = sort_by_occurrence(img_int.flatten())
     # standard color is the 2nd most frequent color
     std_color = color_array[1]
     std_b, mod = divmod(std_color, 256 ** 2)
     std_g, std_r = divmod(mod, 256)
     # noinspection PyTypeChecker
     std_h, std_s, std_v = colors.rgb_to_hsv(np.array([std_r, std_g, std_b]) / 255)
     # print(std_h * 360, std_s * 100, std_v * 100)
     height, width, _ = img.shape
     img_hsv = colors.rgb_to_hsv(img)
     h, s, v = img_hsv[:, :, 0], img_hsv[:, :, 1], img_hsv[:, :, 2]
     h_mask = np.abs(h - std_h) > self.h_tolerance
     s_mask = np.abs(s - std_s) > self.s_tolerance
     delta_v = np.abs(v - std_v)
     v_mask = delta_v > self.v_tolerance
     hsv_mask = np.logical_or(np.logical_or(h_mask, s_mask), v_mask)
     new_img = 1 - delta_v
     new_img[hsv_mask] = 0
     # Three types of grayscale colors in new_img:
     # Type A: 1. Outside noise, or inside point.
     # Type B: between 0 and 1. Outside noise, or contour point.
     # Type C: 0. Inside noise, or background.
     return new_img

    def clean_events(self, events):
        events = events.view(np.recarray)
        events.protocol = self._protocol
        events.montage = self._montage
        events.experiment = self._experiment
        events.exp_version = self._exp_version

        stim_events = np.logical_or(events['type'] == 'STIM', events['type'] == 'STIM_OFF')
        stim_events = np.logical_or(stim_events, events['type'] == 'SHAM')
        stim_event_indices = np.where(stim_events)[0]

        poll_events = np.where(events['type'] == 'NP_POLL')[0]
        first_poll_event = poll_events[0]
        last_poll_event = poll_events[-1]

        # Need the last two events (on/off) before the first np poll and the two events after the last np poll
        stim_before = np.array([index for index in stim_event_indices if index < first_poll_event - 2])
        stim_after = np.array([index for index in stim_event_indices if index > last_poll_event + 2])

        good_range = np.array([index for index in range(len(events)) \
                               if index not in stim_before and index not in stim_after])

        cleaned_events = events[good_range]
        # Remove NP_POLL
        cleaned_events = cleaned_events[cleaned_events['type'] != 'NP_POLL']
        cleaned_events.sort(order='mstime')
        return cleaned_events

def getgeodesicpts(m):
    """
 computes the lat/lon values of the points on the surface of the sphere
 corresponding to a twenty-sided (icosahedral) geodesic.

 @param m: the number of points on the edge of a single geodesic triangle.
 There are 10*(m-1)**2+2 total geodesic points, including the poles.

 @return: C{B{lats, lons}} - rank 1 numpy float32 arrays containing
 the latitudes and longitudes of the geodesic points (in degrees). These
 points are nearly evenly distributed on the surface of the sphere.
    """
    x,y,z = _spherepack.ihgeod(m)
# convert cartesian coords to lat/lon.
    rad2dg = 180./math.pi
    r1 = x*x+y*y
    r = numpy.sqrt(r1+z*z)
    r1 = numpy.sqrt(r1)
    xtmp = numpy.where(numpy.logical_or(x,y),x,numpy.ones(x.shape,numpy.float32))
    ztmp = numpy.where(numpy.logical_or(r1,z),z,numpy.ones(z.shape,numpy.float32))
    lons = rad2dg*numpy.arctan2(y,xtmp)+180.
    lats = rad2dg*numpy.arctan2(r1,ztmp)-90.
    lat = numpy.zeros(10*(m-1)**2+2,numpy.float32)
    lon = numpy.zeros(10*(m-1)**2+2,numpy.float32)
# first two points are poles.
    lat[0] = 90; lat[1] = -90.
    lon[0] = 0.; lon[1] = 0.
    lat[2:] = lats[0:2*(m-1),0:m-1,:].flatten()
    lon[2:] = lons[0:2*(m-1),0:m-1,:].flatten()
    return lat,lon

def diamond_110_001(el, a0, n, crack_surface=[1,1,0], crack_front=[0,0,1],
            skin_x=1.0, skin_y=1.0, vac=5.0):
    nx, ny, nz = n
    third_dir = np.cross(crack_surface, crack_front)
    directions = [ third_dir, crack_surface, crack_front ]
    if np.linalg.det(directions) < 0:
        third_dir = -third_dir
    directions = [ third_dir, crack_surface, crack_front ]
    a = Diamond(el, latticeconstant = a0, size = [ nx,ny,nz ], 
                directions = directions)
    sx, sy, sz = a.get_cell().diagonal()
    a.translate([a0/100,a0/100,a0/100])
    a.set_scaled_positions(a.get_scaled_positions())
    a.center()

    lx  = skin_x*sx/nx
    ly  = skin_y*sy/ny
    r   = a.get_positions()
    g   = np.where(
        np.logical_or(
            np.logical_or(
                np.logical_or(
                    r[:, 0] < lx, r[:, 0] > sx-lx),
                r[:, 1] < ly),
            r[:, 1] > sy-ly),
        np.zeros(len(a), dtype=int),
        np.ones(len(a), dtype=int))
    a.set_array('groups', g)

    a.set_cell([sx+2*vac, sy+2*vac, sz])
    a.translate([vac, vac, 0.0])
    a.set_pbc([False, False, True])

    return a

def filter_params(io, rsh, rs, ee, isc):
    # Function filter_params identifies bad parameter sets. A bad set contains
    # Nan, non-positive or imaginary values for parameters; Rs > Rsh; or data
    # where effective irradiance Ee differs by more than 5% from a linear fit
    # to Isc vs. Ee

    badrsh = np.logical_or(rsh < 0., np.isnan(rsh))
    negrs = rs < 0.
    badrs = np.logical_or(rs > rsh, np.isnan(rs))
    imagrs = ~(np.isreal(rs))
    badio = np.logical_or(~(np.isreal(rs)), io <= 0)
    goodr = np.logical_and(~badrsh, ~imagrs)
    goodr = np.logical_and(goodr, ~negrs)
    goodr = np.logical_and(goodr, ~badrs)
    goodr = np.logical_and(goodr, ~badio)

    matrix = np.vstack((ee / 1000., np.zeros(len(ee)))).T
    eff = np.linalg.lstsq(matrix, isc)[0][0]
    pisc = eff * ee / 1000
    pisc_error = np.abs(pisc - isc) / isc
    # check for departure from linear relation between Isc and Ee
    badiph = pisc_error > .05

    u = np.logical_and(goodr, ~badiph)
    return u

    def natural(self):
        """Make a Natural Colors RGB image composite from
        M-bands only.
        """
        self.check_channels('M05', 'M06', 'M07', 'M10')

        ch1 = self['M10'].check_range()
        ch2 = self['M07'].check_range()
        ch3 = self['M05'].check_range()

        ch2b = self['M06'].check_range()
        ch2 = np.ma.where(ch2.mask, ch2b, ch2)

        common_mask = np.logical_or(ch1.mask, ch2.mask)
        common_mask = np.logical_or(common_mask, ch3.mask)
        ch1.mask = common_mask
        ch2.mask = common_mask
        ch3.mask = common_mask

        img = geo_image.GeoImage((ch1, ch2, ch3),
                                 self.area,
                                 self.time_slot,
                                 fill_value=(0, 0, 0),
                                 mode="RGB",
                                 crange=((0, 90),
                                         (0, 90),
                                         (0, 90)))

        img.enhance(gamma=1.8)

        return img

def get_largest_cc(u,v):
    """
    Return mask with largest connected component in u,v

    """

    if not skimage_available:
        print('*** skimage is not available. get_larget_cc() will not work. ***')
        return np.ones_like(u).astype('bool')
    
    fxx = np.array([[1,-2.0,1.0]])
    fxy = np.array([[-0.25,0,0.25],[0.0,0,0],[0.25,0,-0.25]])
    fyy = fxx.T

    u_ = u.astype('float32')
    v_ = v.astype('float32')
    uxx = cv2.filter2D(u_,-1,fxx)
    uxy = cv2.filter2D(u_,-1,fxy)
    uyy = cv2.filter2D(u_,-1,fyy)

    vxx = cv2.filter2D(v_,-1,fxx)
    vxy = cv2.filter2D(v_,-1,fxy)
    vyy = cv2.filter2D(v_,-1,fyy)

    THRESH=0.1
    ue = np.logical_or(np.logical_or(np.abs(uxx)>THRESH, np.abs(uxy)>THRESH),np.abs(uyy)>THRESH)
    ve = np.logical_or(np.logical_or(np.abs(vxx)>THRESH, np.abs(vxy)>THRESH),np.abs(vyy)>THRESH)
    edg = np.logical_or(ue,ve)
    
    L = measure.label(edg.astype('int32'),neighbors=4)
    
    sums = np.bincount(L.ravel())
    biggest_cc = L==np.argmax(sums)
    return biggest_cc

 def remove_wrongly_sized_connected_components(self, a, min_size, max_size, in_place):
     """
     Adapted from http://github.com/jni/ray/blob/develop/ray/morpho.py
     (MIT License)
     """
     bin_out = self.BinaryOut.value
     
     original_dtype = a.dtype
         
     if not in_place:
         a = a.copy()
     if min_size == 0 and (max_size is None or max_size > numpy.prod(a.shape)): # shortcut for efficiency
         return a
     
     try:
         component_sizes = numpy.bincount( a.ravel() )
     except TypeError:
         # On 32-bit systems, must explicitly convert from uint32 to int
         # (This fix is just for VM testing.)
         component_sizes = numpy.bincount( numpy.asarray(a.ravel(), dtype=int) )
     bad_sizes = component_sizes < min_size
     if max_size is not None:
         numpy.logical_or( bad_sizes, component_sizes > max_size, out=bad_sizes )
     
     bad_locations = bad_sizes[a]
     a[bad_locations] = 0
     if (bin_out):
         # Replace non-zero values with 1
         numpy.place(a,a,1)
     return numpy.array(a, dtype=original_dtype)

 def find_most_similar(self, bids, K_sim=14):
     """Return the bid of the most similar book to parameter bid except the given bid."""
     termv = sparse.csc_matrix((self.M, 1), dtype=int)
     for bid in bids:
         col_num = self.bid_to_col.get(str(bid))
         if col_num is not None:
             termv = termv + self.term_bid_matrix.getcol(col_num)
     if termv.nnz == 0:
         return ()
     termva = termv.toarray()    # Generate a vector for terms
     stop_words_removed = np.logical_and(termva, self.stop_words)
     nonzero = stop_words_removed.nonzero()[0]    # Nonzero indices
     rest_term_rows = self.term_bid_matrix_csr[nonzero]
     docs = np.zeros(self.N, dtype=bool)
     for row in rest_term_rows:
         np.logical_or(docs, row.toarray()[0], docs)
     cols = docs.nonzero()[0]
     matched_matrix = self.term_bid_matrix[:,cols]
     termv.data = self.tf(termv.data) * np.array([self.idf(self.row_to_term[row])
                                                  for row in termv.indices])
     termv = normalize(termv.T, axis=1, copy=False)
     matched_matrix.data = self.tf(matched_matrix.data)
     matched_matrix = normalize(matched_matrix.T, axis=1, copy=False).T
     cos_sims = termv.dot(matched_matrix).toarray()[0]
     found_bids = (self.col_to_bid[col] for col in cols)
     return islice((int(r[1])
                    for r in heapq.nlargest(K_sim, zip(cos_sims, found_bids))
                    if int(r[1]) not in bids),
                    9)

 def relative_error(self,A,B,abs_eps):
     absA = np.abs(A)
     absB = np.abs(B)
     I = np.logical_not(np.logical_or(A==B,np.logical_or(absA < abs_eps, absB < abs_eps)))
     E = np.zeros(A.shape,dtype=A.dtype)
     E[I] = np.abs(A[I]-B[I]) / min(absA[I] + absB[I])
     return E

    def prepareCalculation(self, pic=None, automask=True):
        '''
        prepare data used in calculation
        
        :param pic: str, list of str, or 2d array, if provided, and automask is True, then 
            generate a dynamic mask
        :param automask: bool, if True, and pic is not None, then generate a dynamic mask
        
        :return: None
        '''
        self.staticmask = self.mask.staticMask()
        self.correction = self.calculate.genCorrectionMatrix()
        # if a pic is provided, then generate one-time dynamicmask
        if (pic != None) and automask:
            image = self._getPic(pic)
            # image *= self.correction
            dymask = self.mask.dynamicMask(image)
            if dymask != None:
                dymask = np.logical_or(self.staticmask, dymask)
            else:
                dymask = self.staticmask
            
            if self.config.avgmask:
                mask = self.genAvgMask(image, dymask=dymask)
                self.staticmask = np.logical_or(dymask, mask)

        self.calculate.genIntegrationInds(self.staticmask)
        return

    def bb(self,inFile, xval, yval, zval):
        X_invalid = np.logical_or((inFile.header.min[0] > xval), (inFile.header.max[0] < xval))
        Y_invalid = np.logical_or((inFile.header.min[1] > yval), (inFile.header.max[1] < yval))
        Z_invalid = np.logical_or((inFile.header.min[2] > zval), (inFile.header.max[2] < zval))

        bad_indices = np.where(np.logical_or(X_invalid, Y_invalid, Z_invalid))
        return(bad_indices)

def update_mul_inf_zero(lb1_ub1, lb2_ub2, t):
    if not any(np.isinf(lb1_ub1)) and not any(np.isinf(lb2_ub2)):
        return
        
    t_min, t_max = t
    lb1, ub1 = lb1_ub1
    lb2, ub2 = lb2_ub2
    
    ind1_zero_minus = logical_and(lb1<0, ub1>=0)
    ind1_zero_plus = logical_and(lb1<=0, ub1>0)
    
    ind2_zero_minus = logical_and(lb2<0, ub2>=0)
    ind2_zero_plus = logical_and(lb2<=0, ub2>0)
    
    has_plus_inf_1 = logical_or(logical_and(ind1_zero_minus, lb2==-inf), logical_and(ind1_zero_plus, ub2==inf))
    has_plus_inf_2 = logical_or(logical_and(ind2_zero_minus, lb1==-inf), logical_and(ind2_zero_plus, ub1==inf))
    
    # !!!! lines with zero should be before lines with inf !!!!
    ind = logical_or(logical_and(lb1==-inf, ub2==0), logical_and(lb2==-inf, ub1==0))
    t_max[atleast_1d(logical_and(ind, t_max<0.0))] = 0.0
    
    t_max[atleast_1d(logical_or(has_plus_inf_1, has_plus_inf_2))] = inf
    t_max[atleast_1d(logical_or(logical_and(lb1==0, ub2==inf), logical_and(lb2==0, ub1==inf)))] = inf
    
    has_minus_inf_1 = logical_or(logical_and(ind1_zero_plus, lb2==-inf), logical_and(ind1_zero_minus, ub2==inf))
    has_minus_inf_2 = logical_or(logical_and(ind2_zero_plus, lb1==-inf), logical_and(ind2_zero_minus, ub1==inf))
    # !!!! lines with zero should be before lines with -inf !!!!
    t_min[atleast_1d(logical_or(logical_and(lb1==0, ub2==inf), logical_and(lb2==0, ub1==inf)))] = 0.0
    t_min[atleast_1d(logical_or(logical_and(lb1==-inf, ub2==0), logical_and(lb2==-inf, ub1==0)))] = 0.0
    
    t_min[atleast_1d(logical_or(has_minus_inf_1, has_minus_inf_2))] = -inf

def test_categorical_variables():
    np.random.seed(123)

    def objective(x):
        return np.array(np.sum(x, axis=1).reshape(-1, 1))

    carol_spirits = ['past', 'present', 'yet to come']
    encoding = OneHotEncoding(carol_spirits)
    parameter_space = ParameterSpace([
        ContinuousParameter('real_param', 0.0, 1.0),
        CategoricalParameter('categorical_param', encoding)
    ])

    random_design = LatinDesign(parameter_space)
    x_init = random_design.get_samples(10)

    assert x_init.shape == (10, 4)
    assert np.all(np.logical_or(x_init[:, 1:3] == 0.0, x_init[:, 1:3] == 1.0))

    y_init = objective(x_init)

    gpy_model = GPy.models.GPRegression(x_init, y_init)
    gpy_model.Gaussian_noise.fix(1)
    model = GPyModelWrapper(gpy_model)

    loop = ExperimentalDesignLoop(parameter_space, model)
    loop.run_loop(objective, 5)

    assert len(loop.loop_state.Y) == 15
    assert np.all(np.logical_or(loop.loop_state.X[:, 1:3] == 0.0, loop.loop_state.X[:, 1:3] == 1.0))

def computemergepenalty(ellipses,i,j,L,dfore):
    # compute parameters of merged component
    BWmerge = num.logical_or(L == i+1,L == j+1)
    if not BWmerge.any():
        return (0.,ellipses[i])
    ellipsemerge = weightedregionpropsi(BWmerge,dfore[BWmerge])
    #print 'in computemergepenalty, ellipsemerge is: ' + str(ellipsemerge)
    # see if the major, minor, area are small enough
    if (ellipsemerge.area > params.maxshape.area) or (ellipsemerge.minor > params.maxshape.minor) or (ellipsemerge.major > params.maxshape.major):
        #print 'merged ellipse would be too large'
        return (params.maxpenaltymerge+1,ellipses[i])
    # find pixels that should be foreground according to the ellipse parameters
    (r1,r2,c1,c2) = getboundingboxtight(ellipsemerge,L.shape)
    isforepredmerge = ellipsepixels(ellipsemerge,num.array([r1,r2,c1,c2]))
    # pixels that were foreground
    isforepredi = ellipsepixels(ellipses[i],num.array([r1,r2,c1,c2]))
    isforepredj = ellipsepixels(ellipses[j],num.array([r1,r2,c1,c2]))
    isforepredi = num.logical_or(isforepredi, (L[r1:r2,c1:c2]==i+1))
    # pixels that are now foreground that weren't before
    newforemerge = num.logical_and(isforepredmerge,num.logical_or(isforepredi,isforepredj)==False)
    # compute the total background score for this new region that must be foreground
    dforemerge = dfore[r1:r2,c1:c2].copy()
    dforemerge = 1 - dforemerge[newforemerge]
    dforemerge[dforemerge<0] = 0
    mergepenalty = num.sum(dforemerge)
    #print 'mergepenalty = ' + str(mergepenalty)
    #print 'in computemergepenalty, ellipsemerge is: ' + str(ellipsemerge)
    return (mergepenalty,ellipsemerge)

def split_most_even(values):
    """
    >>> import numpy as np
    >>> T = True; F = False
    >>> values = np.array([1, 2, 2, 3, 3, 4])
    >>> expected = np.array([F, F, F, T, T, T])
    >>> actual = split_most_even(values)
    >>> assert(np.all(expected==actual))
    >>> values = np.array([1, 2, 2, 3, 3, 3, 4])
    >>> expected = np.array([F, F, F, T, T, T, T])
    >>> actual = split_most_even(values)
    >>> assert(np.all(expected==actual))
    >>> values = np.array([2, 3, 3, 3, 3, 3, 4])
    >>> expected = np.array([F, T, T, T, T, T, T])
    >>> actual = split_most_even(values)
    >>> assert(np.all(expected==actual))
    >>> values = np.array([2, 3, 3, 3, 3, 3, 4, 4])
    >>> expected = np.array([F, F, F, F, F, T, T])
    >>> actual = split_most_even(values)
    >>> assert(np.all(expected==actual))
    """
    median = np.median(values)
    below = values < median
    above = values > median
    meds = values == median

    num_below = np.sum(below)
    num_above = np.sum(above)
    if num_below < num_above:
        below = np.logical_or(below, meds)
    else:
        above = np.logical_or(above, meds)

    return below, above, median

def generateBoundingPoints(nx, ny, nz):
    x = np.arange(0, nx)
    y = np.arange(0, ny)
    z = np.arange(0, nz)
    xm, ym, zm = np.meshgrid(x, y, z, indexing='ij')

    xm_l = (xm == 0)        # left
    xm_r = (xm == nx - 1)   # right
    xm_bound = np.logical_or(xm_l, xm_r)

    ym_b = (ym == 0)        # back
    ym_f = (ym == ny - 1)   # front
    ym_bound = np.logical_or(ym_b, ym_f)

    zm_b = (zm == 0)        # bottom
    zm_u = (zm == nz - 1)   # up
    zm_bound = np.logical_or(zm_b, zm_u)

    bound = np.logical_or(xm_bound, ym_bound)
    bound = np.logical_or(bound, zm_bound)

    xb = xm[bound]
    yb = ym[bound]
    zb = zm[bound]

    return np.column_stack((xb.ravel(), yb.ravel(), zb.ravel()))

    def _eucl_max(self, nii1, nii2):
        origdata1 = nii1.get_data()
        origdata1 = np.logical_not(
            np.logical_or(origdata1 == 0, np.isnan(origdata1)))
        origdata2 = nii2.get_data()
        origdata2 = np.logical_not(
            np.logical_or(origdata2 == 0, np.isnan(origdata2)))

        if isdefined(self.inputs.mask_volume):
            maskdata = nb.load(self.inputs.mask_volume).get_data()
            maskdata = np.logical_not(
                np.logical_or(maskdata == 0, np.isnan(maskdata)))
            origdata1 = np.logical_and(maskdata, origdata1)
            origdata2 = np.logical_and(maskdata, origdata2)

        if origdata1.max() == 0 or origdata2.max() == 0:
            return np.NaN

        border1 = self._find_border(origdata1)
        border2 = self._find_border(origdata2)

        set1_coordinates = self._get_coordinates(border1, nii1.affine)
        set2_coordinates = self._get_coordinates(border2, nii2.affine)
        distances = cdist(set1_coordinates.T, set2_coordinates.T)
        mins = np.concatenate(
            (np.amin(distances, axis=0), np.amin(distances, axis=1)))

        return np.max(mins)

    def _run_interface(self, runtime):
        nii1 = nb.load(self.inputs.volume1)
        nii2 = nb.load(self.inputs.volume2)

        origdata1 = np.logical_not(np.logical_or(nii1.get_data() == 0, np.isnan(nii1.get_data())))
        origdata2 = np.logical_not(np.logical_or(nii2.get_data() == 0, np.isnan(nii2.get_data())))

        if isdefined(self.inputs.mask_volume):
            maskdata = nb.load(self.inputs.mask_volume).get_data()
            maskdata = np.logical_not(np.logical_or(maskdata == 0, np.isnan(maskdata)))
            origdata1 = np.logical_and(maskdata, origdata1)
            origdata2 = np.logical_and(maskdata, origdata2)

        for method in ("dice", "jaccard"):
            setattr(self, "_" + method, self._bool_vec_dissimilarity(origdata1, origdata2, method=method))

        self._volume = int(origdata1.sum() - origdata2.sum())

        both_data = np.zeros(origdata1.shape)
        both_data[origdata1] = 1
        both_data[origdata2] += 2

        nb.save(nb.Nifti1Image(both_data, nii1.get_affine(), nii1.get_header()), self.inputs.out_file)

        return runtime

def generateBoundingPoints(nx, ny, nz):
    x = np.arange(0, nx)
    y = np.arange(0, ny)
    z = np.arange(0, nz)
    xv, yv, zv = np.meshgrid(x, y, z, indexing='ij')      

    xv_l = (xv == 0)
    xv_r = (xv == nx-1)
    xv_bound = np.logical_or(xv_l, xv_r)
    
    yv_b = (yv == 0)
    yv_f = (yv == ny-1)
    yv_bound = np.logical_or(yv_b, yv_f)
    
    zv_b = (zv == 0)
    zv_u = (zv == nz-1)
    zv_bound = np.logical_or(zv_b, zv_u)
    
    bound = np.logical_or(xv_bound, yv_bound)
    bound = np.logical_or(bound, zv_bound)
    
    xv = xv[bound]
    yv = yv[bound]
    zv = zv[bound]
    
    return np.column_stack((xv.ravel(), yv.ravel(), zv.ravel()))