Compute horizontal offset between input rasters using normalized cross-correlation (NCC) method
509986
def compute_offset_ncc(dem1, dem2, pad=(9,9), prefilter=False, plot=False): 
    

    #Apply edge detection filter up front - improves results when input DEMs are same resolution
    if prefilter:
        print("Applying LoG edge-detection filter to DEMs")
        sigma = 1
        import scipy.ndimage
        #Note, ndimage alone propagates Nans and greatly reduces valid data area
        #Use the malib.nanfill wrapper to avoid this
        dem1 = malib.nanfill(dem1, scipy.ndimage.filters.gaussian_laplace, sigma) 
        dem2 = malib.nanfill(dem2, scipy.ndimage.filters.gaussian_laplace, sigma) 

    import scipy.signal
    #Compute max offset given dem spatial resolution
    #Should implement arbirary x and y search space
    #xsearch = (20, 41)
    #ysearch = (-10, 1)
    stride = 1
    ref = dem1[::stride,::stride]
    kernel = dem2[pad[0]:-pad[1]:stride, pad[0]:-pad[1]:stride]
    #kernel = dem2[-ysearch[0]:-ysearch[1]:stride, xsearch[0]:-xsearch[1]:stride]

    #Normalize
    ref = (ref - ref.mean()) / ref.std()
    kernel = (kernel - kernel.mean()) / kernel.std()

    #Consider using astropy.convolve here instead of scipy.correlate?

    print("Adding random noise to masked regions")
    #Generate random noise to fill gaps before correlation in frequency domain
    #Normal distribution N(mean, std^2)
    #ref_noise = ref.mask * ref.std() * np.random.rand(*ref.shape) + ref.mean()
    #kernel_noise = kernel.mask * kernel.std() * np.random.rand(*kernel.shape) + kernel.mean()
    #This provides noise in proper range, but noise propagates to m, peak is in different locations!
    #ref_noise = ref.mask * (ref.min() + ref.ptp() * np.random.rand(*ref.shape))
    #kernel_noise = kernel.mask * (kernel.min() + kernel.ptp() * np.random.rand(*kernel.shape))

    #This provides a proper normal distribution with mean=0 and std=1
    ref_noise = ref.mask * (np.random.randn(*ref.shape))
    kernel_noise = kernel.mask * (np.random.randn(*kernel.shape))
    #Add the noise
    ref = ref.filled(0) + ref_noise
    kernel = kernel.filled(0) + kernel_noise

    print("Running 2D correlation with search window (x,y): %i, %i" % (pad[1], pad[0]))
    m = scipy.signal.correlate2d(ref, kernel, 'valid')
    #This has memory issues, but ndimage filters can handle nan
    #m = scipy.ndimage.filters.correlate(ref, kernel)
   
    print("Computing sub-pixel peak")
    int_argmax = np.array(np.unravel_index(m.argmax(), m.shape))
    int_offset = int_argmax*stride - pad
    #int_offset = int_argmax*stride + np.array([ysearch[0], xsearch[0]]) 

    print(m.argmax())
    print(m.shape)
    print(int_argmax)
    print(int_offset)

    #Find sub-pixel peak
    sp_argmax = np.array(find_subpixel_peak_position(m, 'parabolic'))
    #May need to split this into integer and decimal components, multipy stride*int and add decimal
    #sp_offset = int_offset + (sp_argmax - int_argmax)
    sp_offset = sp_argmax - pad
    #sp_offset = sp_argmax + np.array([ysearch[0], xsearch[0]]) 

    print(sp_argmax)
    print(sp_offset)

    if plot: 
        fig, ax = plt.subplots()
        ax.set_title('NCC offset, parabolic SPR')
        ax.imshow(m)
        #plt.scatter(*int_argmax[::-1])
        ax.scatter(*sp_argmax[::-1])
    else:
        fig = None

    return m, int_offset, sp_offset, fig
