estimate background level through finding the most homogeneous area
    and take its average
    
    min_size - relative size of the examined area
772456
def estimateBackgroundLevel(img, image_is_artefact_free=False, 
                            min_rel_size=0.05, max_abs_size=11):
    '''
    
    '''

    s0,s1 = img.shape[:2]
    s = min(max_abs_size, int(max(s0,s1)*min_rel_size))
    arr = np.zeros(shape=(s0-2*s, s1-2*s), dtype=img.dtype)
    
    #fill arr:
    _spatialStd(img, arr, s)
    #most homogeneous area:
    i,j = np.unravel_index(arr.argmin(), arr.shape)
    sub = img[int(i+0.5*s):int(i+s*1.5), 
              int(j+s*0.5):int(j+s*1.5)]

    return np.median(sub)
294565
def backgroundMMD(self, catalog, method='cloud-in-cells', weights=None):
        

        # Select objects in annulus
        cut_annulus = self.roi.inAnnulus(catalog.lon,catalog.lat)
        mag_1 = catalog.mag_1[cut_annulus]
        mag_2 = catalog.mag_2[cut_annulus]

        # Units are (deg^2)
        solid_angle = ugali.utils.binning.take2D(self.solid_angle_mmd, mag_2, mag_1,
                                                 self.roi.bins_mag, self.roi.bins_mag)

        # Weight each object before binning
        # Divide by solid angle and bin size in magnitudes to get number density 
        # [objs / deg^2 / mag^2]
        if weights is None:
            number_density = (solid_angle*self.roi.delta_mag**2)**(-1)
        else:
            number_density = weights*(solid_angle*self.roi.delta_mag**2)**(-1)

        method = str(method).lower()
        if method == 'cloud-in-cells':
            # Apply cloud-in-cells algorithm
            mmd_background = ugali.utils.binning.cloudInCells(mag_2,mag_1,
                                                              [self.roi.bins_mag,self.roi.bins_mag],
                                                              weights=number_density)[0]
        elif method == 'bootstrap':
            # Not implemented
            raise ValueError("Bootstrap method not implemented")
            mag_1 + (mag_1_err * np.random.normal(0, 1., len(mag_1)))
            mag_2 + (mag_2_err * np.random.normal(0, 1., len(mag_2)))

        elif method == 'histogram':
            # Apply raw histogram
            mmd_background = np.histogram2d(mag_1,mag_2,bins=[self.roi.bins_mag,self.roi.bins_mag],
                                            weights=number_density)[0]

        elif method == 'kde':
            # Gridded kernel density estimator
            logger.warning("### KDE not implemented properly")
            mmd_background = ugali.utils.binning.kernelDensity(mag_2,mag_1,
                                                               [self.roi.bins_mag,self.roi.bins_mag],
                                                               weights=number_density)[0]
        elif method == 'uniform':
            logger.warning("### WARNING: Uniform CMD")
            hist = np.histogram2d(mag_1,mag_2,bins=[self.roi.bins_mag,self.roi.bins_mag], 
                                  weights=number_density)[0]
            mmd_background = np.mean(hist)*np.ones(hist.shape)
            observable = (self.solid_angle_mmd > self.minimum_solid_angle)
            mmd_background *= observable
            return mmd_background
        else:
            raise ValueError("Unrecognized method: %s"%method)
        ## Account for the objects that spill out of the observable space
        ## But what about the objects that spill out to red colors??
        #for index_color in range(0, len(self.roi.centers_color)):
        #    for index_mag in range(0, len(self.roi.centers_mag)):
        #        if self.solid_angle_cmd[index_mag][index_color] < self.minimum_solid_angle:
        #            cmd_background[index_mag - 1][index_color] += cmd_background[index_mag][index_color]
        #            cmd_background[index_mag][index_color] = 0.
        #            break

        mmd_area = self.solid_angle_mmd*self.roi.delta_mag**2 # [deg^2 * mag^2]

        # ADW: This accounts for leakage to faint magnitudes
        # But what about the objects that spill out to red colors??
        # Maximum obsevable magnitude index for each color (uses the fact that
        # np.argmin returns first minimum (zero) instance found.
        # NOTE: More complicated maps may have holes causing problems

        observable = (self.solid_angle_mmd > self.minimum_solid_angle)
        index_mag_1 = observable.argmin(axis=0) - 1
        index_mag_2 = np.arange(len(self.roi.centers_mag))
        # Add the cumulative leakage back into the last bin of the CMD
        leakage = (mmd_background * ~observable).sum(axis=0)
        ### mmd_background[[index_mag_1,index_mag_2]] += leakage
        # Zero out all non-observable bins
        ### mmd_background *= observable

        # Avoid dividing by zero by setting empty bins to the value of the 
        # minimum filled bin of the CMD. This choice is arbitrary and 
        # could be replaced by a static minimum, some fraction of the 
        # CMD maximum, some median clipped minimum, etc. However, should 
        # be robust against outliers with very small values.
        min_mmd_background = max(mmd_background[mmd_background > 0.].min(),
                                 1e-4*mmd_background.max())
        mmd_background[observable] = mmd_background[observable].clip(min_mmd_background)

        ### # ADW: This is a fudge factor introduced to renormalize the CMD
        ### # to the number of input stars in the annulus. While leakage
        ### # will still smooth the distribution, it shouldn't result in 
        ### fudge_factor = len(mag) / float((cmd_background*cmd_area).sum())
        ### cmd_background *= fudge_factor

        return mmd_background
