425 lines
9.2 KiB
C
425 lines
9.2 KiB
C
/* @(#) Rank filter.
|
|
* @(#)
|
|
* @(#) int
|
|
* @(#) im_rank( in, out, xsize, ysize, order )
|
|
* @(#) IMAGE *in, *out;
|
|
* @(#) int xsize, ysize;
|
|
* @(#) int order;
|
|
* @(#)
|
|
* @(#) Also: im_rank_raw(). As above, but does not add a black border.
|
|
* @(#)
|
|
* @(#) Returns either 0 (success) or -1 (fail)
|
|
* @(#)
|
|
*
|
|
* Author: JC
|
|
* Written on: 19/8/96
|
|
* Modified on:
|
|
* JC 20/8/96
|
|
* - now uses insert-sort rather than bubble-sort
|
|
* - now works for any non-complex type
|
|
* JC 22/6/01
|
|
* - oops, sanity check on n wrong
|
|
* JC 28/8/03
|
|
* - cleanups
|
|
* - better selection algorithm ... same speed for 3x3, about 3x faster
|
|
* for 5x5, faster still for larger windows
|
|
* - index from zero for consistency with other parts of vips
|
|
* 7/4/04
|
|
* - now uses im_embed() with edge stretching on the input, not
|
|
* the output
|
|
* - sets Xoffset / Yoffset
|
|
* 7/10/04
|
|
* - oops, im_embed() size was wrong
|
|
*/
|
|
|
|
/*
|
|
|
|
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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
|
|
|
*/
|
|
|
|
/*
|
|
|
|
These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
|
|
|
|
*/
|
|
|
|
#ifdef HAVE_CONFIG_H
|
|
#include <config.h>
|
|
#endif /*HAVE_CONFIG_H*/
|
|
#include <vips/intl.h>
|
|
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <assert.h>
|
|
|
|
#include <vips/vips.h>
|
|
|
|
#ifdef WITH_DMALLOC
|
|
#include <dmalloc.h>
|
|
#endif /*WITH_DMALLOC*/
|
|
|
|
/* Global state: save our parameters here.
|
|
*/
|
|
typedef struct {
|
|
IMAGE *in, *out; /* Images we run */
|
|
int xsize, ysize; /* Window size */
|
|
int order; /* Element select */
|
|
int n; /* xsize * ysize */
|
|
} RankInfo;
|
|
|
|
/* Sequence value: just the array we sort in.
|
|
*/
|
|
typedef struct {
|
|
REGION *ir;
|
|
PEL *sort;
|
|
} SeqInfo;
|
|
|
|
/* Free a sequence value.
|
|
*/
|
|
static int
|
|
stop_rank( SeqInfo *seq, IMAGE *in, RankInfo *rnk )
|
|
{
|
|
/* Free attached objects.
|
|
*/
|
|
if( seq->ir ) {
|
|
im_region_free( seq->ir );
|
|
seq->ir = NULL;
|
|
}
|
|
|
|
return( 0 );
|
|
}
|
|
|
|
/* Rank start function.
|
|
*/
|
|
static void *
|
|
start_rank( IMAGE *out, IMAGE *in, RankInfo *rnk )
|
|
{
|
|
SeqInfo *seq = IM_NEW( out, SeqInfo );
|
|
|
|
if( !seq )
|
|
return( NULL );
|
|
|
|
/* Init!
|
|
*/
|
|
seq->ir = NULL;
|
|
seq->sort = NULL;
|
|
|
|
/* Attach region and arrays.
|
|
*/
|
|
seq->ir = im_region_create( in );
|
|
seq->sort = IM_ARRAY( out,
|
|
IM_IMAGE_SIZEOF_ELEMENT( in ) * rnk->n, PEL );
|
|
if( !seq->ir || !seq->sort ) {
|
|
stop_rank( seq, in, rnk );
|
|
return( NULL );
|
|
}
|
|
|
|
return( (void *) seq );
|
|
}
|
|
|
|
#define SWAP( TYPE, A, B ) { \
|
|
TYPE t = (A); \
|
|
(A) = (B); \
|
|
(B) = t; \
|
|
}
|
|
|
|
/* Inner loop for select-sorting TYPE.
|
|
*/
|
|
#define LOOP_SELECT( TYPE ) { \
|
|
TYPE *q = (TYPE *) IM_REGION_ADDR( or, le, y ); \
|
|
TYPE *p = (TYPE *) IM_REGION_ADDR( ir, le, y ); \
|
|
TYPE *sort = (TYPE *) seq->sort; \
|
|
TYPE a; \
|
|
\
|
|
for( x = 0; x < sz; x++ ) { \
|
|
TYPE *d = p + x; \
|
|
\
|
|
/* Copy window into sort[].
|
|
*/ \
|
|
for( k = 0, j = 0; j < rnk->ysize; j++ ) { \
|
|
for( i = 0; i < eaw; i += bands, k++ ) \
|
|
sort[k] = d[i]; \
|
|
d += ls; \
|
|
} \
|
|
\
|
|
/* Rearrange sort[] to make the order-th element the order-th
|
|
* smallest, adapted from Numerical Recipes in C.
|
|
*/ \
|
|
lower = 0; /* Range we know the result lies in */ \
|
|
upper = rnk->n - 1; \
|
|
for(;;) { \
|
|
if( upper - lower < 2 ) { \
|
|
/* 1 or 2 elements left.
|
|
*/ \
|
|
if( upper - lower == 1 && \
|
|
sort[lower] > sort[upper] ) \
|
|
SWAP( TYPE, \
|
|
sort[lower], sort[upper] ); \
|
|
break; \
|
|
} \
|
|
else { \
|
|
/* Pick mid-point of remaining elements.
|
|
*/ \
|
|
mid = (lower + upper) >> 1; \
|
|
\
|
|
/* Sort lower/mid/upper elements, hold
|
|
* midpoint in sort[lower + 1] for
|
|
* partitioning.
|
|
*/ \
|
|
SWAP( TYPE, sort[lower + 1], sort[mid] ); \
|
|
if( sort[lower] > sort[upper] ) \
|
|
SWAP( TYPE, \
|
|
sort[lower], sort[upper] ); \
|
|
if( sort[lower + 1] > sort[upper] ) \
|
|
SWAP( TYPE, \
|
|
sort[lower + 1], sort[upper] );\
|
|
if( sort[lower] > sort[lower + 1] ) \
|
|
SWAP( TYPE, \
|
|
sort[lower], sort[lower + 1] ) \
|
|
\
|
|
i = lower + 1; \
|
|
j = upper; \
|
|
a = sort[lower + 1]; \
|
|
\
|
|
for(;;) { \
|
|
/* Search for out of order elements.
|
|
*/ \
|
|
do \
|
|
i++; \
|
|
while( sort[i] < a ); \
|
|
do \
|
|
j--; \
|
|
while( sort[j] > a ); \
|
|
if( j < i ) \
|
|
break; \
|
|
SWAP( TYPE, sort[i], sort[j] ); \
|
|
} \
|
|
\
|
|
/* Replace mid element.
|
|
*/ \
|
|
sort[lower + 1] = sort[j]; \
|
|
sort[j] = a; \
|
|
\
|
|
/* Move to partition with the kth element.
|
|
*/ \
|
|
if( j >= rnk->order ) \
|
|
upper = j - 1; \
|
|
if( j <= rnk->order ) \
|
|
lower = i; \
|
|
} \
|
|
} \
|
|
\
|
|
q[x] = sort[rnk->order]; \
|
|
} \
|
|
}
|
|
|
|
/* Loop for find max of window.
|
|
*/
|
|
#define LOOP_MAX( TYPE ) { \
|
|
TYPE *q = (TYPE *) IM_REGION_ADDR( or, le, y ); \
|
|
TYPE *p = (TYPE *) IM_REGION_ADDR( ir, le, y ); \
|
|
\
|
|
for( x = 0; x < sz; x++ ) { \
|
|
TYPE *d = &p[x]; \
|
|
TYPE max; \
|
|
\
|
|
max = *d; \
|
|
for( j = 0; j < rnk->ysize; j++ ) { \
|
|
TYPE *e = d; \
|
|
\
|
|
for( i = 0; i < rnk->xsize; i++ ) { \
|
|
if( *e > max ) \
|
|
max = *e; \
|
|
\
|
|
e += bands; \
|
|
} \
|
|
\
|
|
d += ls; \
|
|
} \
|
|
\
|
|
q[x] = max; \
|
|
} \
|
|
}
|
|
|
|
/* Loop for find min of window.
|
|
*/
|
|
#define LOOP_MIN( TYPE ) { \
|
|
TYPE *q = (TYPE *) IM_REGION_ADDR( or, le, y ); \
|
|
TYPE *p = (TYPE *) IM_REGION_ADDR( ir, le, y ); \
|
|
\
|
|
for( x = 0; x < sz; x++ ) { \
|
|
TYPE *d = &p[x]; \
|
|
TYPE min; \
|
|
\
|
|
min = *d; \
|
|
for( j = 0; j < rnk->ysize; j++ ) { \
|
|
TYPE *e = d; \
|
|
\
|
|
for( i = 0; i < rnk->xsize; i++ ) { \
|
|
if( *e < min ) \
|
|
min = *e; \
|
|
\
|
|
e += bands; \
|
|
} \
|
|
\
|
|
d += ls; \
|
|
} \
|
|
\
|
|
q[x] = min; \
|
|
} \
|
|
}
|
|
|
|
#define SWITCH( OPERATION ) \
|
|
switch( rnk->out->BandFmt ) { \
|
|
case IM_BANDFMT_UCHAR: OPERATION( unsigned char ); break; \
|
|
case IM_BANDFMT_CHAR: OPERATION( signed char ); break; \
|
|
case IM_BANDFMT_USHORT: OPERATION( unsigned short ); break; \
|
|
case IM_BANDFMT_SHORT: OPERATION( signed short ); break; \
|
|
case IM_BANDFMT_UINT: OPERATION( unsigned int ); break; \
|
|
case IM_BANDFMT_INT: OPERATION( signed int ); break; \
|
|
case IM_BANDFMT_FLOAT: OPERATION( float ); break; \
|
|
case IM_BANDFMT_DOUBLE: OPERATION( double ); break; \
|
|
\
|
|
default: \
|
|
assert( 0 ); \
|
|
}
|
|
|
|
/* Rank of a REGION.
|
|
*/
|
|
static int
|
|
gen_rank( REGION *or, SeqInfo *seq, IMAGE *in, RankInfo *rnk )
|
|
{
|
|
REGION *ir = seq->ir;
|
|
|
|
Rect *r = &or->valid;
|
|
Rect s;
|
|
int le = r->left;
|
|
int to = r->top;
|
|
int bo = IM_RECT_BOTTOM(r);
|
|
int sz = IM_REGION_N_ELEMENTS( or );
|
|
|
|
int ls;
|
|
int bands = in->Bands;
|
|
int eaw = rnk->xsize * bands; /* elements across window */
|
|
|
|
int x, y;
|
|
int i, j, k;
|
|
int upper, lower, mid;
|
|
|
|
/* Prepare the section of the input image we need. A little larger
|
|
* than the section of the output image we are producing.
|
|
*/
|
|
s = *r;
|
|
s.width += rnk->xsize - 1;
|
|
s.height += rnk->ysize - 1;
|
|
if( im_prepare( ir, &s ) )
|
|
return( -1 );
|
|
ls = IM_REGION_LSKIP( ir ) / IM_IMAGE_SIZEOF_ELEMENT( in );
|
|
|
|
for( y = to; y < bo; y++ ) {
|
|
if( rnk->order == 0 )
|
|
SWITCH( LOOP_MIN )
|
|
else if( rnk->order == rnk->n - 1 )
|
|
SWITCH( LOOP_MAX )
|
|
else
|
|
SWITCH( LOOP_SELECT ) }
|
|
|
|
return( 0 );
|
|
}
|
|
|
|
/* Rank filter.
|
|
*/
|
|
int
|
|
im_rank_raw( IMAGE *in, IMAGE *out, int xsize, int ysize, int order )
|
|
{
|
|
RankInfo *rnk;
|
|
|
|
/* Check parameters.
|
|
*/
|
|
if( !in || in->Coding != IM_CODING_NONE || im_iscomplex( in ) ) {
|
|
im_errormsg( "im_rank: input non-complex uncoded only" );
|
|
return( -1 );
|
|
}
|
|
if( xsize > 1000 || ysize > 1000 || xsize <= 0 || ysize <= 0 ||
|
|
order < 0 || order > xsize * ysize - 1 ) {
|
|
im_errormsg( "im_rank: bad parameters" );
|
|
return( -1 );
|
|
}
|
|
if( im_piocheck( in, out ) )
|
|
return( -1 );
|
|
|
|
/* Save parameters.
|
|
*/
|
|
if( !(rnk = IM_NEW( out, RankInfo )) )
|
|
return( -1 );
|
|
rnk->in = in;
|
|
rnk->out = out;
|
|
rnk->xsize = xsize;
|
|
rnk->ysize = ysize;
|
|
rnk->order = order;
|
|
rnk->n = xsize * ysize;
|
|
|
|
/* Prepare output. Consider a 7x7 window and a 7x7 image --- the output
|
|
* would be 1x1.
|
|
*/
|
|
if( im_cp_desc( out, in ) )
|
|
return( -1 );
|
|
out->Xsize -= xsize - 1;
|
|
out->Ysize -= ysize - 1;
|
|
if( out->Xsize <= 0 || out->Ysize <= 0 ) {
|
|
im_errormsg( "im_rank: image too small for window" );
|
|
return( -1 );
|
|
}
|
|
|
|
/* Set demand hints. FATSTRIP is good for us, as THINSTRIP will cause
|
|
* too many recalculations on overlaps.
|
|
*/
|
|
if( im_demand_hint( out, IM_FATSTRIP, in, NULL ) )
|
|
return( -1 );
|
|
|
|
/* Generate!
|
|
*/
|
|
if( im_generate( out, start_rank, gen_rank, stop_rank, in, rnk ) )
|
|
return( -1 );
|
|
|
|
out->Xoffset = -xsize / 2;
|
|
out->Yoffset = -ysize / 2;
|
|
|
|
return( 0 );
|
|
}
|
|
|
|
/* The above, with a border to make out the same size as in.
|
|
*/
|
|
int
|
|
im_rank( IMAGE *in, IMAGE *out, int xsize, int ysize, int order )
|
|
{
|
|
IMAGE *t1 = im_open_local( out, "im_rank:1", "p" );
|
|
|
|
if( !t1 ||
|
|
im_embed( in, t1, 1,
|
|
xsize/2, ysize/2,
|
|
in->Xsize + xsize - 1, in->Ysize + ysize - 1 ) ||
|
|
im_rank_raw( t1, out, xsize, ysize, order ) )
|
|
return( -1 );
|
|
|
|
out->Xoffset = 0;
|
|
out->Yoffset = 0;
|
|
|
|
return( 0 );
|
|
}
|