start new max hacking

TODO

@@ -1,4 +1,3 @@
-- check libtool version number, should be binary-compat with 7.30 
 
 - quadratic doesn't work for order 3
 
@@ -11,7 +10,11 @@
 
 - seems to leak with vips_check_matrix()?
 
-- test with interpolators
+- test with interpolators, do we add margins correctly?
+
+
+
+- check libtool version number, should be binary-compat with 7.30 
 
 - the operation cache needs to detect invalidate
 

libvips/arithmetic/max.c

@@ -68,21 +68,100 @@
 
 #include "statistic.h"
 
+/* Track max values and position here. We need one of these for each thread,
+ * and one for the main value.
+ *
+ * We will generally only be tracking a small (<10?) number of values, so
+ * simple arrays will be fastest.
+ */
+typedef struct _VipsValues {
+	struct _VipsMax *max;
+
+	/* The max number of values we track.
+	 */
+	int size;
+
+	/* How many values we have in the arrays.
+	 */
+	int n;
+
+	/* Position and values. We track mod**2 for complex and do a sqrt() at
+	 * the end. The three arrays are sorted by values, smallest first.
+	 */
+	double *value;
+	int *x_pos;
+	int *y_pos;
+} VipsValues;
+
 typedef struct _VipsMax {
 	VipsStatistic parent_instance;
 
-	gboolean set;		/* FALSE means no value yet */
-
-	/* The current maximum. When scanning complex images, we keep the
-	 * square of the modulus here and do a single sqrt() right at the end.
+	/* Max number of values we track.
 	 */
-	double max;
+	int max_values;
 
-	/* And its position.
+	/* Global state here.
 	 */
-	int x, y;
+	VipsValues values;
 } VipsMax;
 
+static void
+vips_values_init( VipsValues *values, VipsMax *max )
+{
+	values->max = max;
+
+	values->size = max->max_values;
+	values->n = 0;
+	values->value = VIPS_ARRAY( max, values->size, double );
+	values->x_pos = VIPS_ARRAY( max, values->size, int );
+	values->y_pos = VIPS_ARRAY( max, values->size, int );
+}
+
+/* Add a value. Do nothing if the value is too small.
+ */
+static void
+vips_values_add( VipsValues *values, double v, int x, int y )
+{
+	int i, j;
+
+	/* Find insertion point.
+	 */
+	for( i = 0; i < values->n; i++ )
+		if( values->value[i] > v ) 
+			break;
+
+	/* Array full? 
+	 */
+	if( values->n == values->size ) {
+		if( i > 0 ) {
+			/* We need to move stuff to the left to make space,
+			 * shunting the smallest out.
+			 */
+			for( j = 0; j < i - 1; j++ ) {
+				values->value[j] = values->value[j + 1];
+				values->x_pos[j] = values->x_pos[j + 1];
+				values->y_pos[j] = values->y_pos[j + 1];
+			}
+			values->value[i - 1] = v;
+			values->x_pos[i - 1] = x;
+			values->y_pos[i - 1] = y;
+		}
+	}
+	else {
+		/* Not full, move stuff to the right into empty space.
+		 */
+		for( j = values->n + 1; j > i; j-- ) {
+			values->value[j] = values->value[j - 1];
+			values->x_pos[j] = values->x_pos[j - 1];
+			values->y_pos[j] = values->y_pos[j - 1];
+		}
+		values->value[i] = v;
+		values->x_pos[i] = x;
+		values->y_pos[i] = y;
+		values->n += 1;
+	}
+}
+
 typedef VipsStatisticClass VipsMaxClass;
 
 G_DEFINE_TYPE( VipsMax, vips_max, VIPS_TYPE_STATISTIC );
@@ -93,6 +172,8 @@ vips_max_build( VipsObject *object )
 	VipsStatistic *statistic = VIPS_STATISTIC( object ); 
 	VipsMax *max = (VipsMax *) object;
 
+	vips_values_init( *max->values, max );
+
 	if( VIPS_OBJECT_CLASS( vips_max_parent_class )->build( object ) )
 		return( -1 );
 