1094613
private void findLocalScaleSpaceMax(int []size, int level, int skip) {
		int index0 = spaceIndex;
		int index1 = (spaceIndex + 1) % 3;
		int index2 = (spaceIndex + 2) % 3;

		ImageBorder_F32 inten0 = (ImageBorder_F32)FactoryImageBorderAlgs.value(intensity[index0], 0);
		GrayF32 inten1 = intensity[index1];
		ImageBorder_F32 inten2 = (ImageBorder_F32)FactoryImageBorderAlgs.value(intensity[index2], 0);

		// find local maximums in image 2D space.  Borders need to be ignored since
		// false positives are found around them as an artifact of pixels outside being
		// treated as being zero.
		foundFeatures.reset();
		extractor.setIgnoreBorder(size[level] / (2 * skip));
		extractor.process(intensity[index1],null,null,null,foundFeatures);

		// Can't consider feature which are right up against the border since they might not be a true local
		// maximum when you consider the features on the other side of the ignore border
		int ignoreRadius = extractor.getIgnoreBorder() + extractor.getSearchRadius();
		int ignoreWidth = intensity[index1].width-ignoreRadius;
		int ignoreHeight = intensity[index1].height-ignoreRadius;

		// number of features which can be added
		int numberRemaining;

		// if configured to do so, only select the features with the highest intensity
		QueueCorner features;
		if( sortBest != null ) {
			sortBest.process(intensity[index1],foundFeatures,true);
			features = sortBest.getBestCorners();
			numberRemaining = maxFeaturesPerScale;
		} else {
			features = foundFeatures;
			numberRemaining = Integer.MAX_VALUE;
		}

		int levelSize = size[level];
		int sizeStep = levelSize-size[level-1];

		// see if these local maximums are also a maximum in scale-space
		for( int i = 0; i < features.size && numberRemaining > 0; i++ ) {
			Point2D_I16 f = features.get(i);

			// avoid false positives.  see above comment
			if( f.x < ignoreRadius || f.x >= ignoreWidth || f.y < ignoreRadius || f.y >= ignoreHeight )
				continue;

			float val = inten1.get(f.x,f.y);

			// see if it is a max in scale-space too
			if( checkMax(inten0,val,f.x,f.y) && checkMax(inten2,val,f.x,f.y) ) {

				// find the feature's location to sub-pixel accuracy using a second order polynomial
				// NOTE: In the original paper this was done using a quadratic.  See comments above.
				// NOTE: Using a 2D polynomial for x and y might produce better results.
				float peakX = polyPeak(inten1.get(f.x-1,f.y),inten1.get(f.x,f.y),inten1.get(f.x+1,f.y));
				float peakY = polyPeak(inten1.get(f.x,f.y-1),inten1.get(f.x,f.y),inten1.get(f.x,f.y+1));
				float peakS = polyPeak(inten0.get(f.x,f.y),inten1.get(f.x,f.y),inten2.get(f.x,f.y));

				float interpX = (f.x+peakX)*skip;
				float interpY = (f.y+peakY)*skip;
				float interpS = levelSize+peakS*sizeStep;

				double scale =  1.2*interpS/9.0;
				foundPoints.grow().set(interpX,interpY,scale);
				numberRemaining--;
			}
		}
	}
