Computes two independent 1D histograms from hue and saturation.  Less affects by sparsity, but can produce
worse results since the basic assumption that hue and saturation are decoupled is most of the time false.
773712
public static List<double[]> independentHueSat( List<File> images  ) {
		List<double[]> points = new ArrayList<>();

		// The number of bins is an important parameter.  Try adjusting it
		TupleDesc_F64 histogramHue = new TupleDesc_F64(30);
		TupleDesc_F64 histogramValue = new TupleDesc_F64(30);

		List<TupleDesc_F64> histogramList = new ArrayList<>();
		histogramList.add(histogramHue); histogramList.add(histogramValue);

		Planar<GrayF32> rgb = new Planar<>(GrayF32.class,1,1,3);
		Planar<GrayF32> hsv = new Planar<>(GrayF32.class,1,1,3);

		for( File f : images ) {
			BufferedImage buffered = UtilImageIO.loadImage(f.getPath());
			if( buffered == null ) throw new RuntimeException("Can't load image!");

			rgb.reshape(buffered.getWidth(), buffered.getHeight());
			hsv.reshape(buffered.getWidth(), buffered.getHeight());
			ConvertBufferedImage.convertFrom(buffered, rgb, true);
			ColorHsv.rgbToHsv(rgb, hsv);

			GHistogramFeatureOps.histogram(hsv.getBand(0), 0, 2*Math.PI,histogramHue);
			GHistogramFeatureOps.histogram(hsv.getBand(1), 0, 1, histogramValue);

			// need to combine them into a single descriptor for processing later on
			TupleDesc_F64 imageHist = UtilFeature.combine(histogramList,null);

			UtilFeature.normalizeL2(imageHist); // normalize so that image size doesn't matter

			points.add(imageHist.value);
		}

		return points;
	}
1191713
def mix(color1, color2, pos=0.5):
    
    opp_pos = 1 - pos

    red = color1[0] * pos + color2[0] * opp_pos
    green = color1[1] * pos + color2[1] * opp_pos
    blue = color1[2] * pos + color2[2] * opp_pos
    return int(red), int(green), int(blue)
1045067
private static float hue2rgb(float m1, float m2, float h) {
        if (h < 0) h += 1;
        if (h > 1) h -= 1;
        if (h * 6f < 1) return m1 + (m2 - m1) * 6f * h;
        if (h * 2f < 1) return m2;
        if (h * 3f < 2) return m1 + (m2 - m1) * (2 / 3f - h) * 6f;
        return m1;
    }
459985
function(color1, color2, factor) {
		if (arguments.length < 3) { factor = 0.5; }
		var result = color1.slice();
		for (var i=0;i<3;i++) {
			result[i] = Math.round(result[i] + factor*(color2[i]-color1[i]));
		}
		return result;
	}
9638
public void computeHistogram( T image , RectangleRotate_F32 region ) {

		interpolate.setImage(image);

		c = (float)Math.cos( region.theta );
		s = (float)Math.sin( region.theta );

		for( int i = 0; i < histogram.length; i++ ) {
			histogram[i] = 0;
		}

		// if it is entirely inside the image, interpolate using a faster technique
		if( isInFastBounds(region) ) {
			computeHistogramInside( region);
		} else {
			computeHistogramBorder(image, region);
		}

		normalizeHistogram();
	}
68863
public void process( T image1 , T image2 , ImageFlow output) {

		InputSanityCheck.checkSameShape(image1,image2);

		derivX.reshape(image1.width,image1.height);
		derivY.reshape(image1.width,image1.height);
		derivT.reshape(image1.width,image1.height);

		averageFlow.reshape(output.width,output.height);

		if( resetOutput )
			output.fillZero();

		computeDerivX(image1,image2,derivX);
		computeDerivY(image1,image2,derivY);
		computeDerivT(image1,image2,derivT);

		findFlow(derivX,derivY,derivT,output);
	}
211643
def get_color_dist(c1, c2):
    
    return tuple(abs(v1 - v2) for v1, v2 in itertools.izip(c1, c2))
947672
def hue(p):
    
    min_c = min(p)
    max_c = max(p)
    d = float(max_c - min_c)
    if d == 0:
        return 0

    if max_c == p[0]:
        h = (p[1] - p[2]) / d
    elif max_c == p[1]:
        h = 2 + (p[2] - p[0]) / d
    else:
        h = 4 + (p[0] - p[1]) / d
    h *= 60
    if h < 0:
        h += 360
    return h
1004739
function hsl(hue, saturation, lightness){
  if (1 == arguments.length) {
    utils.assertColor(hue, 'color');
    return hue.hsla;
  } else {
    return hsla(
        hue
      , saturation
      , lightness
      , new nodes.Unit(1));
  }
}
64370
private static int blendColors(int color1, int color2, float ratio) {
    final float inverseRation = 1f - ratio;
    float r = (Color.red(color1) * ratio) + (Color.red(color2) * inverseRation);
    float g = (Color.green(color1) * ratio) + (Color.green(color2) * inverseRation);
    float b = (Color.blue(color1) * ratio) + (Color.blue(color2) * inverseRation);
    return Color.rgb((int) r, (int) g, (int) b);
  }
