From ce6d37189c3a2e217d405882ac14a31539dd5b2e Mon Sep 17 00:00:00 2001 From: John Cupitt Date: Sat, 3 Feb 2018 17:34:49 +0000 Subject: [PATCH] non-interpolatory edge thinning works woo --- libvips/convolution/angle_interpolate_canny.c | 388 ++++++++++++++++++ libvips/convolution/canny.c | 59 +-- 2 files changed, 423 insertions(+), 24 deletions(-) create mode 100644 libvips/convolution/angle_interpolate_canny.c diff --git a/libvips/convolution/angle_interpolate_canny.c b/libvips/convolution/angle_interpolate_canny.c new file mode 100644 index 00000000..4a964020 --- /dev/null +++ b/libvips/convolution/angle_interpolate_canny.c @@ -0,0 +1,388 @@ +/* Canny edge detector + * + * 2/2/18 + * - from vips_canny() + */ + +/* + + This file is part of VIPS. + + VIPS is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA + 02110-1301 USA + + */ + +/* + + These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk + + */ + +/* +#define DEBUG + */ + +#ifdef HAVE_CONFIG_H +#include +#endif /*HAVE_CONFIG_H*/ +#include + +#include +#include +#include + +#include + +typedef struct _VipsCanny { + VipsOperation parent_instance; + + VipsImage *in; + VipsImage *out; + + gdouble sigma; + + /* Need an image vector for start_many. + */ + VipsImage *args[3]; +} VipsCanny; + +typedef VipsOperationClass VipsCannyClass; + +G_DEFINE_TYPE( VipsCanny, vips_canny, VIPS_TYPE_OPERATION ); + +static int +vips_canny_polar_generate( VipsRegion *or, + void *vseq, void *a, void *b, gboolean *stop ) +{ + VipsRegion **in = (VipsRegion **) vseq; + VipsRect *r = &or->valid; + VipsImage *Gx = in[0]->im; + + int x, y, band; + + if( vips_reorder_prepare_many( or->im, in, r ) ) + return( -1 ); + + for( y = 0; y < r->height; y++ ) { + VipsPel *p1 = (VipsPel * restrict) + VIPS_REGION_ADDR( in[0], r->left, r->top + y ); + VipsPel *p2 = (VipsPel * restrict) + VIPS_REGION_ADDR( in[1], r->left, r->top + y ); + VipsPel *q = (VipsPel * restrict) + VIPS_REGION_ADDR( or, r->left, r->top + y ); + + for( x = 0; x < r->width; x++ ) { + for( band = 0; band < Gx->Bands; band++ ) { + int x = p1[band] - 128; + int y = p2[band] - 128; + + /* Faster than hypot() for int args. Scale down + * or we'll clip on very hard edgaes. + * + * Code theta as 0 - 255 for 0 - 360. + */ + q[0] = 0.5 * sqrt( x * x + y * y ); + q[1] = 255.0 * atan2( x, y ) / (2 * VIPS_PI); + + q += 2; + } + + p1 += Gx->Bands; + p2 += Gx->Bands; + } + } + + return( 0 ); +} + +/* Calculate G/theta from Gx/Gy -- rather like rect -> polar, except that we + * code theta as 0 - 255 for 0 - 360 degrees. Scale G down by 0.5 so that we + * don't clip on hard edges. + * + * For a white disc on a black background, the top edge will be 128, right 64, + * bottom 255/0, left 192. + */ +static int +vips_canny_polar( VipsImage **args, VipsImage **out ) +{ + *out = vips_image_new(); + if( vips_image_pipeline_array( *out, + VIPS_DEMAND_STYLE_THINSTRIP, args ) ) + return( -1 ); + (*out)->Bands *= 2; + + if( vips_image_generate( *out, + vips_start_many, vips_canny_polar_generate, vips_stop_many, + args, NULL ) ) + return( -1 ); + + return( 0 ); +} + +static int +vips_canny_nonmax_generate( VipsRegion *or, + void *vseq, void *a, void *b, gboolean *stop ) +{ + VipsRegion *in = (VipsRegion *) vseq; + VipsRect *r = &or->valid; + VipsImage *im = in->im; + int out_bands = or->im->Bands; + + VipsRect rect; + int x, y, band; + int lsk; + int psk; + + /* For each of the 256 angles, the offsets to the four neighbours we + * fetch, and the two weights we use. + */ + int offset_lowa[256]; + int offset_lowb[256]; + int offset_higha[256]; + int offset_highb[256]; + int weighta[256]; + int weightb[256]; + + rect = *r; + rect.width += 2; + rect.height += 2; + if( vips_region_prepare( in, &rect ) ) + return( -1 ); + lsk = VIPS_REGION_LSKIP( in ); + psk = VIPS_IMAGE_SIZEOF_PEL( im ); + + for( angle = 0; angle < 256; angle++ ) { + /* For a white disc on a black bg, angles run anti-clockwise + * from 0 at the bottom. + */ + + } + + + /* For each of the four directions, the offset to get to that pixel + * from the top-left of our 3x3. offseta is the left/up direction, or + * the lower memory address. + */ + offseta[0] = psk; + offsetb[0] = psk + 2 * lsk; + offseta[1] = 2 * psk; + offsetb[1] = 2 * lsk; + offseta[2] = lsk; + offsetb[2] = 2 * psk + lsk; + offseta[3] = 0; + offsetb[3] = 2 * psk + 2 * lsk; + + for( y = 0; y < r->height; y++ ) { + VipsPel *p = (VipsPel * restrict) + VIPS_REGION_ADDR( in, r->left, r->top + y ); + VipsPel *q = (VipsPel * restrict) + VIPS_REGION_ADDR( or, r->left, r->top + y ); + + for( x = 0; x < r->width; x++ ) { + for( band = 0; band < out_bands; band++ ) { + int G = p[lsk + psk]; + int theta = p[lsk + psk + 1]; + VipsPel low = p[offseta[theta]]; + VipsPel high = p[offsetb[theta]]; + + /* Set G to 0 if it's not the local maxima in + * the direction of the gradient. If G is equal + * to the low side, also zero, so wide edges + * with equal gradient move (arbitarilly) left + * and up. + */ + if( G <= low || + G < high ) + G = 0; + + q[band] = G; + + p += 2; + } + + q += out_bands; + } + } + + return( 0 ); +} + +/* Remove non-maximal edges. At each point, compare the G to the G in either + * direction and 0 it if it's not the largest. + */ +static int +vips_canny_nonmax( VipsImage *in, VipsImage **out ) +{ + *out = vips_image_new(); + if( vips_image_pipelinev( *out, + VIPS_DEMAND_STYLE_THINSTRIP, in, NULL ) ) + return( -1 ); + (*out)->Bands /= 2; + (*out)->Xsize -= 2; + (*out)->Ysize -= 2; + + if( vips_image_generate( *out, + vips_start_one, vips_canny_nonmax_generate, vips_stop_one, + in, NULL ) ) + return( -1 ); + + return( 0 ); +} + +static int +vips_canny_build( VipsObject *object ) +{ + VipsCanny *canny = (VipsCanny *) object; + VipsImage **t = (VipsImage **) vips_object_local_array( object, 20 ); + + VipsImage *in; + VipsImage *Gx; + VipsImage *Gy; + + if( VIPS_OBJECT_CLASS( vips_canny_parent_class )->build( object ) ) + return( -1 ); + + in = canny->in; + + if( vips_gaussblur( in, &t[0], canny->sigma, NULL ) ) + return( -1 ); + in = t[0]; + + /* Separated Sobel gives Gx / Gy. Aim for vector path. Scale down by 2 + * to try to avoid clipping. + */ + t[1] = vips_image_new_matrixv( 1, 3, 1.0, 2.0, 1.0 ); + vips_image_set_double( t[1], "scale", 2.0 ); + t[2] = vips_image_new_matrixv( 3, 1, 1.0, 0.0, -1.0 ); + vips_image_set_double( t[2], "offset", 128.0 ); + if( vips_conv( in, &t[3], t[1], + "precision", VIPS_PRECISION_INTEGER, + NULL ) || + vips_conv( t[3], &t[4], t[2], + "precision", VIPS_PRECISION_INTEGER, + NULL ) ) + return( -1 ); + Gx = t[4]; + + t[5] = vips_image_new_matrixv( 3, 1, 1.0, 2.0, 1.0 ); + vips_image_set_double( t[5], "scale", 2.0 ); + t[6] = vips_image_new_matrixv( 1, 3, 1.0, 0.0, -1.0 ); + vips_image_set_double( t[6], "offset", 128.0 ); + if( vips_conv( in, &t[7], t[5], + "precision", VIPS_PRECISION_INTEGER, + NULL ) || + vips_conv( t[7], &t[8], t[6], + "precision", VIPS_PRECISION_INTEGER, + NULL ) ) + return( -1 ); + Gy = t[8]; + + /* Form (G, theta), with theta coded. + */ + canny->args[0] = Gx; + canny->args[1] = Gy; + canny->args[2] = NULL; + if( vips_canny_polar( canny->args, &t[9] ) ) + return( -1 ); + in = t[9]; + + /* Expand by two pixels all around, then thin. + if( vips_embed( in, &t[10], 1, 1, in->Xsize + 2, in->Ysize + 2, + "extend", VIPS_EXTEND_COPY, + NULL ) || + vips_canny_nonmax( t[10], &t[11] ) ) + return( -1 ); + in = t[11]; + */ + + g_object_set( object, "out", vips_image_new(), NULL ); + + if( vips_image_write( in, canny->out ) ) + return( -1 ); + + return( 0 ); +} + +static void +vips_canny_class_init( VipsCannyClass *class ) +{ + GObjectClass *gobject_class = G_OBJECT_CLASS( class ); + VipsObjectClass *object_class = (VipsObjectClass *) class; + VipsOperationClass *operation_class = VIPS_OPERATION_CLASS( class ); + + gobject_class->set_property = vips_object_set_property; + gobject_class->get_property = vips_object_get_property; + + object_class->nickname = "canny"; + object_class->description = _( "gaussian blur" ); + object_class->build = vips_canny_build; + + operation_class->flags = VIPS_OPERATION_SEQUENTIAL; + + VIPS_ARG_IMAGE( class, "in", 1, + _( "Input" ), + _( "Input image" ), + VIPS_ARGUMENT_REQUIRED_INPUT, + G_STRUCT_OFFSET( VipsCanny, in ) ); + + VIPS_ARG_IMAGE( class, "out", 2, + _( "Output" ), + _( "Output image" ), + VIPS_ARGUMENT_REQUIRED_OUTPUT, + G_STRUCT_OFFSET( VipsCanny, out ) ); + + VIPS_ARG_DOUBLE( class, "sigma", 10, + _( "Sigma" ), + _( "Sigma of Gaussian" ), + VIPS_ARGUMENT_OPTIONAL_INPUT, + G_STRUCT_OFFSET( VipsCanny, sigma ), + 0.01, 1000, 1.4 ); + +} + +static void +vips_canny_init( VipsCanny *canny ) +{ + canny->sigma = 1.4; +} + +/** + * vips_canny: (method) + * @in: input image + * @out: (out): output image + * @sigma: how large a mask to use + * @...: %NULL-terminated list of optional named arguments + * + * Optional arguments: + * + * * @sigma: %gdouble, sigma for gaussian blur + * + * See also: vips_gaussblur(). + * + * Returns: 0 on success, -1 on error. + */ +int +vips_canny( VipsImage *in, VipsImage **out, ... ) +{ + va_list ap; + int result; + + va_start( ap, out ); + result = vips_call_split( "canny", ap, in, out ); + va_end( ap ); + + return( result ); +} diff --git a/libvips/convolution/canny.c b/libvips/convolution/canny.c index b21aac22..f87cdcc3 100644 --- a/libvips/convolution/canny.c +++ b/libvips/convolution/canny.c @@ -88,10 +88,13 @@ vips_canny_polar_generate( VipsRegion *or, for( band = 0; band < Gx->Bands; band++ ) { int x = p1[band] - 128; int y = p2[band] - 128; - int a = VIPS_DEG( atan2( x, y ) ); + int a = 360 + VIPS_DEG( atan2( x, y ) ); - q[0] = sqrt( x * x + y * y ); - q[1] = ((a + 360 + 22) / 45) & 0x3; + /* Faster than hypot() for int args. Scale down + * or we'll clip on very hard edges. + */ + q[0] = 0.5 * sqrt( x * x + y * y ); + q[1] = ((a + 22) / 45) & 0x3; q += 2; } @@ -105,16 +108,16 @@ vips_canny_polar_generate( VipsRegion *or, } /* Calculate G/theta from Gx/Gy -- rather like rect -> polar, except that we - * constrain theta to 0/45/90/135, which we code as 0/1/2/3, and we don't care - * about the sign of G. + * code theta as below. Scale G down by 0.5 so that we + * don't clip on hard edges. * - * Theta values: + * For a white disc on a black background, theta is: * - * 3 | 0 | 1 - * --+---+-- - * 2 | x | 2 - * --+---+-- - * 1 | 0 | 3 + * 1 | 0 | 3 + * --+---+-- + * 2 | X | 2 + * --+---+-- + * 3 | 0 | 1 */ static int vips_canny_polar( VipsImage **args, VipsImage **out ) @@ -140,6 +143,7 @@ vips_canny_nonmax_generate( VipsRegion *or, VipsRegion *in = (VipsRegion *) vseq; VipsRect *r = &or->valid; VipsImage *im = in->im; + int out_bands = or->im->Bands; VipsRect rect; int x, y, band; @@ -157,16 +161,17 @@ vips_canny_nonmax_generate( VipsRegion *or, psk = VIPS_IMAGE_SIZEOF_PEL( im ); /* For each of the four directions, the offset to get to that pixel - * from the top-left of our 3x3. + * from the top-left of our 3x3. offseta is the left/up direction, or + * the lower memory address. */ offseta[0] = psk; offsetb[0] = psk + 2 * lsk; - offseta[1] = 2 * psk; - offsetb[1] = 2 * lsk; - offseta[2] = 2 * psk + lsk; - offsetb[2] = lsk; - offseta[3] = 0; - offsetb[3] = 2 * psk + 2 * lsk; + offseta[1] = 0; + offsetb[1] = 2 * psk + 2 * lsk; + offseta[2] = lsk; + offsetb[2] = 2 * psk + lsk; + offseta[3] = 2 * psk; + offsetb[3] = 2 * lsk; for( y = 0; y < r->height; y++ ) { VipsPel *p = (VipsPel * restrict) @@ -175,16 +180,19 @@ vips_canny_nonmax_generate( VipsRegion *or, VIPS_REGION_ADDR( or, r->left, r->top + y ); for( x = 0; x < r->width; x++ ) { - for( band = 0; band < im->Bands; band++ ) { + for( band = 0; band < out_bands; band++ ) { int G = p[lsk + psk]; int theta = p[lsk + psk + 1]; VipsPel low = p[offseta[theta]]; VipsPel high = p[offsetb[theta]]; /* Set G to 0 if it's not the local maxima in - * the direction of the gradient. + * the direction of the gradient. If G is equal + * to the low side, also zero, so wide edges + * with equal gradient move (arbitarilly) left + * and up. */ - if( G < low || + if( G <= low || G < high ) G = 0; @@ -193,7 +201,7 @@ vips_canny_nonmax_generate( VipsRegion *or, p += 2; } - q += im->Bands; + q += out_bands; } } @@ -241,9 +249,11 @@ vips_canny_build( VipsObject *object ) return( -1 ); in = t[0]; - /* Separated Sobel gives Gx / Gy. Aim for vector path. + /* Separated Sobel gives Gx / Gy. Aim for vector path. Scale down by 2 + * to try to avoid clipping. */ t[1] = vips_image_new_matrixv( 1, 3, 1.0, 2.0, 1.0 ); + vips_image_set_double( t[1], "scale", 2.0 ); t[2] = vips_image_new_matrixv( 3, 1, 1.0, 0.0, -1.0 ); vips_image_set_double( t[2], "offset", 128.0 ); if( vips_conv( in, &t[3], t[1], @@ -256,6 +266,7 @@ vips_canny_build( VipsObject *object ) Gx = t[4]; t[5] = vips_image_new_matrixv( 3, 1, 1.0, 2.0, 1.0 ); + vips_image_set_double( t[5], "scale", 2.0 ); t[6] = vips_image_new_matrixv( 1, 3, 1.0, 0.0, -1.0 ); vips_image_set_double( t[6], "offset", 128.0 ); if( vips_conv( in, &t[7], t[5], @@ -278,7 +289,7 @@ vips_canny_build( VipsObject *object ) /* Expand by two pixels all around, then thin. */ - if( vips_embed( in, &t[10], 1, 1, in->Xsize, in->Ysize, + if( vips_embed( in, &t[10], 1, 1, in->Xsize + 2, in->Ysize + 2, "extend", VIPS_EXTEND_COPY, NULL ) || vips_canny_nonmax( t[10], &t[11] ) )