env_c_math_ops.c

Go to the documentation of this file.

/*!@file Envision/env_c_math_ops.c Fixed-point integer math versions of some of our floating-point image functions */
 
// //////////////////////////////////////////////////////////////////// //
// The iLab Neuromorphic Vision C++ Toolkit - Copyright (C) 2000-2005   //
// by the University of Southern California (USC) and the iLab at USC.  //
// See http://iLab.usc.edu for information about this project.          //
// //////////////////////////////////////////////////////////////////// //
// Major portions of the iLab Neuromorphic Vision Toolkit are protected //
// under the U.S. patent ``Computation of Intrinsic Perceptual Saliency //
// in Visual Environments, and Applications'' by Christof Koch and      //
// Laurent Itti, California Institute of Technology, 2001 (patent       //
// pending; application number 09/912,225 filed July 23, 2001; see      //
// http://pair.uspto.gov/cgi-bin/final/home.pl for current status).     //
// //////////////////////////////////////////////////////////////////// //
// This file is part of the iLab Neuromorphic Vision C++ Toolkit.       //
//                                                                      //
// The iLab Neuromorphic Vision C++ Toolkit is free software; you can   //
// redistribute it and/or modify it under the terms of the GNU General  //
// Public License as published by the Free Software Foundation; either  //
// version 2 of the License, or (at your option) any later version.     //
//                                                                      //
// The iLab Neuromorphic Vision C++ Toolkit 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 General Public License for more details.       //
//                                                                      //
// You should have received a copy of the GNU General Public License    //
// along with the iLab Neuromorphic Vision C++ Toolkit; if not, write   //
// to the Free Software Foundation, Inc., 59 Temple Place, Suite 330,   //
// Boston, MA 02111-1307 USA.                                           //
// //////////////////////////////////////////////////////////////////// //
//
// Primary maintainer for this file: Rob Peters <rjpeters at usc dot edu>
// $HeadURL: svn://isvn.usc.edu/software/invt/trunk/saliency/src/Envision/env_c_math_ops.c $
// $Id: env_c_math_ops.c 11331 2009-06-23 17:57:49Z itti $
//
 
#include <jevoisbase/Components/Saliency/env_c_math_ops.h>
 
#include <jevoisbase/Components/Saliency/env_log.h>
#include <jevoisbase/Components/Saliency/env_types.h>
 
// ######################################################################
void env_c_lowpass_5_x_dec_x_fewbits_optim(const intg32* src, const env_size_t w, const env_size_t h,
                                           intg32* dst, const env_size_t w2)
{
  ENV_ASSERT(w2 == w/2);
  
  if (w == 2 || w == 3) //////////////////////////////////////////////////
    for (env_size_t j = 0; j < h; ++j)
    {
      // leftmost point  [ (6^) 4 ] / 10
      *dst++ = (src[0] * 3 + src[1] * 2) / 5;
      src += w;  // src back to same position as dst
    }
  else  ////////////////////////////// general case for width() >= 4
    // *** unfolded version (all particular cases treated) for
    // max speed.
    // notations: in () is the position of dest ptr, and ^ is src ptr
    // ########## horizontal pass
    for (env_size_t j = 0; j < h; ++j)
    {
      const intg32* src2 = src;
      // leftmost point  [ (8^) 4 ] / 12
      *dst++ = (src2[0] * 2 + src2[1]) / 3;
      
      // skip second point
      
      // rest of the line except last 2 points  [ .^ 4 (8) 4 ] / 16
      for (env_size_t i = 0; i < w-3; i += 2)
      {
        *dst++ = (src2[1] + src2[3] + src2[2] * 2) >> 2;
        src2 += 2;
      }
      
      src += w;
    }
}
 
// ######################################################################
void env_c_lowpass_5_y_dec_y_fewbits_optim(const intg32* src, const env_size_t w, const env_size_t h,
                                           intg32* dst, const env_size_t h2)
{
  ENV_ASSERT(h2 == h/2);
  
  // ########## vertical pass  (even though we scan horiz for speedup)
  const env_size_t w2 = w + w, w3 = w2 + w; // speedup
  
  if (h == 2 || h == 3) //////////////////////////////////////////////////
  {
    // topmost points  ( [ (6^) 4 ] / 10 )^T
    for (env_size_t i = 0; i < w; ++i)
    {
      *dst++ = (src[0] * 3 + src[w] * 2) / 5;
      src++;
    }
    src -= w;  // go back to top-left
  }
  else  ///////////////////////////////// general case for height >= 4
  {
    // topmost points  ( [ (8^) 4 ] / 12 )^T
    for (env_size_t k = 0; k < w; ++k)
    {
      *dst++ = (src[ 0] * 2 + src[ w]) / 3;
      src++;
    }
    src -= w;  // go back to top-left
    
    // second point skipped
    
    // rest of the column except last 2 points ( [ .^ 4 (8) 4 ] / 16 )T
    for (env_size_t i = 0; i < h-3; i += 2)
    {
      for (env_size_t k = 0; k < w; ++k)
      {
        *dst++ = (src[ w] + src[w3] + src[w2] * 2) >> 2;
        src++;
      }
      src += w;
    }
  }
}
 