378642
private static double calculateHue(double red, double green, double blue, double max, double chroma) {

    if (chroma == 0) {
      return 0;
    } else {
      double hue;
      if (red == max) {
        hue = (green - blue) / chroma;
      } else if (green == max) {
        hue = (blue - red) / chroma + 2;
      } else {
        hue = (red - green) / chroma + 4;
      }
      hue = hue * 60.0;
      if (hue < 0) {
        hue = hue + Hue.MAX_VALUE;
      }
      return hue;
    }
  }
513630
def mix!(other, pct = 50.0)
      coeff = pct.to_f / 100.0
      new_rgb = to_rgb.zip(other.to_rgb).map{|c1, c2| (c1 * (1 - coeff)) + (c2 * coeff)}
      h, s, l = self.class.rgb_to_hsl(*new_rgb)
      self.hue, self.saturation, self.lightness = 360*h, s, l
    end
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
1054724
def compute_offset_sad(dem1, dem2, pad=(9,9), plot=False):
    
    #This defines the search window size
    #Use half-pixel stride?
    #Note: stride is not properly implemented 
    #stride = 1
    #ref = dem1[::stride,::stride]
    #kernel = dem2[pad[0]:-pad[0]:stride, pad[1]:-pad[1]:stride]
    kernel = dem2[pad[0]:-pad[0], pad[1]:-pad[1]]
    #Want to pad evenly on both sides, so add +1 here
    m = np.zeros((pad[0]*2+1, pad[1]*2+1))
   
    #Find integer pixel offset
    i = j = 0
    for i in range(m.shape[0]):
        print(i)
        for j in range(m.shape[1]):
            print(j)
            ref = dem1[i:i+kernel.shape[0], j:j+kernel.shape[1]]
            diff = ref - kernel
            
            #Remove outliers beyond IQR
            diff_iqr = malib.calcperc(diff, (25,75))
            diff = np.ma.masked_outside(diff, *diff_iqr)
            """ 
            diff_med = np.ma.median(diff)
            diff_mad = malib.mad(diff)
            diff_madr = (diff_med - mad, diff_med + mad)
            diff = np.ma.masked_outside(diff, diff_madr)     
            """
            #Masked areas will decrease sum! Normalize by count of valid pixels
            m[i,j] = np.ma.abs(diff).sum()/diff.count()
    
    #Note, we're dealing with min SAD here, so want to provide -m for sub-pixel refinement 
    m = -m  

    int_argmax = np.array(np.unravel_index(m.argmax(), m.shape))
    int_offset = int_argmax - pad
    
    sp_argmax = np.array(find_subpixel_peak_position(m, 'parabolic'))
    sp_offset = sp_argmax - pad

    if plot:
        plt.figure()
        plt.title('Sum of Absolute Differences')
        plt.imshow(m)
        plt.scatter(*sp_argmax[::-1])
        #plt.show()

    return m, int_offset, sp_offset
1233400
void computeHistogram(int c_x, int c_y, double sigma) {
		int r = (int)Math.ceil(sigma * sigmaEnlarge);

		// specify the area being sampled
		bound.x0 = c_x - r;
		bound.y0 = c_y - r;
		bound.x1 = c_x + r + 1;
		bound.y1 = c_y + r + 1;

		ImageGray rawDX = derivX.getImage();
		ImageGray rawDY = derivY.getImage();

		// make sure it is contained in the image bounds
		BoofMiscOps.boundRectangleInside(rawDX,bound);

		// clear the histogram
		Arrays.fill(histogramMag,0);
		Arrays.fill(histogramX,0);
		Arrays.fill(histogramY,0);

		// construct the histogram
		for( int y = bound.y0; y < bound.y1; y++ ) {
			// iterate through the raw array for speed
			int indexDX = rawDX.startIndex + y*rawDX.stride + bound.x0;
			int indexDY = rawDY.startIndex + y*rawDY.stride + bound.x0;

			for( int x = bound.x0; x < bound.x1; x++ ) {
				float dx = derivX.getF(indexDX++);
				float dy = derivY.getF(indexDY++);

				// edge intensity and angle
				double magnitude = Math.sqrt(dx*dx + dy*dy);
				double theta = UtilAngle.domain2PI(Math.atan2(dy,dx));

				// weight from gaussian
				double weight = computeWeight( x-c_x, y-c_y , sigma );

				// histogram index
				int h = (int)(theta / histAngleBin) % histogramMag.length;

				// update the histogram
				histogramMag[h] += magnitude*weight;
				histogramX[h] += dx*weight;
				histogramY[h] += dy*weight;
			}
		}
	}
562692
public static void histogram( GrayF32 input , float minValue , int histogram[] ) {
		if( BoofConcurrency.USE_CONCURRENT ) {
			ImplImageStatistics_MT.histogram(input,minValue,histogram);
		} else {
			ImplImageStatistics.histogram(input,minValue,histogram);
		}
	}
100471
function inerpolateColors(col1, col2, amount) {
  let col1rgb = parseColor(col1);
  let col2rgb = parseColor(col2);
  return RGBobj2string({
    r: Math.min(255, col1rgb.r * amount + col2rgb.r * (1-amount)),
    g: Math.min(255, col1rgb.g * amount + col2rgb.g * (1-amount)),
    b: Math.min(255, col1rgb.b * amount + col2rgb.b * (1-amount)),
    a: Math.min(1, col1rgb.a * amount + col2rgb.a * (1-amount))
  });
}