263327
def get_sub_area(area, xslice, yslice):
    """"""
    new_area_extent = ((area.pixel_upper_left[0] +
                        (xslice.start - 0.5) * area.pixel_size_x),
                       (area.pixel_upper_left[1] -
                        (yslice.stop - 0.5) * area.pixel_size_y),
                       (area.pixel_upper_left[0] +
                        (xslice.stop - 0.5) * area.pixel_size_x),
                       (area.pixel_upper_left[1] -
                        (yslice.start - 0.5) * area.pixel_size_y))

    return AreaDefinition(area.area_id, area.name,
                          area.proj_id, area.proj_dict,
                          xslice.stop - xslice.start,
                          yslice.stop - yslice.start,
                          new_area_extent)
1225725
protected void detectOctave( II integral , int skip , int ...featureSize ) {

		int w = integral.width/skip;
		int h = integral.height/skip;

		// resize the output intensity image taking in account subsampling
		for( int i = 0; i < intensity.length; i++ ) {
			intensity[i].reshape(w,h);
		}

		// compute feature intensity in each level
		for( int i = 0; i < featureSize.length; i++ ) {
			GIntegralImageFeatureIntensity.hessian(integral,skip,featureSize[i],intensity[spaceIndex]);

			spaceIndex++;
			if( spaceIndex >= 3 )
				spaceIndex = 0;

			// find maximum in scale space
			if( i >= 2 ) {
				findLocalScaleSpaceMax(featureSize,i-1,skip);
			}
		}
	}
533005
def composite(
    background_image,
    foreground_image,
    foreground_width_ratio=0.25,
    foreground_position=(0.0, 0.0),
):
    """"""

    if foreground_width_ratio <= 0:
        return background_image

    composite = background_image.copy()
    width = int(foreground_width_ratio * background_image.shape[1])
    foreground_resized = resize(foreground_image, width)
    size = foreground_resized.shape

    x = int(foreground_position[1] * (background_image.shape[1] - size[1]))
    y = int(foreground_position[0] * (background_image.shape[0] - size[0]))

    # TODO: warn if resulting coordinates are out of bounds?
    composite[y : y + size[0], x : x + size[1]] = foreground_resized

    return composite
514669
private void computeRootStatistics(Area root)
    {
    	for (Box box : root.getBoxes())
    	{
    		int len = letterLength(box.getText());
    		if (len > 0)
    		{
    			int key = colorKey(box.getColor());
    			Integer val = colors.get(key);
    			if (val == null) val = 0;
    			val += len;
    			colors.put(key, val);
    		}
    	}
        for (int i = 0; i < root.getChildCount(); i++)
            computeRootStatistics(root.getChildAt(i));
    }
    
    private int letterLength(String s)
    {
        int len = 0;
        for (int i = 0; i < s.length(); i++)
            if (Character.getType(s.charAt(i)) != Character.SPACE_SEPARATOR)
                len++;
        return len;
    }

    /**
     * Converts a color to a number, 4 bits per color channel so that similar colors have similar numbers.
     */
    private int colorKey(Color color)
    {
        return (color.getRed() / 16) * 256 + (color.getGreen() / 16) * 16 + (color.getBlue() / 16);
    }

}
729559
def threshold(image, block_size=DEFAULT_BLOCKSIZE, mask=None):
    
    if mask is None:
        mask = np.zeros(image.shape[:2], dtype=np.uint8)
        mask[:] = 255

    if len(image.shape) > 2 and image.shape[2] == 4:
        image = cv2.cvtColor(image, cv2.COLOR_BGRA2GRAY)
    res = _calc_block_mean_variance(image, mask, block_size)
    res = image.astype(np.float32) - res.astype(np.float32) + 255
    _, res = cv2.threshold(res, 215, 255, cv2.THRESH_BINARY)
    return res
742201
def estimate_density(DATA_PATH, feature_size):
    """"""
    if not os.path.exists(DATA_PATH):
        raise Exception("Data is not there!")
    density = []
    P = 0.01
    for _ in range(10):
        num_non_zero = 0
        num_sample = 0
        with open(DATA_PATH) as f:
            for line in f:
                if (random.random() < P):
                    num_non_zero += len(line.split(" ")) - 1
                    num_sample += 1
        density.append(num_non_zero * 1.0 / (feature_size * num_sample))
    return sum(density) / len(density)
419406
def step(self):
        """"""
        # Assign attributes to local variables for convenience.
        u = self._u
        gI = self._data
        dgI = self._ddata
        theta = self._theta
        v = self._v
        
        if u is None:
            raise ValueError("the levelset is not set (use set_levelset)")
        
        res = np.copy(u)
        
        # Balloon.
        if v > 0:
            aux = binary_dilation(u, self.structure)
        elif v < 0:
            aux = binary_erosion(u, self.structure)
        if v!= 0:
            res[self._threshold_mask_v] = aux[self._threshold_mask_v]
        
        # Image attachment.
        aux = np.zeros_like(res)
        dres = np.gradient(res)
        for el1, el2 in zip(dgI, dres):
            aux += el1*el2
        res[aux > 0] = 1
        res[aux < 0] = 0
        
        # Smoothing.
        for i in range(self.smoothing):
            res = curvop(res)
        
        self._u = res
