Implement IRAF zscale algorithm
    nsamples=1000 and contrast=0.25 are the IRAF display task defaults
    image is a 2-d numpy array
    returns (z1, z2)
187662
def zscale(image, nsamples=1000, contrast=0.25):
    

    # Sample the image
    samples = zsc_sample(image, nsamples)

    return zscale_samples(samples, contrast=contrast)
651486
def r2z(r):
    
    with np.errstate(invalid='ignore', divide='ignore'):
        return 0.5 * (np.log(1 + r) - np.log(1 - r))
36769
def image2array(image):
    
    nx, ny = image.shape  # find the size of the array
    imgh = np.reshape(image, nx*ny)  # change the shape to be 1d
    return imgh
1321115
def nside2pixarea(nside, degrees=False):
    """"""
    area = nside_to_pixel_area(nside)
    if degrees:
        return area.to(u.deg ** 2).value
    else:
        return area.to(u.sr).value
175551
def _imshow_array2d(self, array, origin='lower', interpolation='none',
                        aspect='auto', **kwargs):
        
        # NOTE: If you change the defaults for this method, please update
        #       the docstring for `imshow` above.

        # calculate extent
        extent = tuple(array.xspan) + tuple(array.yspan)
        if self.get_xscale() == 'log' and extent[0] == 0.:
            extent = (1e-300,) + extent[1:]
        if self.get_yscale() == 'log' and extent[2] == 0.:
            extent = extent[:2] + (1e-300,) + extent[3:]
        kwargs.setdefault('extent', extent)

        return self.imshow(array.value.T, origin=origin, aspect=aspect,
                           interpolation=interpolation, **kwargs)
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)
1338773
int[] getScaledData(int[] srcImage)
    {
        final int[] dstImage = new int[srcImage.length * SCALE * SCALE];

        for (int y = 0; y < height; y++)
        {
            for (int x = 0; x < width; x++)
            {
                process(srcImage, dstImage, x, y);
            }
        }

        return dstImage;
    }
977833
def array2image(array, nx=0, ny=0):
    
    if nx == 0 or ny == 0:
        n = int(np.sqrt(len(array)))
        if n**2 != len(array):
            raise ValueError("lenght of input array given as %s is not square of integer number!" %(len(array)))
        nx, ny = n, n
    image = array.reshape(int(nx), int(ny))
    return image
64024
def display(self, image):
        
        assert(image.mode == self.mode)
        assert(image.size == self.size)
        self._last_image = image.copy()

        # Send zeros to reset, then pixel values then zeros at end
        sz = image.width * image.height * 4
        buf = bytearray(sz * 3)

        m = self._mapping
        for idx, (r, g, b, a) in enumerate(image.getdata()):
            offset = sz + m[idx] * 4
            brightness = (a >> 4) if a != 0xFF else self._brightness
            buf[offset] = (0xE0 | brightness)
            buf[offset + 1] = b
            buf[offset + 2] = g
            buf[offset + 3] = r

        self._serial_interface.data(list(buf))
