DericheGradientImage.C

Go to the documentation of this file.

/*---------------------------------------------------------------------+
 | Library QgarLib, graphics analysis and recognition                  |
 | Copyright (C) 2002  Qgar Project, LORIA                             |
 |                                                                     |
 | This library is free software; you can redistribute it and/or       |
 | modify it under the terms of the GNU Lesser General Public          |
 | License version 2.1, as published by the Free Software Foundation.  |
 |                                                                     |
 | This library 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.         |
 |                                                                     |
 | The GNU Lesser General Public License is included in the file       |
 | LICENSE.LGPL, in the root directory of the Qgar packaging. See      |
 | http://www.gnu.org/licenses/lgpl.html for the terms of the licence. |
 | To receive a paper copy, write to the Free Software Foundation,     |
 | Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.       |
 |                                                                     |
 | Contact Project Qgar for any information:                           |
 |   LORIA - équipe Qgar                                               |
 |   B.P. 239, 54506 Vandoeuvre-lès-Nancy Cedex, France                |
 |   email: qgar-contact@loria.fr                                      |
 |   http://www.qgar.org/                                              |
 *---------------------------------------------------------------------*/


/**
 * @file   DericheGradientImage.C
 * @brief  Implementation of class qgar::DericheGradientImage.
 *
 *         See file DericheGradientImage.H for the interface.
 *
 * @author <a href="mailto:qgar-develop@loria.fr?subject=Qgar fwd Karl Tombre">Karl Tombre</a>
 * @date   July 3, 2001  11:35
 * @since  Qgar 1.0
 */


// STL
#include <cmath>
// QGAR
#include <qgarlib/DericheGradientImage.H>
#include <qgarlib/GenImage.H>


namespace qgar
{


// -------------------------------------------------------------------
// C O N S T R U C T O R S
// -------------------------------------------------------------------


// CONSTRUCT FROM GIVEN GREY-LEVEL IMAGE

DericheGradientImage::DericheGradientImage(const GreyLevelImage& img,
                                           float alpha,
                                           float w,
                                           int normalized)

    : AbstractGradientImage(img),
      _alpha(alpha),
      _w(w)
    
{
  PRIVATEbuildDericheImage(normalized);
}


// CONSTRUCT FROM A FLOAT IMAGE

DericheGradientImage::DericheGradientImage(const FloatImage& img,
                                           float alpha,
                                           float w,
                                           int normalized)

    : AbstractGradientImage(img),
      _alpha(alpha),
      _w(w)