@@ -122,17 +203,17 @@ vips_max_build( VipsObject *object )
 	return( 0 );
 }
 
-/* New sequence value. Make a private VipsMax for this thread.
+/* New sequence value. Make a private VipsValues for this thread.
  */
 static void *
 vips_max_start( VipsStatistic *statistic )
 {
-	VipsMax *max;
+	VipsValues *values;
 
-	max = g_new( VipsMax, 1 );
-	max->set = FALSE;
+	values = g_new( VipsValues, 1 );
+	vips_values_init( values, (VipsMax *) statistic ); 
 
-	return( (void *) max );
+	return( (void *) values );
 }
 
 /* Merge the sequence value back into the per-call state.
@@ -140,18 +221,16 @@ vips_max_start( VipsStatistic *statistic )
 static int
 vips_max_stop( VipsStatistic *statistic, void *seq )
 {
-	VipsMax *global = (VipsMax *) statistic;
-	VipsMax *local = (VipsMax *) seq;
+	VipsMax *max = (VipsMax *) statistic;
+	VipsValues *values = (VipsValues *) seq;
 
-	if( local->set &&
-		(!global->set || local->max > global->max) ) {
-		global->max = local->max;
-		global->x = local->x;
-		global->y = local->y;
-		global->set = TRUE;
-	}
+	int i;
 
-	g_free( local );
+	for( i = 0; i < values->n; i++ )
+		vips_values_add( &max->values, 
+			values->value[i], values->x_pos[i], values->y_pos[i] );
+
+	g_free( values );
 
 	return( 0 );
 }

libvips/resample/quadratic.c

@@ -68,6 +68,20 @@ y = y' + g
   + j y' x'             
   + k y' y' + l x' x'  
 
+input matrix:
+
+  a g
+  --
+  b h
+  c i
+  --
+  d j
+  --
+  e k
+  f l
+
+matrix height may be 1, 3, 4, 6
+
  */
 
 typedef struct _VipsQuadratic {
@@ -134,19 +148,21 @@ vips_quadratic_gen( VipsRegion *or, void *vseq,
 		return( -1 );
 
 	for( yo = ylow; yo < yhigh; yo++ ) {
-		fxi = xlow + vec[0];                /* order 0 */
-		fyi = yo + vec[1];    
-		dx = 1.0;
+		fxi = 0.0;
+		fyi = 0.0;
+		dx = 0.0;
 		dy = 0.0;
+		ddx = 0.0;
+		ddy = 0.0;
 
 		switch( quadratic->order ) {
 		case 3: 
 			fxi += vec[10] * yo * yo + vec[8] * xlow * xlow;
 			fyi += vec[11] * yo * yo + vec[9] * xlow * xlow;
 			dx += vec[8];
-			ddx = vec[8] * 2.0;
+			ddx += vec[8] * 2.0;
 			dy += vec[9];
-			ddy = vec[9] * 2.0;
+			ddy += vec[9] * 2.0;
 
 		case 2: 
 			fxi += vec[6] * xlow * yo;
@@ -161,15 +177,17 @@ vips_quadratic_gen( VipsRegion *or, void *vseq,
 			dy += vec[3];
 
 		case 0: 
-			/* See above for order 0.
-			 */
+			fxi += vec[0];
+			fyi += vec[1];    
 			break;
 
 		default:
 		    	g_assert( 0 );
-		    	return(-7);
+		    	return( -7 );
 		}
 
+		printf( "dx = %g, dy = %g\n", dx, dy );
+
 		q = VIPS_REGION_ADDR( or, xlow, yo );
 
 		for( xo = xlow; xo < xhigh; xo++ ) {
@@ -261,7 +279,7 @@ vips_quadratic_build( VipsObject *object )
 
 	if( !vips_object_argument_isset( object, "interpolator" ) )
 		g_object_set( object, 
-			"interpolator", vips_interpolate_new( "bilinear" ),
+			"interpolate", vips_interpolate_new( "bilinear" ),
 			NULL );
 	interpolate = quadratic->interpolate;