486983
Object[] mosim(double S[][], double SF[][], int SR[], int n, double minn[], double maxn[]) {
        int lenS = S[0].length;
        int lenSF = SF[0].length;

//        int e = S.length;
        int r = lenS + lenSF + 1;
        // Define Simplex .. Simplex = [S SF SR];
        double Simplex[][] = new double[S.length][lenS + lenSF + 1];

        for (int i = 0; i < S.length; i++) {
            for (int j = 0; j < lenS; j++) {
                Simplex[i][j] = S[i][j];
            }
            for (int j = 0; j < lenSF; j++) {
                Simplex[i][lenS + j] = SF[i][j];
            }
            Simplex[i][lenS + lenSF] = (double) SR[i];
        }
        // Sort Simplex
        sort(Simplex, r - 1);
        //Assing function values for worst point in Simplex
        double Fw[] = new double[lenSF];
        for (int i = 0; i < lenSF; i++) {
            Fw[i] = Simplex[Simplex.length - 1][n + i];
        }
        // Assing parameter values worst ranked point in Simplex
        double Sw[] = new double[n];
        for (int i = 0; i < n; i++) {
            Sw[i] = Simplex[Simplex.length - 1][i];
        }
        double Sg[] = new double[n];
        // Compute centroid of Simplex after excluding the worst ranked point
        for (int i = 0; i < n; i++) {
            Sg[i] = 0;
            for (int j = 0; j < Simplex.length - 1; j++) {
                Sg[i] += Simplex[j][i];
            }
            Sg[i] /= (double) (Simplex.length - 1);
        }
        // Compute Reflection step
        double Sref[] = new double[n];
        for (int i = 0; i < n; i++) {
            Sref[i] = 2 * Sg[i] - Sw[i];
        }
        boolean accept = isSampleValid(Sref);
        double Snew[], Fnew[];
        if (!accept) {
            // Compute contraction step
            double Scon[] = new double[n];
            for (int i = 0; i < n; i++) {
                Scon[i] = 0.5 * Sg[i] + 0.5 * Sw[i];
            }
            accept = isSampleValid(Scon);
            //if it fails again, try a mutation
            while (!accept) {
                double newpar[][] = lhsu(this.lowBound, this.upBound, 1);
                Scon = newpar[0];
                accept = isSampleValid(Scon);
            }
            double Fcon[] = compob(Scon);
            Snew = Scon;
            Fnew = Fcon;
        // Update number of function evaluations
        } else {
            // Compute corresponding objective function values
            double Fref[] = compob(Sref);
            double SimplexTmp[][] = new double[Simplex.length][lenSF];
            for (int j = 0; j < lenSF; j++) {
                for (int i = 0; i < Simplex.length - 1; i++) {
                    SimplexTmp[i][j] = Simplex[i][n + j];
                }
                SimplexTmp[Simplex.length - 1][j] = Fref[j];
            }
            // Test for non dominance
            Object ret[] = parrank(SimplexTmp, 2); // 2 ist possible wrong
            int Rref[] = (int[]) ret[0];

            int Rrefmax = Integer.MIN_VALUE;
            for (int i = 0; i < Rref.length - 1; i++) {
                if (Rref[i] > Rrefmax) {
                    Rrefmax = Rref[i];
                }
            }
            if (Rref[Rref.length - 1] <= Rrefmax) {
                Snew = Sref;
                Fnew = Fref;
            } else {
                // Compute contraction step
                double Scon[] = new double[n];
                for (int i = 0; i < n; i++) {
                    Scon[i] = 0.5 * Sg[i] + 0.5 * Sw[i];
                }
                accept = isSampleValid(Scon);
                //if it fails again, try a mutation
                while (!accept) {
                    double newpar[][] = lhsu(this.lowBound, this.upBound, 1);
                    Scon = newpar[0];
                    accept = isSampleValid(Scon);
                }
                // Compute corresponding objective function values
                double Fcon[] = compob(Scon);
                Snew = Scon;
                Fnew = Fcon;
            }
        }
        for (int i = 0; i < n; i++) {
            S[S.length - 1][i] = Snew[i];
        }
        for (int i = 0; i < SF[0].length; i++) {
            SF[SF.length - 1][i] = Fnew[i];
        }
        return new Object[]{S, SF};
    }
597703
def pilscale(self, r):
        
        return r * float(self.upsamplefactor) / float(self.binfactor)
1074625
@Override
	public void process( ImagePyramid<GrayF32> image1 ,
						 ImagePyramid<GrayF32> image2 ) {


		// Process the pyramid from low resolution to high resolution
		boolean first = true;
		for( int i = image1.getNumLayers()-1; i >= 0; i-- ) {
			GrayF32 layer1 = image1.getLayer(i);
			GrayF32 layer2 = image2.getLayer(i);

			// declare memory for this layer
			deriv2X.reshape(layer1.width,layer1.height);
			deriv2Y.reshape(layer1.width,layer1.height);
			warpDeriv2X.reshape(layer1.width,layer1.height);
			warpDeriv2Y.reshape(layer1.width,layer1.height);
			warpImage2.reshape(layer1.width,layer1.height);

			// compute the gradient for the second image
			gradient.process(layer2,deriv2X,deriv2Y);

			if( !first ) {
				// interpolate initial flow from previous layer
				interpolateFlowScale(layer1.width, layer1.height);
			} else {
				// for the very first layer there is no information on flow so set everything to 0
				first = false;

				initFlowX.reshape(layer1.width,layer1.height);
				initFlowY.reshape(layer1.width,layer1.height);
				flowX.reshape(layer1.width,layer1.height);
				flowY.reshape(layer1.width,layer1.height);

				ImageMiscOps.fill(flowX,0);
				ImageMiscOps.fill(flowY,0);
				ImageMiscOps.fill(initFlowX,0);
				ImageMiscOps.fill(initFlowY,0);
			}

			// compute flow for this layer
			processLayer(layer1,layer2,deriv2X,deriv2Y);
		}
	}