// ######################################################################
void env_c_lowpass_9_x_fewbits_optim(const intg32* src, const env_size_t w, const env_size_t h, intg32* dst)
{
  ENV_ASSERT(w >= 9);
  
  // boundary conditions: truncated filter
  for (env_size_t j = 0; j < h; ++j)
  {
    // leftmost points
    *dst++ =                  // [ (72^) 56 28 8 ]
      (src[0] * 72 +
       src[1] * 56 +
       src[2] * 28 +
       src[3] *  8
       ) / 164;
    *dst++ =                  // [ 56^ (72) 56 28 8 ]
      ((src[0] + src[2]) * 56 +
       src[1] * 72 +
       src[3] * 28 +
       src[4] *  8
       ) / 220;
    *dst++ =                  // [ 28^ 56 (72) 56 28 8 ]
      ((src[0] + src[4]) * 28 +
       (src[1] + src[3]) * 56 +
       src[2] * 72 +
       src[5] *  8
       ) / 248;
    
    // far from the borders
    for (env_size_t i = 0; i < w - 6; ++i)
    {
      *dst++ =              // [ 8^ 28 56 (72) 56 28 8 ]
        ((src[0] + src[6]) *  8 +
         (src[1] + src[5]) * 28 +
         (src[2] + src[4]) * 56 +
         src[3] * 72
         ) >> 8;
      ++src;
    }
    
    // rightmost points
    *dst++ =                  // [ 8^ 28 56 (72) 56 28 ]
      (src[0] *  8 +
       (src[1] + src[5]) * 28 +
       (src[2] + src[4]) * 56 +
       src[3] * 72
       ) / 248;
    ++src;
    *dst++ =                  // [ 8^ 28 56 (72) 56 ]
      (src[0] *  8 +
       src[1] * 28 +
       (src[2] + src[4]) * 56 +
       src[3] * 72
       ) / 220;
    ++src;
    *dst++ =                  // [ 8^ 28 56 (72) ]
      (src[0] *  8 +
       src[1] * 28 +
       src[2] * 56 +
       src[3] * 72
       ) / 164;
    src += 4;  // src back to same as dst (start of next line)
  }
}
 
// ######################################################################
void env_c_lowpass_9_y_fewbits_optim(const intg32* src, const env_size_t w, const env_size_t h, intg32* dst)
{
  ENV_ASSERT(h >= 9);
  
  // index computation speedup:
  const env_size_t w2 = w + w, w3 = w2 + w, w4 = w3 + w, w5 = w4 + w, w6 = w5 + w;
  
  // *** vertical pass ***
  for (env_size_t i = 0; i < w; ++i)
  {
    *dst++ =
      (src[ 0] * 72 +
       src[ w] * 56 +
       src[w2] * 28 +
       src[w3] *  8
       ) / 164;
    ++src;
  }
  src -= w; // back to top-left
  for (env_size_t i = 0; i < w; ++i)
  {
    *dst++ =
      ((src[ 0] + src[w2]) * 56 +
       src[ w] * 72 +
       src[w3] * 28 +
       src[w4] *  8
       ) / 220;
    ++src;
  }
  src -= w; // back to top-left
  for (env_size_t i = 0; i < w; ++i)
  {
    *dst++ =
      ((src[ 0] + src[w4]) * 28 +
       (src[ w] + src[w3]) * 56 +
       src[w2] * 72 +
       src[w5] *  8
       ) / 248;
    ++src;
  }
  src -= w; // back to top-left
  
  for (env_size_t j = 0; j < h - 6; j ++)
    for (env_size_t i = 0; i < w; ++i)
    {
      *dst++ =
        ((src[ 0] + src[w6]) *  8 +
         (src[ w] + src[w5]) * 28 +
         (src[w2] + src[w4]) * 56 +
         src[w3]  * 72
         ) >> 8;
      ++src;
    }
  
  for (env_size_t i = 0; i < w; ++i)
  {
    *dst++ =
      (src[ 0] *  8 +
       (src[ w] + src[w5]) * 28 +
       (src[w2] + src[w4]) * 56 +
       src[w3] * 72
       ) / 248;
    ++src;
  }
  for (env_size_t i = 0; i < w; ++i)
  {
    *dst++ =
      (src[ 0] *  8 +
       src[ w] * 28 +
       (src[w2] + src[w4]) * 56 +
       src[w3] * 72
       ) / 220;
    ++src;
  }
  for (env_size_t i = 0; i < w; ++i)
  {
    *dst++ =
      (src[ 0] *  8 +
       src[ w] * 28 +
       src[w2] * 56 +
       src[w3] * 72
       ) / 164;
    ++src;
  }
}
 
// ######################################################################
void env_c_get_min_max(const intg32* src, const env_size_t sz, intg32* xmini, intg32* xmaxi)
{
  ENV_ASSERT(sz > 0);
  *xmini = *xmaxi = src[0];
  for (env_size_t i = 0; i < sz; ++i)
  {
    if (src[i] < *xmini) *xmini = src[i];
    else if (src[i] > *xmaxi) *xmaxi = src[i];
  }
}
 
