284 lines
6.6 KiB
C
284 lines
6.6 KiB
C
/* generate gaussian masks
|
|
*
|
|
* Written on: 30/11/1989 by Nicos
|
|
* Updated on: 6/12/1991
|
|
* 7/8/96 JC
|
|
* - ansified, mem leaks plugged
|
|
* 20/11/98 JC
|
|
* - mask too large check added
|
|
* 18/3/09
|
|
* - bumped max mask size *40
|
|
* - added _sep variant
|
|
* 30/3/09
|
|
* - set scale in _sep variant, why not
|
|
* 21/10/10
|
|
* - gtkdoc
|
|
*/
|
|
|
|
/*
|
|
|
|
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 <vips/vips.h>
|
|
|
|
#ifdef WITH_DMALLOC
|
|
#include <dmalloc.h>
|
|
#endif /*WITH_DMALLOC*/
|
|
|
|
#define IM_MAXMASK 5000
|
|
|
|
/**
|
|
* im_gauss_dmask:
|
|
* @filename: the returned mask has this set as the filename
|
|
* @sigma: standard deviation of mask
|
|
* @min_ampl: minimum amplitude
|
|
*
|
|
* im_gauss_dmask() creates a circularly symmetric Gaussian mask of radius
|
|
* @sigma. The size of the mask is determined by the variable @min_ampl;
|
|
* if for instance the value .1 is entered this means that the produced mask
|
|
* is clipped at values less than 10 percent of the maximum amplitude.
|
|
*
|
|
* The program uses the following equation:
|
|
*
|
|
* H(r) = exp( -(r * r) / (2 * sigma * sigma) )
|
|
*
|
|
* The generated mask has odd size and its maximum value is normalised to
|
|
* 1.0.
|
|
*
|
|
* See also: im_gauss_imask(), im_gauss_imask_sep(), im_log_dmask(), im_conv().
|
|
*
|
|
* Returns: the calculated mask on success, or NULL on error.
|
|
*/
|
|
DOUBLEMASK *
|
|
im_gauss_dmask( const char *filename, double sigma, double min_ampl )
|
|
{
|
|
int x, y, k;
|
|
double distance;
|
|
double temp;
|
|
double *pt1, *pt2, *pt3, *pt4;
|
|
int max_x;
|
|
int xm, ym;
|
|
int xm2, ym2; /* xm2 = xm/2 */
|
|
int offset;
|
|
double *cf, *cfs, *mc;
|
|
DOUBLEMASK *m;
|
|
double sig2, sum; /* sig2 = 2. * sigma * sigma */
|
|
|
|
/* Find the size of the mask depending on the entered data
|
|
*/
|
|
sig2 = 2. * sigma * sigma;
|
|
max_x = 8 * sigma > IM_MAXMASK ? IM_MAXMASK : 8 * sigma ;
|
|
for( x = 0; x < max_x; x++ ) {
|
|
temp = exp( - ((double)(x * x))/sig2 );
|
|
if( temp < min_ampl )
|
|
break;
|
|
}
|
|
if( x == max_x ) {
|
|
im_error( "im_gauss_dmask", "%s", _( "mask too large" ) );
|
|
return( NULL );
|
|
}
|
|
|
|
xm2 = x;
|
|
ym2 = x;
|
|
xm = xm2 * 2 + 1;
|
|
ym = ym2 * 2 + 1;
|
|
|
|
if( !(cfs = IM_ARRAY( NULL, (xm2+1)*(ym2+1), double )) )
|
|
return( NULL );
|
|
|
|
for( k = 0, y = 0; y <= ym2; y++ ) {
|
|
for( x = 0; x <= xm2; x++, k++ ) {
|
|
distance = x*x + y*y;
|
|
cfs[k] = exp( -distance / sig2 );
|
|
}
|
|
}
|
|
|
|
#ifdef PIM_RINT
|
|
for( k = 0, y = 0; y <= ymask_2; y++ ) {
|
|
for( x = 0; x <= xmask_2; x++, k++ )
|
|
fprintf(stderr, "%3.2f ", cfs[k] );
|
|
fprintf(stderr, "\n");
|
|
}
|
|
#endif
|
|
|
|
if( !(m = im_create_dmask( filename, xm, ym )) ) {
|
|
im_free( cfs );
|
|
return( NULL );
|
|
}
|
|
|
|
/* copy the 1/4 cfs into the m
|
|
*/
|
|
cf = cfs;
|
|
offset = xm2 * (xm + 1);
|
|
mc = m->coeff + offset;
|
|
for( y = 0; y <= ym2; y++ ) {
|
|
for( x = 0; x <= xm2; x++ ) {
|
|
pt1 = mc + (y * xm) + x;
|
|
pt2 = mc - (y * xm) + x;
|
|
pt3 = mc + (y * xm) - x;
|
|
pt4 = mc - (y * xm) - x;
|
|
|
|
*pt1 = cf[x];
|
|
*pt2 = cf[x];
|
|
*pt3 = cf[x];
|
|
*pt4 = cf[x];
|
|
}
|
|
|
|
cf += (xm2 + 1);
|
|
}
|
|
im_free( cfs );
|
|
|
|
sum = 0.0;
|
|
for( k = 0, y = 0; y < m->ysize; y++ )
|
|
for( x = 0; x < m->xsize; x++, k++ )
|
|
sum += m->coeff[k];
|
|
m->scale = sum;
|
|
m->offset = 0.0;
|
|
|
|
#ifdef PIM_RINT
|
|
im_print_dmask( m );
|
|
#endif
|
|
return( m );
|
|
}
|
|
|
|
/**
|
|
* im_gauss_dmask_sep:
|
|
* @filename: the returned mask has this set as the filename
|
|
* @sigma: standard deviation of mask
|
|
* @min_ampl: minimum amplitude
|
|
*
|
|
* im_gauss_dmask_sep() works exactly as im_gauss_dmask(), but returns only
|
|
* the central line of the mask. This is useful with im_convsepf().
|
|
*
|
|
* See also: im_gauss_dmask(), im_convsepf().
|
|
*
|
|
* Returns: the calculated mask on success, or NULL on error.
|
|
*/
|
|
DOUBLEMASK *
|
|
im_gauss_dmask_sep( const char *filename, double sigma, double min_ampl )
|
|
{
|
|
DOUBLEMASK *im;
|
|
DOUBLEMASK *im2;
|
|
int i;
|
|
double sum;
|
|
|
|
if( !(im = im_gauss_dmask( filename, sigma, min_ampl )) )
|
|
return( NULL );
|
|
if( !(im2 = im_create_dmask( filename, im->xsize, 1 )) ) {
|
|
im_free_dmask( im );
|
|
return( NULL );
|
|
}
|
|
|
|
sum = 0;
|
|
for( i = 0; i < im->xsize; i++ ) {
|
|
im2->coeff[i] = im->coeff[i + im->xsize * (im->ysize / 2)];
|
|
sum += im2->coeff[i];
|
|
}
|
|
im2->scale = sum;
|
|
|
|
im_free_dmask( im );
|
|
|
|
return( im2 );
|
|
}
|
|
|
|
/**
|
|
* im_gauss_imask:
|
|
* @filename: the returned mask has this set as the filename
|
|
* @sigma: standard deviation of mask
|
|
* @min_ampl: minimum amplitude
|
|
*
|
|
* im_gauss_imask() works exactly as im_gauss_dmask(), but the returned mask
|
|
* is scaled so that it's maximum value it set to 100.
|
|
*
|
|
* See also: im_gauss_dmask(), im_gauss_imask_sep(), im_conv(), im_convsep().
|
|
*
|
|
* Returns: the calculated mask on success, or NULL on error.
|
|
*/
|
|
INTMASK *
|
|
im_gauss_imask( const char *filename, double sigma, double min_ampl )
|
|
{
|
|
DOUBLEMASK *dm;
|
|
INTMASK *im;
|
|
|
|
if( !(dm = im_gauss_dmask( filename, sigma, min_ampl )) )
|
|
return( NULL );
|
|
|
|
if( !(im = im_scale_dmask( dm, dm->filename )) ) {
|
|
im_free_dmask( dm );
|
|
return( NULL );
|
|
}
|
|
im_free_dmask( dm );
|
|
|
|
return( im ) ;
|
|
}
|
|
|
|
/**
|
|
* im_gauss_imask_sep:
|
|
* @filename: the returned mask has this set as the filename
|
|
* @sigma: standard deviation of mask
|
|
* @min_ampl: minimum amplitude
|
|
*
|
|
* im_gauss_imask_sep() works exactly as im_gauss_imask(), but returns only
|
|
* the central line of the mask. This is useful with im_convsep().
|
|
*
|
|
* See also: im_gauss_dmask(), im_convsep().
|
|
*
|
|
* Returns: the calculated mask on success, or NULL on error.
|
|
*/
|
|
INTMASK *
|
|
im_gauss_imask_sep( const char *filename, double sigma, double min_ampl )
|
|
{
|
|
INTMASK *im;
|
|
INTMASK *im2;
|
|
int i;
|
|
int sum;
|
|
|
|
if( !(im = im_gauss_imask( filename, sigma, min_ampl )) )
|
|
return( NULL );
|
|
if( !(im2 = im_create_imask( filename, im->xsize, 1 )) ) {
|
|
im_free_imask( im );
|
|
return( NULL );
|
|
}
|
|
|
|
sum = 0;
|
|
for( i = 0; i < im->xsize; i++ ) {
|
|
im2->coeff[i] = im->coeff[i + im->xsize * (im->ysize / 2)];
|
|
sum += im2->coeff[i];
|
|
}
|
|
im2->scale = sum;
|
|
|
|
im_free_imask( im );
|
|
|
|
return( im2 );
|
|
}
|