00001
00002 /***************************************************************************
00003  *  laplace.cpp - Implementation of a laplace filter
00004  *
00005  *  Created: Thu Jun 16 16:30:23 2005
00006  *  Copyright  2005-2007  Tim Niemueller [www.niemueller.de]
00007  *
00008  ****************************************************************************/
00009
00010 /*  This program is free software; you can redistribute it and/or modify
00011  *  it under the terms of the GNU General Public License as published by
00012  *  the Free Software Foundation; either version 2 of the License, or
00013  *  (at your option) any later version. A runtime exception applies to
00014  *  this software (see LICENSE.GPL_WRE file mentioned below for details).
00015  *
00016  *  This program is distributed in the hope that it will be useful,
00017  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00018  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019  *  GNU Library General Public License for more details.
00020  *
00021  *  Read the full text in the LICENSE.GPL_WRE file in the doc directory.
00022  */
00023
00024 #include <filters/laplace.h>
00025
00026 #include <core/exception.h>
00027
00028 #include <ippi.h>
00029 #include <cmath>
00030 #include <cstdlib>
00031 
00032 /** @class FilterLaplace <filters/laplace.h>
00033  * Laplacian filter.
00034  * Laplacian of Gaussian filter.
00035  * @author Tim Niemueller
00036  */
00037 
00038 /** Constructor. */
00039 FilterLaplace::FilterLaplace()
00040   : Filter("FilterLaplace")
00041 {
00042   kernel = NULL;
00043 }
00044
00045 
00046 /** Constructor.
00047  * @param sigma sigma for Laplacian
00048  * @param size size of kernel
00049  * @param scale scale factor
00050  */
00051 FilterLaplace::FilterLaplace(float sigma, unsigned int size, float scale)
00052   : Filter("FilterLaplace")
00053 {
00054   kernel_size = size;
00055   kernel = (int *)malloc( size * size * sizeof(int) );
00056   calculate_kernel( kernel, sigma, size, scale );
00057 }
00058
00059 
00060 /** Destructor. */
00061 FilterLaplace::~FilterLaplace()
00062 {
00063   if ( kernel != NULL ) {
00064     free( kernel );
00065   }
00066 }
00067
00068
00069 void
00070 FilterLaplace::apply()
00071 {
00072   IppiSize size;
00073   size.width = src_roi[0]->width - kernel_size;
00074   size.height = src_roi[0]->height - kernel_size;
00075
00076   IppStatus status;
00077
00078   if ( kernel == NULL ) {
00079     //                                   base + number of bytes to line y              + pixel bytes
00080     status = ippiFilterLaplace_8u_C1R( src[0] + (src_roi[0]->start.y * src_roi[0]->line_step) + (src_roi[0]->start.x * src_roi[0]->pixel_step), src_roi[0]->line_step,
00081                                        dst + (dst_roi->start.y * dst_roi->line_step) + (dst_roi->start.x * dst_roi->pixel_step), dst_roi->line_step,
00082                                        size, ippMskSize5x5 );
00083   } else {
00084     IppiSize ksize = { kernel_size, kernel_size };
00085     IppiPoint kanchor = { (kernel_size + 1) / 2, (kernel_size + 1) / 2 };
00086
00087     /*
00088     std::cout << "steps:   " << src_roi[0]->line_step << "   " << dst_roi->line_step << std::endl
00089               << "ksize:   " << ksize.width << " x " << ksize.height << std::endl
00090               << "kanchor: " << kanchor.x << "," << kanchor.y << std::endl;
00091     */
00092
00093     status = ippiFilter_8u_C1R( src[0] + ((src_roi[0]->start.y + kernel_size / 2) * src_roi[0]->line_step) + ((src_roi[0]->start.x + kernel_size / 2) * src_roi[0]->pixel_step), src_roi[0]->line_step,
00094                                 dst + ((dst_roi->start.y + kernel_size / 2) * dst_roi->line_step) + ((dst_roi->start.x + kernel_size / 2) * dst_roi->pixel_step), dst_roi->line_step,
00095                                 size, kernel, ksize, kanchor, 1 );
00096
00097   }
00098
00099   if ( status != ippStsNoErr ) {
00100     throw fawkes::Exception("Laplace filter failed with %i\n", status);
00101   }
00102
00103   /*
00104   std::cout << "FilterLaplace: ippiFilterLaplace exit code: " << std::flush;
00105   switch (status) {
00106   case ippStsNoErr:
00107     std::cout << "ippStsNoErr";
00108     break;
00109   case ippStsNullPtrErr:
00110     std::cout << "ippStsNullPtrErr";
00111     break;
00112   case ippStsSizeErr:
00113     std::cout << "ippStsSizeErr";
00114     break;
00115   case ippStsStepErr:
00116     std::cout << "ippStsStepErr";
00117     break;
00118   case ippStsMaskSizeErr:
00119     std::cout << "ippStsMaskSizeErr";
00120     break;
00121   default:
00122     std::cout << "Unknown status " << status;
00123   }
00124   std::cout << std::endl;
00125   */
00126 }
00127
00128 
00129 /** Calculate a Laplacian of Gaussian kernel.
00130  * The kernel is calculated with the formula
00131  * \f[
00132  *   roundf( \frac{-1}{\pi * \sigma^4} * 
00133  *           ( 1 - \frac{w^2 + h^2}{2 * \sigma^2} )
00134  *           * e^{-\frac{w^2 + h^2}{2 * \sigma^2}} * \mathtt{scale} )
00135  * \f]                     
00136  *
00137  * @param kernel buffer contains kernel upon return
00138  * @param sigma sigma for formula
00139  * @param size kernel is of quadratic size \f$\mathtt{size} \times \mathtt{size}\f$
00140  * @param scale scale parameter in formula
00141  */
00142 void
00143 FilterLaplace::calculate_kernel(int *kernel, float sigma, unsigned int size, float scale)
00144 {
00145   //  title "LoGFUNC__________________________________________"
00146
00147   /*
00148   std::cout.precision( 5 );
00149   std::cout.width( 10 );
00150 
00151   std::cout << "Discrete Laplacian kernel for sigma=" << sigma
00152             << " quadratic size of " << size
00153             << " scaled by " << scale << std::endl;
00154   */
00155   for (int h = (-(int)(size / 2)); h <= (int)((size - 1) / 2); ++h) {
00156     for (int w = (-(int)(size / 2)); w <= (int)((size - 1) / 2); ++w) {
00157       //float v = ( (w*w + h*h - 2 * sigma * sigma) / sigma * sigma * sigma * sigma )
00158       //* exp( -( (w*w + h*h) / (2 * sigma * sigma) ));
00159       int v =  (int)roundf( - 1/( M_PI * sigma * sigma * sigma * sigma ) *
00160                            ( 1 - ( (w*w + h*h) / (2 * sigma * sigma) ) )
00161                            * exp( -( (w*w + h*h) / (2 * sigma * sigma) )) * scale  );
00162       // std::cout << "   " << v << std::flush;
00163       kernel[ (h + (size / 2)) * size + (w + (size / 2)) ] = v;
00164     }
00165     //std::cout << std::endl;
00166   }
00167
00168   /*
00169   for (int h = 0; h < size; ++h) {
00170     for (int w = 0; w < size; ++w) {
00171       std::cout << "   " << kernel[ h * size + w ] << std::flush;
00172     }
00173     std::cout << std::endl;
00174   }
00175   */
00176
00177 }