libvips/libsrc/freq_filt/fmaskcir.c
2007-08-29 16:23:50 +00:00

665 lines
16 KiB
C

/* @(#) Typical filter function
* @(#) va_list is flag, filter parameters
* @(#)
* @(#) The following masks are implemented in this file
* @(#) flag, filter shape, parameters
* @(#) band pass ring reject filters
* @(#) 12 -\> idealbpf, parameters: frequency cutoff, width
* @(#) 13 -\> idealbrf, parameters: frequency cutoff, width
* @(#) 14 -\> butbpf, parameters: order, freq cutoff, width, ampl cutoff
* @(#) 15 -\> butbrf, parameters: order, freq cutoff, width, ampl cutoff
* @(#) 16 -\> gaussianbpf, parameters: frequency cutoff, width, ampl cutoff
* @(#) 17 -\> gaussianbrf, parameters: frequency cutoff, width, ampl cutoff
* @(#)
* @(#) The whole mask is created at once and written into the image file
* @(#)
* @(#) The following functions are contained within this file:
* @(#) Details are preceding the source code of each function
* @(#)
* @(#) int im__fmaskcir( out, flag, ap)
* @(#) IMAGE *out;
* @(#) enum mask_type flag;
* @(#) va_list ap;
* @(#)
*
* Copyright: N. Dessipris, 1991
* Written on: Nov 1991
* Updated on: Dec 1991
* 20/9/95 JC
* - modernised
*/
/*
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 <math.h>
#include <stdarg.h>
#include <vips/vips.h>
#include <vips/fmask.h>
#ifdef WITH_DMALLOC
#include <dmalloc.h>
#endif /*WITH_DMALLOC*/
/************************************************************************
* malloc space and create normalised coefficients accross
* the x (horizontal) and y (vertical) direction.
* xs, ys are the image sizes
* xd and yd are the scrambled distributions of x and y in the rotated
* Fourier transform
* xplusd is the non scrambled distribution of (x+x0)*(x+x0) centred at 0
* xminus is the non scrambled distribution of (x-x0)*(x-x0) centred at 0
* similar for yplusd and yminusd
************************************************************************/
static int
alloc( IMAGE *out,
int xs, int ys,
int **xd, int **yd,
int **xplusd, int **xminusd, int **yplusd, int **yminusd,
int x0, int y0,
float **line )
{
int i;
int *x, *y, *xp, *xm, *yp, *ym;
int *pp, *pm;
float *l;
x = IM_ARRAY( out, xs, int );
y = IM_ARRAY( out, ys, int );
xp = IM_ARRAY( out, xs, int );
xm = IM_ARRAY( out, xs, int );
yp = IM_ARRAY( out, ys, int );
ym = IM_ARRAY( out, ys, int );
l = IM_ARRAY( out, xs, float );
if( !x || !y || !xp || !xm || !yp || !ym || !l )
return( -1 );
/* if ys = 8 then y = {0,1,2,3,-4,-3,-2,-1}.
*/
for( i = 0; i < ys/2; i++ ) {
y[i] = i;
y[i+ys/2] = -ys/2 + i;
}
for( i = 0; i < xs/2; i++ ) {
x[i] = i;
x[i+xs/2] = -xs/2 + i;
}
*xd = x;
*yd = y;
pp = yp + ys/2;
pm = ym + ys/2;
for( i = -ys/2; i < ys/2; i++ ) {
pp[i] = (i + y0)*(i + y0);
pm[i] = (i - y0)*(i - y0);
}
*yplusd = yp + ys/2;
*yminusd = ym + ys/2;
pp = xp + xs/2;
pm = xm + xs/2;
for( i = -xs/2; i < xs/2; i++ ) {
pp[i] = (i+x0)*(i+x0);
pm[i] = (i-x0)*(i-x0);
}
*xplusd = xp + xs/2;
*xminusd = xm + xs/2;
*line = l;
return( 0 );
}
/************************************************************************/
/* FLAG = 12 */
/* Creates an ideal band pass filter mask */
/* The band is two CIRCLEs of radius r centred */
/* at (fcx, fcy) and (-fcx, -fcy) */
/************************************************************************/
static int
ideal_bpf( IMAGE *out, double fcx, double fcy, double r )
{
int x, y;
int xs = out->Xsize;
int ys = out->Ysize;
float *line, *cpline;
int *xd, *yd, *xplusx0d, *xminusx0d, *yplusy0d, *yminusy0d;
int x0, y0, d1_2, d2_2, r2;
int y2plus, y2minus;
if( xs != ys ) {
im_errormsg( "ideal_bpf: bad sizes" );
return( -1 );
}
if( fabs(fcx) <= 1.0 && fabs(fcy) < 1.0 && r > 0.0 && r < 1.0 ) {
x0 = fcx*xs / 2.0;
y0 = fcy*ys / 2.0;
r2 = r*r*xs / 4.0;
}
else if( fabs(fcx) < xs/2 && fabs(fcy) < ys/2 && r >= 1.0 ) {
x0 = fcx;
y0 = fcy;
r2 = r*r;
}
else {
im_errormsg( "ideal_bpf: bad args" );
return( -1 );
}
if( alloc( out, xs, ys,
&xd, &yd, &xplusx0d, &xminusx0d,
&yplusy0d, &yminusy0d, x0, y0, &line ) )
return( -1 );
for( y = 0; y < ys; y++ ) {
cpline = line;
y2plus = yplusy0d[yd[y]];
y2minus = yminusy0d[yd[y]];
for( x = 0; x < xs; x++ ) {
d1_2 = xminusx0d[xd[x]] + y2minus;
d2_2 = xplusx0d[xd[x]] + y2plus;
if( d1_2 <= r2 )
*cpline = 1.0;
else if( d2_2 <= r2 )
*cpline = 1.0;
else
*cpline = 0.0;
if( x == 0 && y == 0 )
*cpline = 1.0; /* allow the dc component */
cpline++;
}
if( im_writeline( y, out, (PEL *) line ) )
return( -1 );
}
return( 0 );
}
/************************************************************************/
/* FLAG = 13 */
/* Creates an ideal band reject filter mask */
/* The band is a CIRCLE of internal radius fc-df and external radius fc+df*/
/************************************************************************/
static int
ideal_brf( IMAGE *out, double fcx, double fcy, double r )
{
int x, y;
int xs = out->Xsize;
int ys = out->Ysize;
float *line, *cpline;
int *xd, *yd, *xplusx0d, *xminusx0d, *yplusy0d, *yminusy0d;
int x0, y0, d1_2, d2_2, r2;
int y2plus, y2minus;
if( xs != ys ) {
im_errormsg( "ideal_brf: bad args" );
return( -1 );
}
if( fabs(fcx) <= 1.0 && fabs(fcy) <= 1.0 && r > 0.0 && r < 1.0 ) {
x0 = fcx*xs / 2.0;
y0 = fcy*ys / 2.0;
r2 = r*r*xs / 4.0;
}
else if( fabs(fcx) < xs/2 && fabs(fcy) < ys/2 && r >= 1.0 ) {
x0 = fcx;
y0 = fcy;
r2 = r*r;
}
else {
im_errormsg( "ideal_brf: bad args" );
return( -1 );
}
if( alloc( out, xs, ys,
&xd, &yd, &xplusx0d, &xminusx0d,
&yplusy0d, &yminusy0d, x0, y0, &line ) )
return( -1 );
for( y = 0; y < ys; y++ ) {
cpline = line;
y2plus = yplusy0d[yd[y]];
y2minus = yminusy0d[yd[y]];
for( x = 0; x < xs; x++ ) {
d1_2 = xminusx0d[xd[x]] + y2minus;
d2_2 = xplusx0d[xd[x]] + y2plus;
if( d1_2 <= r2 )
*cpline = 0.0;
else if( d2_2 <= r2 )
*cpline = 0.0;
else
*cpline = 1.0;
if( x == 0 && y == 0 )
*cpline = 1.0;
cpline++;
}
if( im_writeline( y, out, (PEL *) line ) )
return( -1 );
}
return( 0 );
}
/************************************************************************/
/* FLAG = 14 */
/* Creates a butterworth band pass filter mask */
/* The band is two CIRCLES centred at (fcx, fcy) and (-fcx, -fcy) */
/* The program assummes that the peaks of the 2d mask are at the */
/* centres above and are set to 1.0. The amplitude of both circle mask */
/* are added and the cuttof frequency is calculated on the plane */
/* which passes though the centre of the circles and the 0 point */
/************************************************************************/
static int
butterworth_bpf( IMAGE *out,
double order, double fcx, double fcy, double r, double ac )
{
int x, y;
int xs = out->Xsize;
int ys = out->Ysize;
float *line, *cpline;
int *xd, *yd, *xplusx0d, *xminusx0d, *yplusy0d, *yminusy0d;
int x0, y0;
double cnst, cnsta, d1_2, d2_2, nr2; /* nr2 is new r squared */
int y2plus, y2minus;
if( xs != ys || order < 1.0 ) {
im_errormsg( "butterworth_bpf: bad sizes" );
return( -1 );
}
if( fabs(fcx) <= 1.0 && fabs(fcy) <= 1.0 && r > 0.0 && r < 1.0 ) {
x0 = fcx*xs / 2.0;
y0 = fcy*ys / 2.0;
nr2 = r*r*xs*xs / 4.0;
}
else if( fabs(fcx) < xs/2 && fabs(fcy) < ys/2 && r >= 1.0 ) {
x0 = fcx;
y0 = fcy;
nr2 = r*r;
}
else {
im_errormsg( "butterworth_bpf: bad args" );
return( -1 );
}
if( ac >= 1.0 || ac < 0.0) {
im_errormsg( "butterworth_bpf: bad args" );
return( -1 );
}
if( alloc( out, xs, ys,
&xd, &yd, &xplusx0d, &xminusx0d,
&yplusy0d, &yminusy0d, x0, y0, &line ) )
return( -1 );
/* Filter shape: radius d0, centres at (x0, y0), (-x0,-y0)
* H(d) = H1(d) + H2(d)
* H(d) = cnst1/(1 + cnst2 * pow((d-d0)/d0, 2*order)) +
* cnst1/(1 + cnst2 * pow((d+d0)/d0, 2*order));
* for d=+d0 H(+d0) = 1.0; for d=-d0 H(-d0) = 1.0;
* for d=+da H(+da) = ampl_cutof; for d=-da H1(-da) = ampl_cutof;
* da = (xa, ya)
* xa = x0*(1 - radius/sqrt(x0*x0+y0*y0))
* ya = y0*(1 - radius/sqrt(x0*x0+y0*y0))
*/
cnst = (1.0/ac) - 1.0;
/* normalise the amplitude at (x0,y0) to 1.0
*/
cnsta = 1.0 / (1.0 + 1.0 /
(1.0 + cnst*pow( 4.0*(x0*x0 + y0*y0)/nr2, order )));
for( y = 0; y < ys; y++ ) {
cpline = line;
y2plus = yplusy0d[yd[y]];
y2minus = yminusy0d[yd[y]];
for( x = 0; x < xs; x++ ) {
d1_2 = xminusx0d[xd[x]] + y2minus;
d2_2 = xplusx0d[xd[x]] + y2plus;
*cpline = cnsta * (
1.0 / (1.0 + cnst * pow( d1_2/nr2, order )) +
1.0 / (1.0 + cnst * pow( d2_2/nr2, order )) );
if( x == 0 && y == 0 )
*cpline = 1.0;
cpline++;
}
if( im_writeline( y, out, (PEL *) line ) )
return( -1 );
}
return( 0 );
}
/************************************************************************/
/* FLAG = 15 */
/* Creates a butterworth band pass filter mask */
/* The band is a the 1-H(f) of above */
/************************************************************************/
static int
butterworth_brf( IMAGE *out,
double order, double fcx, double fcy, double r, double ac )
{
int x, y;
int xs = out->Xsize;
int ys = out->Ysize;
float *line, *cpline;
int *xd, *yd, *xplusx0d, *xminusx0d, *yplusy0d, *yminusy0d;
int x0, y0;
double cnst, cnsta, d1_2, d2_2, nr2; /* nr2 is new r squared */
int y2plus, y2minus;
if( xs != ys || order < 1.0 ) {
im_errormsg( "butterworth_brf: bad sizes" );
return( -1 );
}
if( fabs(fcx) <= 1.0 && fabs(fcy) <= 1.0 && r > 0.0 && r < 1.0 ) {
x0 = fcx * xs / 2.0;
y0 = fcy * ys / 2.0;
nr2 = r*r*xs*xs / 4.0;
}
else if( fabs(fcx) < xs/2 && fabs(fcy) < ys/2 && r >= 1.0 ) {
x0 = fcx;
y0 = fcy;
nr2 = r*r;
}
else {
im_errormsg( "butterworth_brf: bad args" );
return( -1 );
}
if( ac >= 1.0 || ac < 0.0) {
im_errormsg( "butterworth_brf: bad args" );
return( -1 );
}
if( alloc( out, xs, ys,
&xd, &yd, &xplusx0d, &xminusx0d,
&yplusy0d, &yminusy0d, x0, y0, &line ) )
return( -1 );
cnst = (1.0/ac) - 1.0;
/* normalise the amplitude at (x0,y0) to 1.0
*/
cnsta = 1.0 / (1.0 + 1.0 / (1.0 +
cnst * pow( 4.0*(x0*x0 + y0*y0)/nr2, order )));
for( y = 0; y < ys; y++ ) {
cpline = line;
y2plus = yplusy0d[yd[y]];
y2minus = yminusy0d[yd[y]];
for( x = 0; x < xs; x++ ) {
d1_2 = xminusx0d[xd[x]] + y2minus;
d2_2 = xplusx0d[xd[x]] + y2plus;
if( d1_2 == 0.0 || d2_2 == 0.0 )
*cpline = 0;
else
*cpline = 1.0 - cnsta *
( 1.0/(1.0 + cnst*pow( d1_2/nr2, order )) +
1.0/(1.0 + cnst*pow( d2_2/nr2, order )) );
if( x == 0 && y == 0 )
*cpline = 1.0;
cpline++;
}
if( im_writeline( y, out, (PEL *) line ) )
return( -1 );
}
return( 0 );
}
/************************************************************************/
/* FLAG = 16 */
/* Creates a gaussian band pass filter mask */
/* The band is a RING of internal radius fc-df and external radius fc+df*/
/************************************************************************/
static int
gaussian_bpf( IMAGE *out, double fcx, double fcy, double r, double ac )
{
int x, y;
int xs = out->Xsize;
int ys = out->Ysize;
float *line, *cpline;
int *xd, *yd, *xplusx0d, *xminusx0d, *yplusy0d, *yminusy0d;
int x0, y0;
double cnst, cnsta, d1_2, d2_2, nr2; /* nr2 is new r squared */
int y2plus, y2minus;
if( xs != ys ) {
im_errormsg( "gauss_bpf: bad sizes" );
return( -1 );
}
if( fabs(fcx) <= 1.0 && fabs(fcy) <= 1.0 && r > 0.0 && r < 1.0 ) {
x0 = fcx*xs / 2.0;
y0 = fcy*ys / 2.0;
nr2 = r*r*xs*xs / 4.0;
}
else if( fabs(fcx) < xs/2 && fabs(fcy) < ys/2 && r >= 1.0 ) {
x0 = fcx;
y0 = fcy;
nr2 = r*r;
}
else {
im_errormsg( "gauss_bpf: bad args (f)" );
return( -1 );
}
if( ac >= 1.0 || ac < 0.0 ) {
im_errormsg( "gauss_bpf: bad args (ac)" );
return( -1 );
}
if( alloc( out, xs, ys,
&xd, &yd, &xplusx0d, &xminusx0d,
&yplusy0d, &yminusy0d, x0, y0, &line ) )
return( -1 );
cnst = -log( ac );
/* normalise the amplitude at (x0,y0) to 1.0
*/
cnsta = 1.0/(1.0 + exp( - cnst * 4.0 * (x0*x0+y0*y0) / nr2 ));
for( y = 0; y < ys; y++ ) {
cpline = line;
y2plus = yplusy0d[yd[y]];
y2minus = yminusy0d[yd[y]];
for( x = 0; x < xs; x++ ) {
d1_2 = xminusx0d[xd[x]] + y2minus;
d2_2 = xplusx0d[xd[x]] + y2plus;
*cpline = cnsta *
(exp( -cnst * d1_2/nr2 ) +
exp( -cnst * d2_2/nr2 ));
if( x == 0 && y == 0 )
*cpline = 1.0;
cpline++;
}
if( im_writeline( y, out, (PEL *) line ) )
return( -1 );
}
return( 0 );
}
/************************************************************************/
/* FLAG = 17 */
/* Creates a gaussian band reject filter mask */
/* The band is a RING of internal radius fc-df and external radius fc+df*/
/************************************************************************/
static int
gaussian_brf( IMAGE *out, double fcx, double fcy, double r, double ac )
{
int x, y;
int xs = out->Xsize;
int ys = out->Ysize;
float *line, *cpline;
int *xd, *yd, *xplusx0d, *xminusx0d, *yplusy0d, *yminusy0d;
int x0, y0;
double cnst, cnsta, d1_2, d2_2, nr2; /* nr2 is new r squared */
int y2plus, y2minus;
if( xs != ys ) {
im_errormsg( "gauss_brf: bad sizes" );
return( -1 );
}
if( fabs(fcx) <= 1.0 && fabs(fcy) <= 1.0 && r > 0.0 && r < 1.0 ) {
x0 = fcx*xs / 2.0;
y0 = fcy*ys / 2.0;
nr2 = r*r*xs*xs / 4.0;
}
else if( fabs(fcx) < xs/2 && fabs(fcy) < ys/2 && r >= 1.0 ) {
x0 = fcx;
y0 = fcy;
nr2 = r * r;
}
else {
im_errormsg( "gauss_brf: bad args" );
return( -1 );
}
if( ac >= 1.0 || ac < 0.0 ) {
im_errormsg( "gauss_brf: bad args" );
return( -1 );
}
if( alloc( out, xs, ys,
&xd, &yd, &xplusx0d, &xminusx0d,
&yplusy0d, &yminusy0d, x0, y0, &line ) )
return( -1 );
cnst = -log( ac );
/* normalise the amplitude at (x0,y0) to 1.0
*/
cnsta = 1.0/(1.0 + exp( - cnst * 4.0 * (x0*x0+y0*y0) / nr2 ));
for( y = 0; y < ys; y++ ) {
cpline = line;
y2plus = yplusy0d[yd[y]];
y2minus = yminusy0d[yd[y]];
for( x = 0; x < xs; x++ ) {
d1_2 = xminusx0d[xd[x]] + y2minus;
d2_2 = xplusx0d[xd[x]] + y2plus;
*cpline = 1.0 - cnsta *
(exp( -cnst * d1_2/nr2 ) +
exp( -cnst * d2_2/nr2 ));
if( x == 0 && y == 0 )
*cpline = 1.0;
cpline++;
}
if( im_writeline( y, out, (PEL *) line ) )
return( -1 );
}
return( 0 );
}
/* Creates bandpass filter masks
* All arithmetic internally has been carried out in double precision;
* however all masks are written as floats with maximum value normalised to 1.0
*/
int
im__fmaskcir( IMAGE *out, MaskType flag, va_list ap )
{
/* May be fewer than 5 args ... but extract them all anyway. Should be
* safe.
*/
double p0 = va_arg( ap, double );
double p1 = va_arg( ap, double );
double p2 = va_arg( ap, double );
double p3 = va_arg( ap, double );
double p4 = va_arg( ap, double );
switch( flag ) {
/* Band pass - band reject.
*/
case MASK_IDEAL_BANDPASS:
return( ideal_bpf( out, p0, p1, p2 ) );
case MASK_IDEAL_BANDREJECT:
return( ideal_brf( out, p0, p1, p2 ) );
case MASK_BUTTERWORTH_BANDPASS:
return( butterworth_bpf( out, p0, p1, p2, p3, p4 ) );
case MASK_BUTTERWORTH_BANDREJECT:
return( butterworth_brf( out, p0, p1, p2, p3, p4 ) );
case MASK_GAUSS_BANDPASS:
return( gaussian_bpf( out, p0, p1, p2, p3 ) );
case MASK_GAUSS_BANDREJECT:
return( gaussian_brf( out, p0, p1, p2, p3 ) );
default:
im_errormsg( "im__fmaskcir: unimplemented mask" );
return( -1 );
}
/*NOTREACHED*/
}