1355980
def dcmsort(folder, copy_series=False, verbose=False):
    ''' 
    '''

    # list files in the input folder
    files = os.listdir(folder)

    srs = {}
    for f in files:
        if os.path.isfile(os.path.join(folder, f)) and f.endswith(dcmext):
            dhdr = dcm.read_file(os.path.join(folder, f))
            #---------------------------------
            # image size
            imsz = np.zeros(2, dtype=np.int64)
            if [0x028,0x010] in dhdr:
                imsz[0] = dhdr[0x028,0x010].value
            if [0x028,0x011] in dhdr:
                imsz[1] = dhdr[0x028,0x011].value
            # voxel size
            vxsz = np.zeros(3, dtype=np.float64)
            if [0x028,0x030] in dhdr and [0x018,0x050] in dhdr:
                pxsz =  np.array([float(e) for e in dhdr[0x028,0x030].value])
                vxsz[:2] = pxsz
                vxsz[2] = float(dhdr[0x018,0x050].value)
            # orientation
            ornt = np.zeros(6, dtype=np.float64)
            if [0x020,0x037] in dhdr:
                ornt = np.array([float(e) for e in dhdr[0x20,0x37].value])
            # seires description, time and study time
            srs_dcrp = dhdr[0x0008, 0x103e].value
            srs_time = dhdr[0x0008, 0x0031].value[:6]
            std_time = dhdr[0x0008, 0x0030].value[:6]

            if verbose:
                print 'series desciption:', srs_dcrp
                print 'series time:', srs_time
                print 'study  time:', std_time
                print '---------------------------------------------------'

            #----------
            # series for any category (can be multiple scans within the same category)
            recognised_series = False
            srs_k = srs.keys()
            for s in srs_k:
                if  np.array_equal(srs[s]['imorient'],  ornt) and \
                    np.array_equal(srs[s]['imsize'],    imsz) and \
                    np.array_equal(srs[s]['voxsize'],   vxsz) and \
                    srs[s]['tseries'] == srs_time:
                    recognised_series = True
                    break
            # if series was not found, create one
            if not recognised_series:
                s = srs_dcrp + '_' + srs_time
                srs[s] = {}
                srs[s]['imorient']  = ornt
                srs[s]['imsize']    = imsz
                srs[s]['voxsize']   = vxsz
                srs[s]['tseries']   = srs_time
            # append the file name
            if 'files' not in srs[s]: srs[s]['files'] = []
            if copy_series:
                srsdir = os.path.join(folder, s)
                create_dir( srsdir )
                shutil.copy(os.path.join(folder, f), srsdir)
                srs[s]['files'].append( os.path.join(srsdir, f) )
            else:
                srs[s]['files'].append( os.path.join(folder, f) )

    return srs
164269
public void computeNccDescriptor( NccFeature f , float x0 , float y0 , float x1 , float y1 ) {
		double mean = 0;
		float widthStep = (x1-x0)/15.0f;
		float heightStep = (y1-y0)/15.0f;

		// compute the mean value
		int index = 0;
		for( int y = 0; y < 15; y++ ) {
			float sampleY = y0 + y*heightStep;
			for( int x = 0; x < 15; x++ ) {
				mean += f.value[index++] = interpolate.get_fast(x0 + x * widthStep, sampleY);
			}
		}
		mean /= 15*15;

		// compute the variance and save the difference from the mean
		double variance = 0;
		index = 0;
		for( int y = 0; y < 15; y++ ) {
			for( int x = 0; x < 15; x++ ) {
				double v = f.value[index++] -= mean;
				variance += v*v;
			}
		}
		variance /= 15*15;
		f.mean = mean;
		f.sigma = Math.sqrt(variance);
	}