304163
def positionToIntensityUncertaintyForPxGroup(image, std, y0, y1, x0, x1):
    '''
    
    '''
    fy, fx = y1 - y0, x1 - x0
    if fy != fx:
        raise Exception('averaged area need to be square ATM')
    image = _coarsenImage(image, fx)
    k = _kSizeFromStd(std)
    y0 = int(round(y0 / fy))
    x0 = int(round(x0 / fx))
    arr = image[y0 - k:y0 + k, x0 - k:x0 + k]
    U = positionToIntensityUncertainty(arr, std / fx, std / fx)
    return U[k:-k, k:-k]
267127
def standardDeviation2d(img, ksize=5, blurred=None):
    '''
    
    '''
    if ksize not in (list, tuple):
        ksize = (ksize,ksize)

    if blurred is None:
        blurred = gaussian_filter(img, ksize)
    else:
        assert blurred.shape == img.shape
    
    std = np.empty_like(img)
    
    _calc(img, ksize[0], ksize[1], blurred, std)
    
    return std
474798
function (inData, outData, width, height) {
        var i, n = width * height * 4,
            gr1, gr2, gg1, gg2, gb1, gb2,
            data1 = [],
            data2 = [];

        convolve3x3(inData, data1, width, height,
            [
                [-1, 0, 1],
                [-2, 0, 2],
                [-1, 0, 1]
            ]);

        convolve3x3(inData, data2, width, height,
            [
                [-1, -2, -1],
                [ 0, 0, 0],
                [ 1, 2, 1]
            ]);

        for (i = 0; i < n; i += 4) {
            gr1 = data1[i];
            gr2 = data2[i];
            gg1 = data1[i + 1];
            gg2 = data2[i + 1];
            gb1 = data1[i + 2];
            gb2 = data2[i + 2];

            if (gr1 < 0) {
                gr1 = -gr1;
            }
            if (gr2 < 0) {
                gr2 = -gr2;
            }
            if (gg1 < 0) {
                gg1 = -gg1;
            }
            if (gg2 < 0) {
                gg2 = -gg2;
            }
            if (gb1 < 0) {
                gb1 = -gb1;
            }
            if (gb2 < 0) {
                gb2 = -gb2;
            }

            outData[i] = 255 - (gr1 + gr2) * 0.8;
            outData[i + 1] = 255 - (gg1 + gg2) * 0.8;
            outData[i + 2] = 255 - (gb1 + gb2) * 0.8;
            outData[i + 3] = inData[i + 3];
        }
    }
