Gaussian blur

 @param [Float] sigma Spread of Gauss bell.
 @param [Float] max_error Error of approximated filter.

 @return [Node] Result of filter operation.
349748
def gauss_blur( sigma, max_error = 1.0 / 0x100 )
      filter_type = DFLOAT.align typecode
      filter = Sequence[ *Array.gauss_blur_filter( sigma, max_error / dimension ) ].
        to_type filter_type
      ( dimension - 1 ).downto( 0 ).inject self do |retval,i|
        retval.convolve filter
      end
    end
1155134
def calc_fwhm_gaussian(self, arr1d, medv=None, gauss_fn=None):
        
        if gauss_fn is None:
            gauss_fn = self.gaussian

        N = len(arr1d)
        X = np.array(list(range(N)))
        Y = arr1d
        # Fitting works more reliably if we do the following
        # a. subtract sky background
        if medv is None:
            medv = get_median(Y)
        Y = Y - medv
        maxv = Y.max()
        # b. clamp to 0..max (of the sky subtracted field)
        Y = Y.clip(0, maxv)

        # Fit a gaussian
        p0 = [0, N - 1, maxv]              # Inital guess
        # Distance to the target function
        errfunc = lambda p, x, y: gauss_fn(x, p) - y  # noqa
        # Least square fit to the gaussian
        with self.lock:
            # NOTE: without this mutex, optimize.leastsq causes a fatal error
            # sometimes--it appears not to be thread safe.
            # The error is:
            # "SystemError: null argument to internal routine"
            # "Fatal Python error: GC object already tracked"
            p1, success = optimize.leastsq(errfunc, p0[:], args=(X, Y))

        if not success:
            raise IQCalcError("FWHM gaussian fitting failed")

        mu, sdev, maxv = p1
        self.logger.debug("mu=%f sdev=%f maxv=%f" % (mu, sdev, maxv))

        # Now that we have the sdev from fitting, we can calculate FWHM
        fwhm = 2.0 * np.sqrt(2.0 * np.log(2.0)) * sdev
        # some routines choke on numpy values and need "pure" Python floats
        # e.g. when marshalling through a remote procedure interface
        fwhm = float(fwhm)
        mu = float(mu)
        sdev = float(sdev)
        maxv = float(maxv)

        res = Bunch.Bunch(fwhm=fwhm, mu=mu, sdev=sdev, maxv=maxv,
                          fit_fn=gauss_fn, fit_args=[mu, sdev, maxv])
        return res
63206
public int updateMixture( float[] pixelValue , float[] dataRow , int modelIndex ) {

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

		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;
				bestIndex = index;
			}
		}

		// Update the model for the best gaussian
		if( bestIndex != -1 ) {
			// If there is a good fit update the model
			float weight = dataRow[bestIndex];
			float variance = dataRow[bestIndex+1];

			weight += learningRate*(1f-weight);
			dataRow[bestIndex]   = 1; // set to one so that it can't possible go negative

			float sumDeltaSq = 0;
			for (int i = 0; i < numBands; i++) {
				float mean = dataRow[bestIndex+2+i];
				float delta = pixelValue[i]-mean;
				dataRow[bestIndex+2+i] = mean + delta*learningRate/weight;
				sumDeltaSq += delta*delta;
			}
			sumDeltaSq /= numBands;
			dataRow[bestIndex+1] = variance + (learningRate /weight)*(sumDeltaSq*1.2F - variance);
			// 1.2f is a fudge factor. Empirical testing shows that the above equation is biased. Can't be bothered
			// to verify the derivation and see if there's a mistake

			// Update Gaussian weights and prune models
			updateWeightAndPrune(dataRow, modelIndex, ng, bestIndex, weight);

			return weight >= significantWeight ? 0 : 1;
		} else if( ng < maxGaussians ) {
			// if there is no good fit then create a new model, if there is room

			bestIndex = modelIndex + ng*gaussianStride;
			dataRow[bestIndex]   = 1; // weight is changed later or it's the only model
			dataRow[bestIndex+1] = initialVariance;
			for (int i = 0; i < numBands; i++) {
				dataRow[bestIndex+2+i] = pixelValue[i];
			}
			// There are no models. Return unknown
			if( ng == 0 )
				return unknownValue;

			updateWeightAndPrune(dataRow, modelIndex, ng+1, bestIndex, learningRate);
			return 1;
		} else {
			// didn't match any models and can't create a new model
			return 1;
		}
	}