1258580
function crop(data, rows, cols, nRows, nCols) {

    if (rows === nRows && cols === nCols) {
        //Do nothing. Returns the same input!!! Be careful
        return data;
    }

    var output = new Array(nCols * nRows);

    var shiftR = Math.floor((rows - nRows) / 2);
    var shiftC = Math.floor((cols - nCols) / 2);
    var destinyRow, sourceRow, i, j;
    for (i = 0; i < nRows; i++) {
        destinyRow = i * nCols;
        sourceRow = (i + shiftR) * cols;
        for (j = 0; j < nCols; j++) {
            output[destinyRow + j] = data[sourceRow + (j + shiftC)];
        }
    }

    return output;
}
529293
def cross_validate(self, ax):
        '''
        

        '''

        # The CDPP to beat
        cdpp_opt = self.get_cdpp_arr()

        # Loop over all chunks
        for b, brkpt in enumerate(self.breakpoints):

            log.info("Cross-validating chunk %d/%d..." %
                     (b + 1, len(self.breakpoints)))

            # Mask for current chunk
            m = self.get_masked_chunk(b)

            # Mask transits and outliers
            time = self.time[m]
            flux = self.fraw[m]
            ferr = self.fraw_err[m]
            med = np.nanmedian(self.fraw)

            # Setup the GP
            gp = GP(self.kernel, self.kernel_params, white=False)
            gp.compute(time, ferr)

            # The masks
            masks = list(Chunks(np.arange(0, len(time)),
                                len(time) // self.cdivs))

            # The pre-computed matrices
            pre_v = [self.cv_precompute(mask, b) for mask in masks]

            # Initialize with the nPLD solution
            log_lam_opt = np.log10(self.lam[b])
            scatter_opt = self.validation_scatter(
                log_lam_opt, b, masks, pre_v, gp, flux, time, med)
            log.info("Iter 0/%d: " % (self.piter) +
                     "logL = (%s), s = %.3f" %
                     (", ".join(["%.3f" % l for l in log_lam_opt]),
                      scatter_opt))

            # Do `piter` iterations
            for p in range(self.piter):

                # Perturb the initial condition a bit
                log_lam = np.array(
                    np.log10(self.lam[b])) * \
                    (1 + self.ppert * np.random.randn(len(self.lam[b])))
                scatter = self.validation_scatter(
                    log_lam, b, masks, pre_v, gp, flux, time, med)
                log.info("Initializing at: " +
                         "logL = (%s), s = %.3f" %
                         (", ".join(["%.3f" % l for l in log_lam]), scatter))

                # Call the minimizer
                log_lam, scatter, _, _, _, _ = \
                    fmin_powell(self.validation_scatter, log_lam,
                                args=(b, masks, pre_v, gp, flux, time, med),
                                maxfun=self.pmaxf, disp=False,
                                full_output=True)

                # Did it improve the CDPP?
                tmp = np.array(self.lam[b])
                self.lam[b] = 10 ** log_lam
                self.compute()
                cdpp = self.get_cdpp_arr()[b]
                self.lam[b] = tmp
                if cdpp < cdpp_opt[b]:
                    cdpp_opt[b] = cdpp
                    log_lam_opt = log_lam

                # Log it
                log.info("Iter %d/%d: " % (p + 1, self.piter) +
                         "logL = (%s), s = %.3f" %
                         (", ".join(["%.3f" % l for l in log_lam]), scatter))

            # The best solution
            log.info("Found minimum: logL = (%s), s = %.3f" %
                     (", ".join(["%.3f" % l for l in log_lam_opt]),
                      scatter_opt))
            self.lam[b] = 10 ** log_lam_opt

        # We're just going to plot lambda as a function of chunk number
        bs = np.arange(len(self.breakpoints))
        color = ['k', 'b', 'r', 'g', 'y']
        for n in range(self.pld_order):
            ax[0].plot(bs + 1, [np.log10(self.lam[b][n])
                                for b in bs], '.', color=color[n])
            ax[0].plot(bs + 1, [np.log10(self.lam[b][n])
                                for b in bs], '-', color=color[n], alpha=0.25)
        ax[0].set_ylabel(r'$\log\Lambda$', fontsize=5)
        ax[0].margins(0.1, 0.1)
        ax[0].set_xticks(np.arange(1, len(self.breakpoints) + 1))
        ax[0].set_xticklabels([])

        # Now plot the CDPP
        cdpp_arr = self.get_cdpp_arr()
        ax[1].plot(bs + 1, cdpp_arr, 'b.')
        ax[1].plot(bs + 1, cdpp_arr, 'b-', alpha=0.25)
        ax[1].margins(0.1, 0.1)
        ax[1].set_ylabel(r'Scatter (ppm)', fontsize=5)
        ax[1].set_xlabel(r'Chunk', fontsize=5)
        ax[1].set_xticks(np.arange(1, len(self.breakpoints) + 1))
667488
protected void findAdjacentRegions(GrayS32 pixelToRegion) {
		// -------- Do the inner pixels first
		if( connect.length == 4 )
			adjacentInner4(pixelToRegion);
		else if( connect.length == 8 ) {
			adjacentInner8(pixelToRegion);
		}

		adjacentBorder(pixelToRegion);
	}
873369
def compute_offset_nuth(dh, slope, aspect, min_count=100, remove_outliers=True, plot=True):
    
    import scipy.optimize as optimization

    if dh.count() < min_count:
        sys.exit("Not enough dh samples")
    if slope.count() < min_count:
        sys.exit("Not enough slope/aspect samples")

    #mean_dh = dh.mean()
    #mean_slope = slope.mean()
    #c_seed = (mean_dh/np.tan(np.deg2rad(mean_slope))) 
    med_dh = malib.fast_median(dh)
    med_slope = malib.fast_median(slope)
    c_seed = (med_dh/np.tan(np.deg2rad(med_slope))) 

    x0 = np.array([0.0, 0.0, c_seed])
  
    print("Computing common mask")
    common_mask = ~(malib.common_mask([dh, aspect, slope]))

    #Prepare x and y data
    xdata = aspect[common_mask].data
    ydata = (dh[common_mask]/np.tan(np.deg2rad(slope[common_mask]))).data

    print("Initial sample count:")
    print(ydata.size)

    if remove_outliers:
        print("Removing outliers")
        #print("Absolute dz filter: %0.2f" % max_dz)
        #diff = np.ma.masked_greater(diff, max_dz)
        #print(diff.count())

        #Outlier dz filter
        f = 3
        sigma, u = (ydata.std(), ydata.mean())
        #sigma, u = malib.mad(ydata, return_med=True)
        rmin = u - f*sigma
        rmax = u + f*sigma
        print("3-sigma filter: %0.2f - %0.2f" % (rmin, rmax))
        idx = (ydata >= rmin) & (ydata <= rmax)
        xdata = xdata[idx]
        ydata = ydata[idx]
        print(ydata.size)

    #Generate synthetic data to test curve_fit
    #xdata = np.arange(0,360,0.01)
    #ydata = f(xdata, 20.0, 130.0, -3.0) + 20*np.random.normal(size=len(xdata))
    
    #Limit sample size
    #n = 10000
    #idx = random.sample(range(xdata.size), n)
    #xdata = xdata[idx]
    #ydata = ydata[idx]

    #Compute robust statistics for 1-degree bins
    nbins = 360
    bin_range = (0., 360.)
    bin_width = 1.0
    bin_count, bin_edges, bin_centers = malib.bin_stats(xdata, ydata, stat='count', nbins=nbins, bin_range=bin_range)
    bin_med, bin_edges, bin_centers = malib.bin_stats(xdata, ydata, stat='median', nbins=nbins, bin_range=bin_range)
    #Needed to estimate sigma for weighted lsq
    #bin_mad, bin_edges, bin_centers = malib.bin_stats(xdata, ydata, stat=malib.mad, nbins=nbins, bin_range=bin_range)
    #Started implementing this for more generic binning, needs testing
    #bin_count, x_bin_edges, y_bin_edges = malib.get_2dhist(xdata, ydata, \
    #        xlim=bin_range, nbins=(nbins, nbins), stat='count')

    """
    #Mask bins in grid directions, can potentially contain biased stats
    #Especially true for SGM algorithm
    #badbins = [0, 90, 180, 270, 360]
    badbins = [0, 45, 90, 135, 180, 225, 270, 315, 360]
    bin_stat = np.ma.masked_where(np.around(bin_edges[:-1]) % 45 == 0, bin_stat)
    bin_edges = np.ma.masked_where(np.around(bin_edges[:-1]) % 45 == 0, bin_edges)
    """

    #Remove any bins with only a few points
    min_bin_sample_count = 9
    idx = (bin_count.filled(0) >= min_bin_sample_count) 
    bin_count = bin_count[idx].data
    bin_med = bin_med[idx].data
    #bin_mad = bin_mad[idx].data
    bin_centers = bin_centers[idx]

    fit = None
    fit_fig = None

    #Want a good distribution of bins, at least 1/4 to 1/2 of sinusoid, to ensure good fit
    #Need at least 3 valid bins to fit 3 parameters in nuth_func
    #min_bin_count = 3
    min_bin_count = 90 
    
    #Not going to help if we have a step function between two plateaus, but better than nothing
    #Calculate bin aspect spread
    bin_ptp = np.cos(np.radians(bin_centers)).ptp()
    min_bin_ptp = 1.0 

    #Should iterate here, if not enough bins, increase bin width
    if len(bin_med) >= min_bin_count and bin_ptp >= min_bin_ptp:

        print("Computing fit")
        #Unweighted fit
        fit = optimization.curve_fit(nuth_func, bin_centers, bin_med, x0)[0]

        #Weight by observed spread in each bin 
        #sigma = bin_mad
        #fit = optimization.curve_fit(nuth_func, bin_centers, bin_med, x0, sigma, absolute_sigma=True)[0]

        #Weight by bin count
        #sigma = bin_count.max()/bin_count
        #fit = optimization.curve_fit(nuth_func, bin_centers, bin_med, x0, sigma, absolute_sigma=False)[0]

        print(fit)

        if plot:
            print("Generating Nuth and Kaab plot")
            bin_idx = np.digitize(xdata, bin_edges)
            output = []
            for i in np.arange(1, len(bin_edges)):
                output.append(ydata[bin_idx==i])
            #flierprops={'marker':'.'}
            lw = 0.25
            whiskerprops={'linewidth':lw}
            capprops={'linewidth':lw}
            boxprops={'facecolor':'k', 'linewidth':0}
            medianprops={'marker':'o', 'ms':1, 'color':'r'}
            fit_fig, ax = plt.subplots(figsize=(6,6))
            #widths = (bin_width/2.0)
            widths = 2.5*(bin_count/bin_count.max())
            #widths = bin_count/np.percentile(bin_count, 50)
            #Stride
            s=3
            #This is inefficient, but we have list of arrays with different length, need to filter
            #Reduntant with earlier filter, should refactor
            bp = ax.boxplot(np.array(output)[idx][::s], positions=bin_centers[::s], widths=widths[::s], showfliers=False, \
                    patch_artist=True, boxprops=boxprops, whiskerprops=whiskerprops, capprops=capprops, \
                    medianprops=medianprops)
            bin_ticks = [0, 45, 90, 135, 180, 225, 270, 315, 360]
            ax.set_xticks(bin_ticks)
            ax.set_xticklabels(bin_ticks)
            """
            #Can pull out medians from boxplot
            #We are computing multiple times, inefficient
            bp_bin_med = []
            for medline in bp['medians']:
                bp_bin_med.append(medline.get_ydata()[0])
            """

            #Plot the fit
            f_a = nuth_func(bin_centers, fit[0], fit[1], fit[2])
            nuth_func_str = r'$y=%0.2f*cos(%0.2f-x)+%0.2f$' % tuple(fit)
            ax.plot(bin_centers, f_a, 'b', label=nuth_func_str)

            ax.set_xlabel('Aspect (deg)')
            ax.set_ylabel('dh/tan(slope) (m)')
            ax.axhline(color='gray', linewidth=0.5)

            ax.set_xlim(*bin_range)
            ylim = ax.get_ylim()
            abs_ylim = np.max(np.abs(ylim))
            #abs_ylim = np.max(np.abs([ydata.min(), ydata.max()]))
            #pad = 0.2 * abs_ylim 
            pad = 0
            ylim = (-abs_ylim - pad, abs_ylim + pad)
            minylim = (-10,10)
            if ylim[0] > minylim[0]:
                ylim = minylim
            ax.set_ylim(*ylim)
            ax.legend(prop={'size':8})

    return fit, fit_fig
447675
@Override
	public void process(GrayF32 input , GrayU8 output ) {
		inputPow2.reshape(input.width,input.height);
		inputMean.reshape(input.width,input.height);
		inputMeanPow2.reshape(input.width,input.height);
		inputPow2Mean.reshape(input.width,input.height);
		stdev.reshape(input.width,input.height);
		tmp.reshape(input.width,input.height);
		inputPow2.reshape(input.width,input.height);

		int radius = width.computeI(Math.min(input.width,input.height))/2;

		// mean of input image = E[X]
		BlurImageOps.mean(input, inputMean, radius, tmp, work);

		// standard deviation = sqrt( E[X^2] + E[X]^2)
		PixelMath.pow2(input, inputPow2);
		BlurImageOps.mean(inputPow2,inputPow2Mean,radius,tmp, work);
		PixelMath.pow2(inputMean,inputMeanPow2);
		PixelMath.subtract(inputPow2Mean, inputMeanPow2, stdev);
		PixelMath.sqrt(stdev, stdev);

		float R = ImageStatistics.max(stdev);

		if( down ) {
			BoofConcurrency.loopFor(0, input.height, y -> {
				int i = y * stdev.width;
				int indexIn = input.startIndex + y * input.stride;
				int indexOut = output.startIndex + y * output.stride;

				for (int x = 0; x < input.width; x++, i++) {
					// threshold = mean.*(1 + k * ((deviation/R)-1));
					float threshold = inputMean.data[i] * (1.0f + k * (stdev.data[i] / R - 1.0f));
					output.data[indexOut++] = (byte) (input.data[indexIn++] <= threshold ? 1 : 0);
				}
			});
		} else {
			BoofConcurrency.loopFor(0, input.height, y -> {
				int i = y * stdev.width;
				int indexIn = input.startIndex + y * input.stride;
				int indexOut = output.startIndex + y * output.stride;

				for (int x = 0; x < input.width; x++, i++) {
					// threshold = mean.*(1 + k * ((deviation/R)-1));
					float threshold = inputMean.data[i] * (1.0f + k * (stdev.data[i] / R - 1.0f));
					output.data[indexOut++] = (byte) (input.data[indexIn++] >= threshold ? 1 : 0);
				}
			});
		}
	}
1324379
public static void equalizeLocalInner( GrayU8 input , int radius , GrayU8 output ,
										 IWorkArrays workArrays ) {

		int width = 2*radius+1;
		int area = width*width;
		int maxValue = workArrays.length()-1;

		//CONCURRENT_BELOW BoofConcurrency.loopBlocks(radius,input.height-radius,(y0,y1)->{
		int y0 = radius, y1 = input.height-radius;
		int[] histogram = workArrays.pop();
		for( int y = y0; y < y1; y++ ) {
			localHistogram(input,0,y-radius,width,y+radius+1,histogram);

			// compute equalized pixel value using the local histogram
			int inputValue = input.unsafe_get(radius, y);
			int sum = 0;
			for( int i = 0; i <= inputValue; i++ ) {
				sum += histogram[i];
			}

			output.set(radius,y, (sum*maxValue)/area );

			// start of old and new columns in histogram region
			int indexOld = input.startIndex + y*input.stride;
			int indexNew = indexOld+width;

			// index of pixel being examined
			int indexIn = input.startIndex + y*input.stride+radius+1;
			int indexOut = output.startIndex + y*output.stride+radius+1;

			for( int x = radius+1; x < input.width-radius; x++ ) {

				// update local histogram by removing the left column
				for( int i = -radius; i <= radius; i++ ) {
					histogram[input.data[indexOld + i*input.stride] & 0xFF]--;
				}

				// update local histogram by adding the right column
				for( int i = -radius; i <= radius; i++ ) {
					histogram[input.data[indexNew + i*input.stride] & 0xFF]++;
				}

				// compute equalized pixel value using the local histogram
				inputValue =  input.data[indexIn++] & 0xFF;
				sum = 0;
				for( int i = 0; i <= inputValue; i++ ) {
					sum += histogram[i];
				}

				output.data[indexOut++] = (byte)((sum*maxValue)/area);

				indexOld++;
				indexNew++;
			}
		}
		workArrays.recycle(histogram);
		//CONCURRENT_ABOVE });
	}
401793
def roi_align(featuremap, boxes, resolution):
    
    # sample 4 locations per roi bin
    ret = crop_and_resize(
        featuremap, boxes,
        tf.zeros([tf.shape(boxes)[0]], dtype=tf.int32),
        resolution * 2)
    ret = tf.nn.avg_pool(ret, [1, 1, 2, 2], [1, 1, 2, 2], padding='SAME', data_format='NCHW')
    return ret
1200471
def multilevel_roi_align(features, rcnn_boxes, resolution):
    
    assert len(features) == 4, features
    # Reassign rcnn_boxes to levels
    level_ids, level_boxes = fpn_map_rois_to_levels(rcnn_boxes)
    all_rois = []

    # Crop patches from corresponding levels
    for i, boxes, featuremap in zip(itertools.count(), level_boxes, features):
        with tf.name_scope('roi_level{}'.format(i + 2)):
            boxes_on_featuremap = boxes * (1.0 / cfg.FPN.ANCHOR_STRIDES[i])
            all_rois.append(roi_align(featuremap, boxes_on_featuremap, resolution))

    # this can fail if using TF<=1.8 with MKL build
    all_rois = tf.concat(all_rois, axis=0)  # NCHW
    # Unshuffle to the original order, to match the original samples
    level_id_perm = tf.concat(level_ids, axis=0)  # A permutation of 1~N
    level_id_invert_perm = tf.invert_permutation(level_id_perm)
    all_rois = tf.gather(all_rois, level_id_invert_perm)
    return all_rois
359963
private BufferedImage rasterizeConvoLayers(@NonNull List<INDArray> tensors3D, BufferedImage sourceImage) {
        long width = 0;
        long height = 0;

        int border = 1;
        int padding_row = 2;
        int padding_col = 80;

        /*
            We determine height of joint output image. We assume that first position holds maximum dimensionality
         */
        val shape = tensors3D.get(0).shape();
        val numImages = shape[0];
        height = (shape[2]);
        width = (shape[1]);
        //        log.info("Output image dimensions: {height: " + height + ", width: " + width + "}");
        int maxHeight = 0; //(height + (border * 2 ) + padding_row) * numImages;
        int totalWidth = 0;
        int iOffset = 1;

        Orientation orientation = Orientation.LANDSCAPE;
        /*
            for debug purposes we'll use portait only now
         */
        if (tensors3D.size() > 3) {
            orientation = Orientation.PORTRAIT;
        }



        List<BufferedImage> images = new ArrayList<>();
        for (int layer = 0; layer < tensors3D.size(); layer++) {
            INDArray tad = tensors3D.get(layer);
            int zoomed = 0;

            BufferedImage image = null;
            if (orientation == Orientation.LANDSCAPE) {
                maxHeight = (int) ((height + (border * 2) + padding_row) * numImages);
                image = renderMultipleImagesLandscape(tad, maxHeight, (int) width, (int) height);
                totalWidth += image.getWidth() + padding_col;
            } else if (orientation == Orientation.PORTRAIT) {
                totalWidth = (int) ((width + (border * 2) + padding_row) * numImages);
                image = renderMultipleImagesPortrait(tad, totalWidth, (int) width, (int) height);
                maxHeight += image.getHeight() + padding_col;
            }

            images.add(image);
        }

        if (orientation == Orientation.LANDSCAPE) {
            // append some space for arrows
            totalWidth += padding_col * 2;
        } else if (orientation == Orientation.PORTRAIT) {
            maxHeight += padding_col * 2;
            maxHeight += sourceImage.getHeight() + (padding_col * 2);
        }

        BufferedImage output = new BufferedImage(totalWidth, maxHeight, BufferedImage.TYPE_INT_RGB);
        Graphics2D graphics2D = output.createGraphics();

        graphics2D.setPaint(bgColor);
        graphics2D.fillRect(0, 0, output.getWidth(), output.getHeight());

        BufferedImage singleArrow = null;
        BufferedImage multipleArrows = null;

        /*
            We try to add nice flow arrow here
         */
        try {

            if (orientation == Orientation.LANDSCAPE) {
                try {
                    ClassPathResource resource = new ClassPathResource("arrow_sing.PNG");
                    ClassPathResource resource2 = new ClassPathResource("arrow_mul.PNG");

                    singleArrow = ImageIO.read(resource.getInputStream());
                    multipleArrows = ImageIO.read(resource2.getInputStream());
                } catch (Exception e) {
                }

                graphics2D.drawImage(sourceImage, (padding_col / 2) - (sourceImage.getWidth() / 2),
                                (maxHeight / 2) - (sourceImage.getHeight() / 2), null);

                graphics2D.setPaint(borderColor);
                graphics2D.drawRect((padding_col / 2) - (sourceImage.getWidth() / 2),
                                (maxHeight / 2) - (sourceImage.getHeight() / 2), sourceImage.getWidth(),
                                sourceImage.getHeight());

                iOffset += sourceImage.getWidth();

                if (singleArrow != null)
                    graphics2D.drawImage(singleArrow, iOffset + (padding_col / 2) - (singleArrow.getWidth() / 2),
                                    (maxHeight / 2) - (singleArrow.getHeight() / 2), null);
            } else {
                try {
                    ClassPathResource resource = new ClassPathResource("arrow_singi.PNG");
                    ClassPathResource resource2 = new ClassPathResource("arrow_muli.PNG");

                    singleArrow = ImageIO.read(resource.getInputStream());
                    multipleArrows = ImageIO.read(resource2.getInputStream());
                } catch (Exception e) {
                }

                graphics2D.drawImage(sourceImage, (totalWidth / 2) - (sourceImage.getWidth() / 2),
                                (padding_col / 2) - (sourceImage.getHeight() / 2), null);

                graphics2D.setPaint(borderColor);
                graphics2D.drawRect((totalWidth / 2) - (sourceImage.getWidth() / 2),
                                (padding_col / 2) - (sourceImage.getHeight() / 2), sourceImage.getWidth(),
                                sourceImage.getHeight());

                iOffset += sourceImage.getHeight();
                if (singleArrow != null)
                    graphics2D.drawImage(singleArrow, (totalWidth / 2) - (singleArrow.getWidth() / 2),
                                    iOffset + (padding_col / 2) - (singleArrow.getHeight() / 2), null);

            }
            iOffset += padding_col;
        } catch (Exception e) {
            // if we can't load images - ignore them
        }



        /*
            now we merge all images into one big image with some offset
        */


        for (int i = 0; i < images.size(); i++) {
            BufferedImage curImage = images.get(i);
            if (orientation == Orientation.LANDSCAPE) {
                // image grows from left to right
                graphics2D.drawImage(curImage, iOffset, 1, null);
                iOffset += curImage.getWidth() + padding_col;

                if (singleArrow != null && multipleArrows != null) {
                    if (i < images.size() - 1) {
                        // draw multiple arrows here
                        if (multipleArrows != null)
                            graphics2D.drawImage(multipleArrows,
                                            iOffset - (padding_col / 2) - (multipleArrows.getWidth() / 2),
                                            (maxHeight / 2) - (multipleArrows.getHeight() / 2), null);
                    } else {
                        // draw single arrow
                        //    graphics2D.drawImage(singleArrow, iOffset - (padding_col / 2) - (singleArrow.getWidth() / 2), (maxHeight / 2) - (singleArrow.getHeight() / 2), null);
                    }
                }
            } else if (orientation == Orientation.PORTRAIT) {
                // image grows from top to bottom
                graphics2D.drawImage(curImage, 1, iOffset, null);
                iOffset += curImage.getHeight() + padding_col;

                if (singleArrow != null && multipleArrows != null) {
                    if (i < images.size() - 1) {
                        // draw multiple arrows here
                        if (multipleArrows != null)
                            graphics2D.drawImage(multipleArrows, (totalWidth / 2) - (multipleArrows.getWidth() / 2),
                                            iOffset - (padding_col / 2) - (multipleArrows.getHeight() / 2), null);
                    } else {
                        // draw single arrow
                        //   graphics2D.drawImage(singleArrow, (totalWidth / 2) - (singleArrow.getWidth() / 2),  iOffset - (padding_col / 2) - (singleArrow.getHeight() / 2) , null);
                    }
                }
            }
        }

        return output;
    }
484198
private void process(int[] srcImage, int[] dstImage, int x, int y)
    {
        final int a = getSourcePixel(srcImage, x - 1, y - 1);
        final int b = getSourcePixel(srcImage, x, y - 1);
        final int c = getSourcePixel(srcImage, x + 1, y - 1);
        final int d = getSourcePixel(srcImage, x - 1, y);
        final int e = getSourcePixel(srcImage, x, y);
        final int f = getSourcePixel(srcImage, x + 1, y);
        final int g = getSourcePixel(srcImage, x - 1, y + 1);
        final int h = getSourcePixel(srcImage, x, y + 1);
        final int i = getSourcePixel(srcImage, x + 1, y + 1);
        int e0 = e;
        int e1 = e;
        int e2 = e;
        int e3 = e;
        int e4 = e;
        int e5 = e;
        int e6 = e;
        int e7 = e;
        int e8 = e;

        if (b != h && d != f)
        {
            e0 = computeE0(b, d, e);
            e1 = computeE1(a, b, c, d, e, f);
            e2 = computeE2(b, e, f);
            e3 = computeE3(a, b, d, e, g, h);
            e4 = e;
            e5 = computeE5(b, c, e, f, h, i);
            e6 = computeE6(d, e, h);
            e7 = computeE7(d, e, f, g, h, i);
            e8 = computeE8(e, f, h);
        }

        setDestPixel(dstImage, x * SCALE, y * SCALE, e0);
        setDestPixel(dstImage, x * SCALE + 1, y * SCALE, e1);
        setDestPixel(dstImage, x * SCALE + 2, y * SCALE, e2);
        setDestPixel(dstImage, x * SCALE, y * SCALE + 1, e3);
        setDestPixel(dstImage, x * SCALE + 1, y * SCALE + 1, e4);
        setDestPixel(dstImage, x * SCALE + 2, y * SCALE + 1, e5);
        setDestPixel(dstImage, x * SCALE, y * SCALE + 2, e6);
        setDestPixel(dstImage, x * SCALE + 1, y * SCALE + 2, e7);
        setDestPixel(dstImage, x * SCALE + 2, y * SCALE + 2, e8);
    }
777607
public static void equalizeLocalNaive( GrayU8 input , int radius , GrayU8 output ,
										   IWorkArrays workArrays )
	{
		int width = 2*radius+1;
		int maxValue = workArrays.length()-1;

		//CONCURRENT_BELOW BoofConcurrency.loopBlocks(0,input.height,(idx0,idx1)->{
		int idx0 = 0, idx1 = input.height;
		int[] histogram = workArrays.pop();
		for( int y = idx0; y < idx1; y++ ) {
			// make sure it's inside the image bounds
			int y0 = y-radius;
			int y1 = y+radius+1;
			if( y0 < 0 ) {
				y0 = 0; y1 = width;
				if( y1 > input.height )
					y1 = input.height;
			} else if( y1 > input.height ) {
				y1 = input.height;
				y0 = y1 - width;
				if( y0 < 0 )
					y0 = 0;
			}

			// pixel indexes
			int indexIn = input.startIndex + y*input.stride;
			int indexOut = output.startIndex + y*output.stride;

			for( int x = 0; x < input.width; x++ ) {
				// make sure it's inside the image bounds
				int x0 = x-radius;
				int x1 = x+radius+1;
				if( x0 < 0 ) {
					x0 = 0; x1 = width;
					if( x1 > input.width )
						x1 = input.width;
				} else if( x1 > input.width ) {
					x1 = input.width;
					x0 = x1 - width;
					if( x0 < 0 )
						x0 = 0;
				}

				// compute the local histogram
				localHistogram(input,x0,y0,x1,y1,histogram);

				// only need to compute up to the value of the input pixel
				int inputValue =  input.data[indexIn++] & 0xFF;
				int sum = 0;
				for( int i = 0; i <= inputValue; i++ ) {
					sum += histogram[i];
				}

				int area = (y1-y0)*(x1-x0);
				output.data[indexOut++] = (byte)((sum*maxValue)/area);
			}
		}
		workArrays.recycle(histogram);
		//CONCURRENT_ABOVE });
	}
705357
public static void horizontal(Kernel1D_F64 kernel, GrayF64 src, GrayF64 dst ) {
		InputSanityCheck.checkSameShape(src, dst);

		boolean processed = BOverrideConvolveImageNormalized.invokeNativeHorizontal(kernel,src,dst);
		
		if( !processed ) {
			if( kernel.width >= src.width ) {
				ConvolveNormalizedNaive_SB.horizontal(kernel,src,dst);
			} else {
				if( Math.abs(kernel.computeSum() - 1.0f) > 1e-4f ) {
					Kernel1D_F64 k = kernel.copy();
					KernelMath.normalizeSumToOne(k);
					kernel = k;
				}
				ConvolveImageNoBorder.horizontal(kernel,src,dst);
				ConvolveNormalized_JustBorder_SB.horizontal(kernel,src,dst);
			}
		}
	}
644072
private void horizontallyFromSrcToWorkGray(BufferedImage srcImg, byte[][] workPixels, int start, int delta) {
    final
    int[]
      tempPixels =
      new int[srcWidth];   // Used if we work on int based bitmaps, later used to keep channel values
    final byte[] srcPixels = new byte[srcWidth]; // create reusable row to minimize memory overhead

    for (int k = start; k < srcHeight; k = k + delta) {
      ImageUtils.getPixelsBGR(srcImg, k, srcWidth, srcPixels, tempPixels);

      for (int i = dstWidth - 1; i >= 0; i--) {
        final int max = horizontalSubsamplingData.arrN[i];

        float sample0 = 0.0f;
        int index = i * horizontalSubsamplingData.numContributors;
        for (int j = max - 1; j >= 0; j--) {
          float arrWeight = horizontalSubsamplingData.arrWeight[index];
          int pixelIndex = horizontalSubsamplingData.arrPixel[index];

          sample0 += (srcPixels[pixelIndex] & 0xff) * arrWeight;
          index++;
        }

        workPixels[k][i] = toByte(sample0);
      }
    }
  }
195468
public static void horizontal(GrayF32 input, GrayF32 output, int radius) {

		boolean processed = BOverrideConvolveImageMean.invokeNativeHorizontal(input, output, radius);

		if( !processed ) {
			Kernel1D_F32 kernel = FactoryKernel.table1D_F32(radius, true);
			if (kernel.width > input.width) {
				ConvolveImageNormalized.horizontal(kernel, input, output);
			} else {
				InputSanityCheck.checkSameShape(input, output);
				ConvolveNormalized_JustBorder_SB.horizontal(kernel, input, output);
				if(BoofConcurrency.USE_CONCURRENT) {
					ImplConvolveMean_MT.horizontal(input, output, radius);
				} else {
					ImplConvolveMean.horizontal(input, output, radius);
				}
			}
		}
	}