665 lines
16 KiB
C
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*/
|
||
|
|
}
|