943203
def nan_maximum_filter(arr, ksize):
    '''
    
    '''
    out = np.empty_like(arr)
    _calc(arr, out, ksize//2)
    return out
994634
def freq_filt(bma):
    
    """
    Want to fit linear function to artifact line in freq space,
    Then mask everything near that line at distances of ~5-200 pixels, 
    Or whatever the maximum CCD artifact dimension happens to be, 
    This will depend on scaling - consult CCD map for interval
    """

    #Fill ndv with random data
    bf = malib.randomfill(bma)

    import scipy.fftpack
    f = scipy.fftpack.fft2(bf)
    ff = scipy.fftpack.fftshift(f)

    #Ben suggested a Hahn filter here, remove the low frequency, high amplitude information
    #Then do a second fft?

    #np.log(np.abs(ff))

    #perc = malib.calcperc(np.real(ff), perc=(80, 95))
    #malib.iv(numpy.real(ff), clim=perc)

    #See http://scipy-lectures.github.io/advanced/image_processing/

    #Starting at a,b, compute argmax along vertical axis for restricted range
    #Fit line to the x and y argmax values
    #Mask [argmax[y]-1:argmax[y]+1]

    #Create radial mask
    ff_dim = np.array(ff.shape)
    a,b = ff_dim/2
    n = ff_dim.max()
    y,x = np.ogrid[-a:n-a, -b:n-b]
    r1 = 40 
    r2 = 60 
    ff_mask = np.ma.make_mask(ff)
    radial_mask = (r1**2 <= x**2 + y**2) & (x**2 + y**2 < r2**2)
    #Note issues with rounding indices here
    #Hacked in +1 for testing
    ff_mask[:] = radial_mask[a-ff_dim[0]/2:a+ff_dim[0], b-ff_dim[1]/2:b+1+ff_dim[1]/2]

    #Combine radial and line mask

    #Convert mask to 0-1, then feather

    fm = ff * ff_mask

    #Inverse fft
    bf_filt = scipy.fftpack.ifft2(scipy.fftpack.ifftshift(fm))

    #Apply original mask
    bf_filt = np.ma.masked_array(bf_filt, bma.mask)
1227466
function guess_errors() {
		var el = 0,
			b = [0, 0, 0, 0, 0],
			t = [];

		var deg_lambda = 0,
			lambda = [1, 0, 0, 0, 0]; // error+erasure locator poly

		// Berlekamp-Massey algorithm to determine error+erasure locator polynomial

		for (var r = 0; r < 4; r++) {
			var discr = calc_discrepancy(lambda, r + 1); // Compute discrepancy at the r-th step in poly-form

			if (discr !== 0) {
				deg_lambda = 0;

				for (var i = 0; i < 5; i++) {
					t[i] = lambda[i] ^ gmult(discr, b[i]);

					if (t[i]) deg_lambda = i;
				}

				if (2 * el <= r) {
					el = r + 1 - el;

					for (i = 0; i < 5; i++) {
						b[i] = gmult(lambda[i], ginv(discr));
					}
				}

				lambda = t.slice(); // copy
			}

			b.unshift(0); // shift => mul by x
		}

		// Find roots of the locator polynomial.

		var errloc = find_errors(lambda);

		var errors = errloc.length;

		if (errors < 1 || errors > 2) return false;

		if (deg_lambda !== errors) return false; // deg(lambda) unequal to number of roots => uncorrectable error

		// Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo x**(4)). Also find deg(omega).

		var omega = [0, 0, 0, 0, 0];

		for (var i = 0; i < 4; i++) {
			var t = 0;

			for (var j = 0; j < i; j++) {
				t ^= gmult(syndrome[i + 1 - j], lambda[j]);
			}

			omega[i] = t;
		}

		// Compute error values in poly-form.

		for (r = 0; r < errors; r++) {
			var t = 0;
			var pos = errloc[r];
			var root = 31 - pos;

			for (i = 0; i < 4; i++) // evaluate Omega at alpha^(-i)
			{
				t ^= gmult(omega[i], gexp[(root * i) % 31]);
			}

			if (t) // evaluate Lambda' (derivative) at alpha^(-i); all odd powers disappear
			{
				var denom = gmult(lambda[1], 1) ^ gmult(lambda[3], gexp[(root * 2) % 31]);

				if (denom === 0) return false;

				if (pos > 12) pos -= 14;

				codeword[pos] ^= gmult(t, ginv(denom));
			}
		}

		return true;
	}
47445
def _get_vqsr_training(filter_type, vrn_files, gatk_type):
    
    params = []
    for name, train_info, fname in _get_training_data(vrn_files)[filter_type]:
        if gatk_type == "gatk4":
            params.extend(["--resource:%s,%s" % (name, train_info), fname])
            if filter_type == "INDEL":
                params.extend(["--max-gaussians", "4"])
        else:
            params.extend(["-resource:%s,VCF,%s" % (name, train_info), fname])
            if filter_type == "INDEL":
                params.extend(["--maxGaussians", "4"])
    return params
1366201
def _twofilter_smoothing_ON2(self, t, ti, info, phi, lwinfo):
        
        sp, sw = 0., 0.
        upb = lwinfo.max() + self.wgt[t].lw.max()
        if hasattr(self.model, 'upper_bound_trans'):
            upb += self.model.upper_bound_trans(t + 1)
        # Loop over n, to avoid having in memory a NxN matrix
        for n in range(self.N):
            omegan = np.exp(lwinfo + self.wgt[t].lw[n] - upb
                            + self.model.logpt(t + 1, self.X[t][n],
                                               info.hist.X[ti]))
            sp += np.sum(omegan * phi(self.X[t][n], info.hist.X[ti]))
            sw += np.sum(omegan)
        return sp / sw
884310
def _twofilter_smoothing_ON(self, t, ti, info, phi, lwinfo, return_ess,
                               modif_forward, modif_info):
        
        if modif_info is not None:
            lwinfo += modif_info
        Winfo = rs.exp_and_normalise(lwinfo)
        I = rs.multinomial(Winfo)
        if modif_forward is not None:
            lw = self.wgt[t].lw + modif_forward
            W = rs.exp_and_normalise(lw)
        else:
            W = self.wgt[t].W
        J = rs.multinomial(W)
        log_omega = self.model.logpt(t + 1, self.X[t][J], info.hist.X[ti][I])
        if modif_forward is not None:
            log_omega -= modif_forward[J]
        if modif_info is not None:
            log_omega -= modif_info[I]
        Om = rs.exp_and_normalise(log_omega)
        est = np.average(phi(self.X[t][J], info.hist.X[ti][I]), axis=0,
                         weights=Om)
        if return_ess:
            return (est, 1. / np.sum(Om**2))
        else:
            return est
1021625
def gauss_fltr_opencv(dem, size=3, sigma=1):
    
    import cv2
    dem = malib.checkma(dem)
    dem_cv = cv2.GaussianBlur(dem.filled(np.nan), (size, size), sigma)
    out = np.ma.fix_invalid(dem_cv)
    out.set_fill_value(dem.fill_value)
    return out
246951
def run_box_to_gaussian(logdir, verbose=False):
  
  if verbose:
    logger.info('--- Starting run: box_to_gaussian')

  tf.compat.v1.reset_default_graph()
  tf.compat.v1.set_random_seed(0)

  image = get_image(verbose=verbose)
  blur_radius = tf.compat.v1.placeholder(shape=(), dtype=tf.int32)
  with tf.name_scope('filter'):
    blur_side_length = blur_radius * 2 + 1
    pixel_filter = tf.ones((blur_side_length, blur_side_length))
    pixel_filter = (pixel_filter
                    / tf.cast(tf.size(input=pixel_filter), tf.float32))  # normalize

  iterations = 4
  images = [tf.cast(image, tf.float32) / 255.0]
  for _ in xrange(iterations):
    images.append(convolve(images[-1], pixel_filter))
  with tf.name_scope('convert_to_uint8'):
    images = tf.stack(
        [tf.cast(255 * tf.clip_by_value(image_, 0.0, 1.0), tf.uint8)
         for image_ in images])

  summ = image_summary.op(
      'box_to_gaussian', images, max_outputs=iterations,
      display_name='Gaussian blur as a limit process of box blurs',
      description=('Demonstration of forming a Gaussian blur by '
                   'composing box blurs, each of which can be expressed '
                   'as a 2D convolution.\n\n'
                   'A Gaussian blur is formed by convolving a Gaussian '
                   'kernel over an image. But a Gaussian kernel is '
                   'itself the limit of convolving a constant kernel '
                   'with itself many times. Thus, while applying '
                   'a box-filter convolution just once produces '
                   'results that are noticeably different from those '
                   'of a Gaussian blur, repeating the same convolution '
                   'just a few times causes the result to rapidly '
                   'converge to an actual Gaussian blur.\n\n'
                   'Here, the step value controls the blur radius, '
                   'and the image sample controls the number of times '
                   'that the convolution is applied (plus one). '
                   'So, when *sample*=1, the original image is shown; '
                   '*sample*=2 shows a box blur; and a hypothetical '
                   '*sample*=&infin; would show a true Gaussian blur.\n\n'
                   'This is one ingredient in a recipe to compute very '
                   'fast Gaussian blurs. The other pieces require '
                   'special treatment for the box blurs themselves '
                   '(decomposition to dual one-dimensional box blurs, '
                   'each of which is computed with a sliding window); '
                   'we don&rsquo;t perform those optimizations here.\n\n'
                   '[Here are some slides describing the full process.]'
                   '(%s)\n\n'
                   '%s'
                   % ('http://elynxsdk.free.fr/ext-docs/Blur/Fast_box_blur.pdf',
                      IMAGE_CREDIT)))

  with tf.compat.v1.Session() as sess:
    sess.run(image.initializer)
    writer = tf.summary.FileWriter(os.path.join(logdir, 'box_to_gaussian'))
    writer.add_graph(sess.graph)
    for step in xrange(8):
      if verbose:
        logger.info('--- box_to_gaussian: step: %s' % step)
        feed_dict = {blur_radius: step}
      run_options = tf.compat.v1.RunOptions(trace_level=tf.compat.v1.RunOptions.FULL_TRACE)
      run_metadata = config_pb2.RunMetadata()
      s = sess.run(summ, feed_dict=feed_dict,
                   options=run_options, run_metadata=run_metadata)
      writer.add_summary(s, global_step=step)
      writer.add_run_metadata(run_metadata, 'step_%04d' % step)
    writer.close()
1296908
def filter(self, s, method='exact', order=None):
        
        if self._modulation_first:
            return super(Modulation, self).filter(s, method='exact')
        else:
            # The dot product with each modulated kernel is equivalent to the
            # GFT, as for the localization and the IGFT.
            y = np.empty((self.G.n_vertices, self.G.n_vertices))
            for i in range(self.G.n_vertices):
                x = s * self._kernels.localize(i)
                y[i] = np.sqrt(self.G.n_vertices) * self.G.gft(x)
            return y
848038
def normcdf(x, mu, sigma):
    
    t = x-mu
    y = 0.5*erfcc(-t/(sigma*math.sqrt(2.0)))
    if y > 1.0:
        y = 1.0
    return y
404893
static BufferedImage gaussianBlur(BufferedImage src, BufferedImage dst, int radius) {
        int width  = src.getWidth();
        int height = src.getHeight();

        if (dst == null || dst.getWidth() != width || dst.getHeight() != height || src.getType() != dst.getType()) {
            dst = createColorModelCompatibleImage(src);
        }

        float[] kernel = createGaussianKernel(radius);

        if (src.getType() == BufferedImage.TYPE_INT_ARGB) {
            int[] srcPixels = new int[width * height];
            int[] dstPixels = new int[width * height];

            getPixels(src, 0, 0, width, height, srcPixels);

            // horizontal pass
            blur(srcPixels, dstPixels, width, height, kernel, radius);

            // vertical pass
            // noinspection SuspiciousNameCombination
            blur(dstPixels, srcPixels, height, width, kernel, radius);

            // the result is now stored in srcPixels due to the 2nd pass
            setPixels(dst, 0, 0, width, height, srcPixels);
        } else if (src.getType() == BufferedImage.TYPE_BYTE_GRAY) {
            byte[] srcPixels = new byte[width * height];
            byte[] dstPixels = new byte[width * height];

            getPixels(src, 0, 0, width, height, srcPixels);

            // horizontal pass
            blur(srcPixels, dstPixels, width, height, kernel, radius);

            // vertical pass
            // noinspection SuspiciousNameCombination
            blur(dstPixels, srcPixels, height, width, kernel, radius);

            // the result is now stored in srcPixels due to the 2nd pass
            setPixels(dst, 0, 0, width, height, srcPixels);
        } else {
            throw new IllegalArgumentException("EffectUtils.gaussianBlur() src image is not a supported type, type=[" + src.getType()
                                               + "]");
        }

        return dst;
    }
91265
def normal_errors(axes, covariance_matrix, **kwargs):
    
    level = kwargs.pop('level',1)
    traditional_layout = kwargs.pop('traditional_layout',True)
    d = N.diagonal(covariance_matrix)
    ell = ellipse(**kwargs)

    if axes[2,2] < 0:
        axes *= -1

    # Not sure where this factor comes from but it
    # seems to make things work better
    c1 = 2
    axis_lengths = d[:2]
    f = N.linalg.norm(
        ell*axis_lengths,axis=1)

    e0 = -ell.T*d[2]*c1
    e = N.vstack((e0,f))

    _ = dot(e.T,axes).T

    if traditional_layout:
        lon,lat = stereonet_math.cart2sph(_[2],_[0],-_[1])
    else:
        lon,lat = stereonet_math.cart2sph(-_[1],_[0],_[2])
    return list(zip(lon,lat))
785939
function BlurXFilter()
{
    core.AbstractFilter.call(this,
        // vertex shader
        fs.readFileSync(__dirname + '/blurX.vert', 'utf8'),
        // fragment shader
        fs.readFileSync(__dirname + '/blur.frag', 'utf8'),
        // set the uniforms
        {
            strength: { type: '1f', value: 1 }
        }
    );

    /**
     * Sets the number of passes for blur. More passes means higher quaility bluring.
     *
     * @member {number}
     * @memberof BlurXFilter#
     * @default 1
     */
    this.passes = 1;

    this.strength = 4;
}
932696
def fastFilter(arr, ksize=30, every=None, resize=True, fn='median',
                  interpolation=cv2.INTER_LANCZOS4,
                  smoothksize=0,
                  borderMode=cv2.BORDER_REFLECT):
    '''
    
    '''
    if every is None:
        every = max(ksize//3, 1)
    else:
        assert ksize >= 3*every
    s0,s1 = arr.shape[:2]
    
    ss0 = s0//every
    every = s0//ss0
    ss1 = s1//every
    
    out = np.full((ss0+1,ss1+1), np.nan)
    
    c = {'median':_calcMedian,
         'nanmedian':_calcNanMedian,
         'nanmean':_calcNanMean,
         'mean':_calcMean,
         }[fn]
    ss0,ss1 = c(arr, out, ksize, every)
    out = out[:ss0,:ss1]
    
    if smoothksize:
        out = gaussian_filter(out, smoothksize)
        
    
    if not resize:
        return out
    return cv2.resize(out, arr.shape[:2][::-1],
               interpolation=interpolation)