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