// ######################################################################
void env_c_inplace_rectify(intg32* dst, const env_size_t sz)
{
  for (env_size_t i = 0; i < sz; ++i) if (dst[i] < 0) dst[i] = 0;
}
 
// ######################################################################
void env_c_inplace_normalize(intg32* const dst, const env_size_t sz, const intg32 nmin, const intg32 nmax,
                             intg32* const actualmin_p, intg32* const actualmax_p, const intg32 rangeThresh)
{
  ENV_ASSERT(sz > 0);
  ENV_ASSERT(nmax >= nmin);
  
  intg32 mi, ma;
  env_c_get_min_max(dst, sz, &mi, &ma);
  const intg32 old_scale = ma - mi;
  if (old_scale == 0 || old_scale < rangeThresh) // input image is uniform
  { for (env_size_t i = 0; i < sz; ++i) dst[i] = 0; return; }
  const intg32 new_scale = nmax - nmin;
  
  if (new_scale == 0) // output range is uniform
  { for (env_size_t i = 0; i < sz; ++i) dst[i] = nmin; return; }
  
  
  intg32 actualmin, actualmax;
  
  if (old_scale == new_scale)
  {
    const intg32 add = nmin - mi;
    if (add != 0) for (env_size_t i = 0; i < sz; ++i) dst[i] += add;
    
    actualmin = nmin;
    actualmax = nmax;
  }
#if defined(ENV_INTG64_TYPE)
  else
  {
    for (env_size_t i = 0; i < sz; ++i)
      dst[i] = nmin + ((((intg64)(dst[i] - mi)) * new_scale) / old_scale);
    
    actualmin = nmin;
    actualmax = nmin + ((((intg64)(old_scale)) * new_scale) / old_scale);
  }
#else
  else if (old_scale > new_scale)
  {
    // we want to do new = (old*new_scale)/oscale; however, the old*new_scale part might overflow if old or new_scale is
    // too large, so let's reduce both new_scale and oscale until new_scale*old won't overflow
    
    const intg32 nscale_max = INTG32_MAX / old_scale;
    
    intg32 new_scale_reduced = new_scale;
    intg32 old_scale_reduced = old_scale;
    
    while (new_scale_reduced > nscale_max)
    {
      new_scale_reduced /= 2;
      old_scale_reduced /= 2;
    }
    
    ENV_ASSERT(new_scale_reduced >= 1);
    ENV_ASSERT(old_scale_reduced >= 1);
    
    for (env_size_t i = 0; i < sz; ++i)
      dst[i] = nmin + (((dst[i] - mi) * new_scale_reduced) / (old_scale_reduced+1));
    
    actualmin = nmin;
    actualmax = nmin + ((old_scale * new_scale_reduced) / (old_scale_reduced+1));
  }
  else // (old_scale < new_scale)
  {
    const intg32 mul = new_scale / old_scale;
    for (env_size_t i = 0; i < sz; ++i) dst[i] = nmin + ((dst[i] - mi) * mul);
    
    actualmin = nmin;
    actualmax = nmin + (old_scale * mul);
  }
#endif // !defined(ENV_INTG64_TYPE)
  
  // Don't assign to the pointers until the very end, in case
  // the user passes pointers to nmin,nmax as the
  // actualmin/actualmax pointers
  if (actualmin_p) *actualmin_p = actualmin;
  if (actualmax_p) *actualmax_p = actualmax;
}
 
// ######################################################################
void env_c_luminance_from_byte(const struct env_rgb_pixel* const src, const env_size_t sz,
                               const env_size_t nbits, intg32* const dst)
{
  if (nbits > 8)
    for (env_size_t i = 0; i < sz; ++i)
      dst[i] = ((src[i].p[0] + src[i].p[1] + src[i].p[2]) / 3) << (nbits - 8);
  else if (nbits < 8)
    for (env_size_t i = 0; i < sz; ++i)
      dst[i] = ((src[i].p[0] + src[i].p[1] + src[i].p[2]) / 3) >> (8 - nbits);
  else // nbits == 8
    for (env_size_t i = 0; i < sz; ++i)
      dst[i] = (src[i].p[0] + src[i].p[1] + src[i].p[2]) / 3;
}
 
// ######################################################################
void env_c_image_div_scalar(const intg32* const a, const env_size_t sz, intg32 val, intg32* const dst)
{
  for (env_size_t i = 0; i < sz; ++i) dst[i] = a[i] / val;
}
 
// ######################################################################
void env_c_image_div_scalar_accum(const intg32* const a, const env_size_t sz, intg32 val, intg32* const dst)
{
  for (env_size_t i = 0; i < sz; ++i) dst[i] += a[i] / val;
}
 
// ######################################################################
void env_c_image_minus_image(const intg32* const a, const intg32* const b, const env_size_t sz, intg32* const dst)
{
  for (env_size_t i = 0; i < sz; ++i) dst[i] = a[i] - b[i];
}