1260546
def dImage_wrt_2dVerts(observed, visible, visibility, barycentric, image_width, image_height, num_verts, f):
    """"""

    n_channels = np.atleast_3d(observed).shape[2]
    shape = visibility.shape

    # Step 1: get the structure ready, ie the IS and the JS
    IS = np.tile(col(visible), (1, 2*f.shape[1])).ravel()
    JS = col(f[visibility.ravel()[visible]].ravel())
    JS = np.hstack((JS*2, JS*2+1)).ravel()

    pxs = np.asarray(visible % shape[1], np.int32)
    pys = np.asarray(np.floor(np.floor(visible) / shape[1]), np.int32)

    if n_channels > 1:
        IS = np.concatenate([IS*n_channels+i for i in range(n_channels)])
        JS = np.concatenate([JS for i in range(n_channels)])

    # Step 2: get the data ready, ie the actual values of the derivatives
    ksize=1
    sobel_normalizer = cv2.Sobel(np.asarray(np.tile(row(np.arange(10)), (10, 1)), np.float64), cv2.CV_64F, dx=1, dy=0, ksize=ksize)[5,5]
    xdiff = -cv2.Sobel(observed, cv2.CV_64F, dx=1, dy=0, ksize=ksize) / sobel_normalizer
    ydiff = -cv2.Sobel(observed, cv2.CV_64F, dx=0, dy=1, ksize=ksize) / sobel_normalizer

    xdiff = np.atleast_3d(xdiff)
    ydiff = np.atleast_3d(ydiff)

    datas = []

    # The data is weighted according to barycentric coordinates
    bc0 = col(barycentric[pys, pxs, 0])
    bc1 = col(barycentric[pys, pxs, 1])
    bc2 = col(barycentric[pys, pxs, 2])
    for k in range(n_channels):
        dxs = xdiff[pys, pxs, k]
        dys = ydiff[pys, pxs, k]
        if f.shape[1] == 3:
            datas.append(np.hstack((col(dxs)*bc0,col(dys)*bc0,col(dxs)*bc1,col(dys)*bc1,col(dxs)*bc2,col(dys)*bc2)).ravel())
        else:
            datas.append(np.hstack((col(dxs)*bc0,col(dys)*bc0,col(dxs)*bc1,col(dys)*bc1)).ravel())

    data = np.concatenate(datas)

    ij = np.vstack((IS.ravel(), JS.ravel()))
    result = sp.csc_matrix((data, ij), shape=(image_width*image_height*n_channels, num_verts*2))

    return result
622248
def Z_device(self, filter_order=None, window_size=None, tol=0.05):
        '''
        
        '''
        ind = mlab.find(self.fb_resistor >= 0)
        Z1 = np.empty(self.fb_resistor.shape)
        Z1.fill(np.nan)
        # convert to masked array
        Z1 = np.ma.masked_invalid(Z1)

        R2 = self.calibration.R_fb[self.fb_resistor[ind]]
        C2 = self.calibration.C_fb[self.fb_resistor[ind]]
        Z1[ind] = compute_from_transfer_function(self.calibration.hw_version
                                                 .major, 'Z1',
                                                 V1=self.V_total()[ind],
                                                 V2=self.V_fb[ind], R2=R2,
                                                 C2=C2, f=self.frequency)

        Z1 = np.ma.masked_invalid(pd.Series(Z1, pd.to_datetime(self.time, unit='s')
            ).interpolate(method='time').values)
        Z1.fill_value = np.nan
        Z1.data[Z1.mask] = Z1.fill_value

        # if we're filtering and we don't have a window size specified,
        # automatically determine one
        if filter_order and window_size is None:
            window_size = self._get_window_size(tol)

        # if the filter_order or window size is None or if the window size is
        # smaller than filter_order + 2, don't filter
        if (filter_order is None or window_size is None or window_size < filter_order + 2):
            pass
        else:
            # if the window size is less than half the sample length
            if window_size and window_size < len(Z1) / 2:
                # suppress polyfit warnings
                with warnings.catch_warnings():
                    warnings.simplefilter("ignore")
                    Z1 = savgol_filter(Z1, window_size, filter_order)
            else: # fit a line
                result = self.mean_velocity(tol=tol)
                if result['dt'] and \
                    result['dt'] > 0.1 * self.time[-1] and result['p'][0] > 0:
                    if self.calibration._c_drop:
                        c_drop = self.calibration.c_drop(self.frequency)
                    else:
                        c_drop = self.capacitance()[-1] / self.area
                    if self.calibration._c_filler:
                        c_filler = self.calibration.c_filler(self.frequency)
                    else:
                        c_filler = 0
                    x = result['p'][0]*self.time + result['p'][1]
                    C = self.area * (x * (c_drop - c_filler) / \
                                     np.sqrt(self.area) + c_filler)
                    Z1 = 1.0 / (2.0 * math.pi * self.frequency * C)
                    Z1[mlab.find(self.time==result['t_end'])[0]+1:] = \
                        Z1[mlab.find(self.time==result['t_end'])[0]]
                else:
                    Z1 = np.mean(Z1)*np.ones(Z1.shape)
        return Z1