706398
private void findHomogeneousLeaves(AreaImpl root)
    {
        if (root.getChildCount() > 1)
            findSuperAreas(root, 1);
        for (int i = 0; i < root.getChildCount(); i++)
            findHomogeneousLeaves((AreaImpl) root.getChildAt(i));
    }
1220974
def _get_area_def_uniform_sampling(self, lon0, channel):
        """"""
        logger.debug('Computing area definition')

        if lon0 is not None:
            # Define proj4 projection parameters
            proj_dict = {'a': EQUATOR_RADIUS,
                         'b': POLE_RADIUS,
                         'lon_0': lon0,
                         'h': ALTITUDE,
                         'proj': 'geos',
                         'units': 'm'}

            # Calculate maximum scanning angles
            xmax, ymax = get_geostationary_angle_extent(
                namedtuple('area', ['proj_dict'])(proj_dict))

            # Derive area extent using small angle approximation (maximum
            # scanning angle is ~8.6 degrees)
            llx, lly, urx, ury = ALTITUDE * np.array([-xmax, -ymax, xmax, ymax])
            area_extent = [llx, lly, urx, ury]

            # Original image is oversampled. Create pyresample area definition
            # with uniform sampling in N-S and E-W direction
            if self._is_vis(channel):
                sampling = SAMPLING_NS_VIS
            else:
                sampling = SAMPLING_NS_IR
            pix_size = ALTITUDE * sampling
            area_def = pyresample.geometry.AreaDefinition(
                'goes_geos_uniform',
                '{} geostationary projection (uniform sampling)'.format(self.platform_name),
                'goes_geos_uniform',
                proj_dict,
                np.rint((urx - llx) / pix_size).astype(int),
                np.rint((ury - lly) / pix_size).astype(int),
                area_extent)

            return area_def
        else:
            return None
712859
public int checkBackground( float[] pixelValue , float[] dataRow , int modelIndex ) {

		// see which gaussian is the best fit based on Mahalanobis distance
		int index = modelIndex;
		float bestDistance = maxDistance*numBands;
		float bestWeight = 0;

		int ng; // number of gaussians in use
		for (ng = 0; ng < maxGaussians; ng++, index += gaussianStride) {
			float variance = dataRow[index + 1];
			if (variance <= 0) {
				break;
			}

			float mahalanobis = 0;
			for (int i = 0; i < numBands; i++) {
				float mean = dataRow[index + 2+i];
				float delta = pixelValue[i] - mean;
				mahalanobis += delta * delta / variance;
			}

			if (mahalanobis < bestDistance) {
				bestDistance = mahalanobis;
				bestWeight = dataRow[index];
			}
		}

		if( ng == 0 ) // There are no models. Return unknown
			return unknownValue;
		return bestWeight >= significantWeight ? 0 : 1;
	}
687027
function(poly1, poly2)
{
    var mini1 = findMSA(poly2, poly1.tPlanes);
    if(mini1 == -1) return NONE;
    var min1 = last_MSA_min;
    
    var mini2 = findMSA(poly1, poly2.tPlanes);
    if(mini2 == -1) return NONE;
    var min2 = last_MSA_min;
    
    // There is overlap, find the penetrating verts
    if(min1 > min2)
        return findVerts(poly1, poly2, poly1.tPlanes[mini1].n, min1);
    else
        return findVerts(poly1, poly2, vneg(poly2.tPlanes[mini2].n), min2);
}
905753
def get_unit_primitive_area(self):
        

        A_ads = self.surface_area
        A_clean = self.clean_entry.surface_area
        n = (A_ads / A_clean)
        return n
431512
def _areadist(ax, v, xr, c, bins=100, by=None, alpha=1, label=None):
    
    y, x = np.histogram(v[~np.isnan(v)], bins)
    x = x[:-1]

    if by is None:
        by = np.zeros((bins,))

    ax.fill_between(x, y, by, facecolor=c, alpha=alpha, label=label)
    return y
663270
def _correctArtefacts(self, image, threshold):
        '''
        
        '''
        image = np.nan_to_num(image)
        medianThreshold(image, threshold, copy=False)
        return image