  {
    PRIVATEbuildDericheImage(normalized);
  }


// -------------------------------------------------------------------
// A U X I L I A R I E S
// -------------------------------------------------------------------


void
DericheGradientImage::PRIVATEbuildDericheImage(int normalized)
{
  float c0 = _alpha * _alpha + _w * _w;

  // normalization coefficient for smoothing
  float k = (1 - 2 * (float) (exp((double) (- _alpha)) * cos((double) _w))
             + (float) exp((double) (-2 * _alpha))) * c0 /
    (2 * _alpha * (float) (exp((double) (- _alpha)) * sin((double) _w))
     + _w * (1 - (float) exp((double) (-2 * _alpha))));

  // normalization coefficient for derivation
  a0 = k * _w / c0;
  float c;
  if (normalized)
    c = k * c0;
  else
    c = (1 - 2 * (float) (exp((double)(-_alpha)) * cos((double) _w))
         + (float) exp((double) (-2 * _alpha))) /
      ((float) (exp((double) (-_alpha)) * sin((double) _w)));

  float c1 = k * _alpha / c0;

  // Derivation coefficients
  b1 = (-2.f) * (float) (exp((double) (-_alpha)) * cos((double) _w));
  b2 = (float) exp((double) (-2. * _alpha));

  // Smoothing coefficients
  a1 = (- a0 * (float) cos((double) _w) + c1 * (float) sin((double) _w)) * 
    (float) exp((double) (-_alpha));
  a2 = a1 - a0 * b1;
  a3 = (- a0) * b2;

  a = (- c) * (float) (exp((double) (- _alpha)) * sin((double) _w));

  // Filter in x
  dfilterX();
  // Filter in y
  dfilterY();

  // maybe set some lines and columns to 0 ???
}



// X FILTER
// ========


void
DericheGradientImage::dfilterX()
{
  float* mem = new float [_width];
  float* buff1 = new float [_width];
  float* buff2 = new float [_width];
  float* qstart = _dxImg.pPixMap();
  int i;
  int j;

  // Derive along x axis
  for (i = 0 ; i < _height ; ++i, qstart += _width)
    {
      // From left to right
      float* p = buff1;
      // Remember that _dxImg and _dyImg are initialized
      // with original values
      float* q = qstart;
      float* pmem = mem;

      // First two pixels
      *p = (*q);
      *pmem++ = a * (*q);
      ++p;
      ++q;
      *p = (*q);
      *pmem++ = a * (*q);
      ++p;
      ++q;

      // And loop on rest with recurrence formula
      for (j = 2 ; j < _width ; ++j, ++p, ++q, ++pmem)
        {
          *pmem = a * (*q);
          *p = (*(pmem - 1)) - b1 * (*(p-1)) - b2 * (*(p-2));
        }

      // From right to left
      p = buff2 + _width - 1;
      q = qstart + _width - 1;
      
      // Two last
      *p-- = (*q--);
      *p-- = (*q--);
      
      // And loop on rest with recurrence formula
      pmem = mem + _width - 3;
      for (j = _width - 3 ; j >= 0 ;  --j, --p, --pmem)
        {
          *p = (- *(pmem + 1)) - b1 * (*(p+1)) - b2 * (*(p+2));
        }
      
      // Sum the two results
      float* r = qstart;
      p = buff1;
      q = buff2;
      for (j = 0 ; j < _width ; ++j)
        {
          *r++ = (*p++) + *q++;
        }
    }

  delete [] mem;
  delete [] buff1;
  delete [] buff2;

  buff1 = new float [_height];
  buff2 = new float [_height];

  // Smooth along y axis
  qstart = _dxImg.pPixMap();
  for (j = 0 ; j < _width ; j++, qstart++)
    {
      // From top to bottom
      float* p = buff1;
      float* q = qstart;

      // First two
      *p = (*q);
      p++;
      q += _width;
      *p = (*q);
      p++;
      q += _width;

      // Loop with recurrence formula
      for (i = 2 ; i < _height ; i++, p++, q += _width)
        *p = a0 * (*q) + a1 * *(q - _width) - b1 * (*(p-1)) - b2 * (*(p-2));

      // From bottom to top
      p = buff2 + _height - 1;
      q = qstart + (_height - 1) * _width;

      // Two last
      *p = (*q);
      p--; q -= _width;
      *p = (*q);
      p--; q -= _width;
      
      // Loop on rest with recurrence formula
      for (i = _height - 3 ; i >= 0 ; i--, p--, q -= _width)
        *p = a2 * (*(q + _width)) + a3 * (*(q + (_width + _width)))
          - b1 * (*(p+1)) - b2 * (*(p+2));

      // Sum up the two results
      float* r = qstart;
      p = buff1;
      q = buff2;
      for (i = 0 ; i < _height ; i++, r += _width)
        *r = (*p++) + *q++;
    }

  delete [] buff1;
  delete [] buff2;
}


// Y FILTER
// ========


void
DericheGradientImage::dfilterY()
{
  float* mem = new float [_height];
  float* buff1 = new float [_height];
  float* buff2 = new float [_height];
  int i;
  int j;

  // Derive along y axis
  float* qstart = _dyImg.pPixMap();

  for (j = 0 ; j < _width ; ++j, ++qstart)
    {
      // From top to bottom
      float* p = buff1;
      // Remember that _dxImg and _dyImg are initialized
      // with original values
      float* q = qstart;

      float* pmem = mem;

      // First two
      *p = (*q);
      *pmem++ = a * (*q);
      ++p;
      q += _width;
      *p = (*q);
      *pmem++ = a * (*q);
      ++p;
      q += _width;

      // Loop with recurrence formula
      for (i = 2 ; i < _height ; ++i, ++p, q += _width, ++pmem)
        {
          *pmem = a * (*q);
          *p = (*(pmem - 1)) - b1 * (*(p-1)) - b2 * (*(p-2));
        }

      // From bottom to top
      p = buff2 + _height - 1;
      q = qstart + (_height - 1) * _width;

      // Two last
      *p = (*q);
      --p;
      q -= _width;
      *p = (*q);
      --p;
      q -= _width;
      
      // Loop on rest with recurrence formula
      pmem = mem + _height - 3;
      for (i = _height - 3 ; i >= 0 ;--i, --p, --pmem)
        {
          *p = (- *(pmem + 1)) - b1 * (*(p+1)) - b2 * (*(p+2));
        }

      // Sum up the two results
      float* r = qstart;
      p = buff1;
      q = buff2;
      for (i = 0 ; i < _height ; ++i, r += _width)
        {
          *r = (*p++) + *q++;
        }
    } // END for j

  delete [] mem;
  delete [] buff1;
  delete [] buff2;

  buff1 = new float [_width];
  buff2 = new float [_width];

  qstart = _dyImg.pPixMap();
  // Smooth along x axis
  for (i = 0 ; i < _height ; ++i, qstart += _width)
    {
      // From left to right
      float* p = buff1;
      float* q = qstart;

      // First two pixels
      *p++ = (*q++);
      *p++ = (*q++);

      // And loop on rest with recurrence formula
      for (j = 2 ; j < _width ; ++j, ++p, ++q)
        {
          *p = a0 * (*q) + a1 * (*(q - 1)) - b1 * (*(p-1)) - b2 * (*(p-2));
        }

      // From right to left
      p = buff2 + _width - 1;
      q = qstart + _width - 1;
      
      // Two last
      *p-- = (*q--);
      *p-- = (*q--);
      
      // And loop on rest with recurrence formula
      for (j = _width - 3 ; j >= 0 ;  --j, --p, --q)
        {
          *p = a2 * (*(q + 1)) + a3 * (*(q + 2)) - b1 * (*(p+1)) - b2 * (*(p+2));
        }

      // Sum the two results
      float* r = qstart;
      p = buff1;
      q = buff2;
      for (j = 0 ; j < _width ; ++j)
        {
          *r++ = (*p++) + *q++;
        }}


  delete [] buff1;
  delete [] buff2;
}

// -------------------------------------------------------------------

} // namespace qgar