add CLAHE

vips_hist_local() has a new param, max_slope, which sets the maximum
amount that the local contrast can be boosted by -- CLAHE

see https://github.com/jcupitt/libvips/issues/268

also fix a couple of small range problems

- scale by 255, not 256, to avoid an overflow
- cum hist includes the current value, so sum to <= target, not < target

ChangeLog

@@ -27,6 +27,7 @@
   aferrero2707
 - kick load operations from cache on read error, thanks gaillard 
 - fix return from C++ assignment operator overloads (+=, -= etc) 
+- add @max_slope to vips_local_hist() to implement CLAHE
 
 8/12/16 started 8.4.5
 - allow libgsf-1.14.26 to help centos, thanks tdiprima

libvips/histogram/hist_local.c

@@ -23,6 +23,9 @@
  * 	- any number of bands
  * 20/1/17
  * 	- add contrast limit
+ * 	- sum to <= target, not < target, since cumulative hists include the
+ * 	  current value
+ * 	- scale result by 255, not 256, to avoid overflow
  */
 
 /*
@@ -73,15 +76,7 @@ typedef struct _VipsHistLocal {
 	int width;
 	int height;
 
-	double clip_limit;
-
-	/* width * height, ie. the area of the histogram we make.
-	 */
-	int area;
-
-	/* clip_limit * area, ie. relative to our maximum height.
-	 */
-	int clip;
+	int max_slope;
 
 } VipsHistLocal;
 
@@ -189,11 +184,10 @@ vips_hist_local_generate( VipsRegion *or,
 			memset( seq->hist[b], 0, 256 * sizeof( unsigned int ) );
 		p1 = p;
 		for( j = 0; j < local->height; j++ ) {
-			i = 0; 
-			for( x = 0; x < local->width; x++ )
-				for( b = 0; b < bands; b++ )
-					seq->hist[b][p1[i++]] += 1;
-			
+			for( i = 0, x = 0; x < local->width; x++ )
+				for( b = 0; b < bands; b++, i++ )
+					seq->hist[b][p1[i]] += 1;
+
 			p1 += lsk;
 		}
 
@@ -204,7 +198,8 @@ vips_hist_local_generate( VipsRegion *or,
 				/* Sum histogram up to current pel.
 				 */
 				unsigned int *hist = seq->hist[b]; 
-				int target = p[centre];
+				const int target = p[centre];
+				const int max_slope = local->max_slope;
 
 				int sum;
 
@@ -212,31 +207,52 @@ vips_hist_local_generate( VipsRegion *or,
 
 				/* For CLAHE we need to limit the height of the
 				 * hist to limit the amount we boost the
-				 * contrast by. If there's no limit, we can
-				 * take a shorter path.
+				 * contrast by. 
 				 */
-				if( local->clip_limit < 1.0 ) { 
+				if( max_slope > 0 ) {
 					int sum_over;
 
 					sum_over = 0;
-					for( i = 0; i < 256; i++ ) {
-						if( hist[i] > local->clip ) 
-							sum_over += hist[i];
+
+					/* Must be <= target, since a cum hist
+					 * always includes the current element.
+					 */
+					for( i = 0; i <= target; i++ ) {
+						if( hist[i] > max_slope ) {
+							sum_over += hist[i] - 
+								max_slope;
+							sum += max_slope;
+						}
 						else 
 							sum += hist[i];
 					}
 
-					sum += sum_over / 256;
+					for( ; i < 256; i++ ) {
+						if( hist[i] > max_slope ) 
+							sum_over += hist[i] - 
+								max_slope;
+					}
 
+					/* The extra clipped off bit from the
+					 * top of the hist is spread over all
+					 * bins equally, then summed to target.
+					 */
+					sum += (target + 1) * sum_over / 256;
 				}
 				else {
 					sum = 0;
-					for( i = 0; i < target; i++ )
+					for( i = 0; i <= target; i++ )
 						sum += hist[i];
-				
 				}
 
-				*q++ = 256 * sum / local->area;
+				/* This can't overflow, even in
+				 * contrast-limited mode.
+				 *
+				 * Scale by 255, not 256, or we'll get
+				 * overflow.
+				 */
+				*q++ = 255 * sum / 
+					(local->width * local->height);
 
 				/* Adapt histogram --- remove the pels from 
 				 * the left hand column, add in pels for a 
@@ -285,10 +301,6 @@ vips_hist_local_build( VipsObject *object )
 		return( -1 );
 	}
 
-	local->area = local->width * local->height;
-
-	local->clip = local->clip_limit * local->area;
-
 	/* Expand the input. 
 	 */
 	if( vips_embed( in, &t[1], 
@@ -364,19 +376,18 @@ vips_hist_local_class_init( VipsHistLocalClass *class )
 		G_STRUCT_OFFSET( VipsHistLocal, height ),
 		1, VIPS_MAX_COORD, 1 );
 
-	VIPS_ARG_DOUBLE( class, "clip_limit", 6, 
-		_( "Clip limit" ), 
-		_( "Clip the histogram at this height (CLAHE)" ),
+	VIPS_ARG_INT( class, "max_slope", 6, 
+		_( "Max slope" ), 
+		_( "Maximum slope (CLAHE)" ),
 		VIPS_ARGUMENT_OPTIONAL_INPUT,
-		G_STRUCT_OFFSET( VipsHistLocal, clip_limit ),
-		0.0, 1.0, 1.0 );
+		G_STRUCT_OFFSET( VipsHistLocal, max_slope ),
+		0, 100, 0 );
 
 }
 
 static void
 vips_hist_local_init( VipsHistLocal *local )
 {
-	local->clip_limit = 1.0;
 }
 
 /**
@@ -387,11 +398,20 @@ vips_hist_local_init( VipsHistLocal *local )
  * @height: height of region
  * @...: %NULL-terminated list of optional named arguments
  *
+ * Optional arguments:
+ *
+ * * @max_slope: maximum brightening
+ *
  * Performs local histogram equalisation on @in using a
- * window of size @xwin by @ywin centered on the input pixel. 
+ * window of size @width by @height centered on the input pixel. 
  *
  * The output image is the same size as the input image. The edge pixels are
- * created by copy edge pixels of the input image outwards.
+ * created by mirroring the input image outwards.
+ *
+ * If @max_slope is greater than 0, it sets the maximum value for the slope of
+ * the cumulative histogram, that is, the maximum brightening that is
+ * performed. A value of 3 is often used. Local histogram equalization with
+ * contrast limiting is usually called CLAHE.
  *
  * See also: vips_hist_equal().
  *