segmentImageMC2.C

Go to the documentation of this file.

/*!@file VFAT/segmentImageMC2.C Basic image segmenter blob finder using color */

// //////////////////////////////////////////////////////////////////// //
// The iLab Neuromorphic Vision C++ Toolkit - Copyright (C) 2001 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: T. Nathan Mundhenk <mundhenk@usc.edu>
// $HeadURL: svn://isvn.usc.edu/software/invt/trunk/saliency/src/VFAT/segmentImageMC2.C $
// $Id: segmentImageMC2.C 10794 2009-02-08 06:21:09Z itti $
//

// ############################################################
// ############################################################
// ##### --- VFAT ---
// ##### Vision Feature Analysis Tool:
// ##### T. Nathan Mundhenk nathan@mundhenk.com
// ##### Laurent Itt itti@pollux.usc.edu
// #####
// ############################################################
// ############################################################

#ifndef SEGMENTIMAGEMC2_C_DEFINED
#define SEGMENTIMAGEMC2_C_DEFINED

#include "Util/Assert.H"
#include "VFAT/segmentImageMC2.H"
#include <iostream>
#include <vector>
#include <cstdio>
#include <cstdlib>

/*********************************************************************/
/*********************************************************************/
/*        PUBLIC methods                                             */
/*********************************************************************/
/*********************************************************************/

/*********************************************************************/

template SI_TEMPLATE_CLASS
segmentImageMC2<SI_TEMPLATE>::segmentImageMC2()
{
  SI_set1 = false; SI_set2 = false; SI_set3 = false; SI_set4 = false;
  LINFO("CREATED");
  SI_lowThresh.resize(SI_channels,0);
  SI_highThresh.resize(SI_channels,0);
  SI_maxIDVal = 0;
  SI_killVal  = 2;
  Image<FLOAT> timage;
  SI_infeatureMaps.resize(SI_channels,timage);
  SI_useCandidateBandPass = true;
  SI_removeSingles        = true;
  SI_set4                 = true;
}

/*********************************************************************/

template SI_TEMPLATE_CLASS
segmentImageMC2<SI_TEMPLATE>::~segmentImageMC2()
{}

/* set color segmentation values for color and
   threshold
*/

/*********************************************************************/

template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIsetVal(const typename
                                           std::vector<FLOAT> &val,
                                           const typename
                                           std::vector<FLOAT> &thresh)
{
  typename std::vector<FLOAT>::iterator ilowThresh    = SI_lowThresh.begin();
  typename std::vector<FLOAT>::iterator ihighThresh   = SI_highThresh.begin();
  typename std::vector<FLOAT>::const_iterator ival    = val.begin();
  typename std::vector<FLOAT>::const_iterator ithresh = thresh.begin();

  for(INT i = 0; i < val.size(); i++, ++ilowThresh, ++ihighThresh,
        ++ival, ++ithresh)
  {
    *ilowThresh  = *ival - *ithresh;
    *ihighThresh = *ival + *ithresh;
  }
  SI_set1 = true;
}

/*********************************************************************/

template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIsetValThresh(const typename
                                                 std::vector<FLOAT> &high,
                                                 const typename
                                                 std::vector<FLOAT> &low)
{
  SI_lowThresh  = low;
  SI_highThresh = high;

  typename std::vector<FLOAT>::iterator ilow  = SI_lowThresh.begin();
  typename std::vector<FLOAT>::iterator ihigh = SI_highThresh.begin();

  // Check to make sure low is less than high for each feature
  while(ilow != SI_lowThresh.end())
  {
    if(*ilow >= *ihigh)
    {
      LFATAL("Low threshold must be less than high theshold. Got L %f H %f",*ilow,*ihigh);
    }
    ++ilow; ++ihigh;
  }

  SI_set1 = true;
}
/*********************************************************************/

template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIresetCandidates(const bool whichWay)
{
  Image<bool>::iterator
    candidatePixelsIter    = SI_candidatePixels.beginw();

  Image<bool>::iterator
    preCandidatePixelsIter = SI_preCandidatePixels.beginw();

  // set initial candidates for each pixel

  while(candidatePixelsIter != SI_candidatePixels.endw())
  {
    *candidatePixelsIter    = false;
    *preCandidatePixelsIter = whichWay;
    ++candidatePixelsIter; ++preCandidatePixelsIter;
  }
}

/*********************************************************************/
/* set size of window frame to inspect in image */

template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIsetFrame(int *x, int *y)
{
  SI_masterVec.resize(*x*(*y),-1);
  SI_reOrderVec.resize(*x*(*y));
  SI_reverseOrderVec.resize(*x*(*y));
  SI_centerX.resize(*x*(*y));
  SI_centerY.resize(*x*(*y));
  SI_Xsum.resize(*x*(*y));
  SI_Ysum.resize(*x*(*y));
  SI_mass.resize(*x*(*y));
  SI_xmin.resize(*x*(*y));
  SI_xmax.resize(*x*(*y));
  SI_ymin.resize(*x*(*y));
  SI_ymax.resize(*x*(*y));
  SI_reset.resize(*x*(*y));
  LINFO("SETTING WIDTH %d, HEIGHT %d",*x,*y);
  SI_blobID.resize(*x,*y,-1);
  SI_candidatePixels.resize(*x,*y,false);
  SI_preCandidatePixels.resize(*x,*y,false);
  SI_set2 = true;
}

/*********************************************************************/

template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIsetAvg(const INT doAvg)
{
  typename std::vector<FLOAT> temp(SI_channels,0);

  SI_avg.resize(doAvg,temp);
  SI_std.resize(doAvg,temp);
  SI_N.resize(doAvg,0);
  SI_tempAvg.resize(doAvg,0);
  SI_tempStd.resize(doAvg,0);
  SI_iter = doAvg;
  SI_count = 0;
  SI_set3 = true;
}

/*********************************************************************/

template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIresetAvg()
{
  typename std::vector<FLOAT> temp(SI_channels,0);

  for(typename std::vector<std::vector<FLOAT> >::iterator
        iavg = SI_avg.begin(); iavg != SI_avg.end(); ++iavg)
    *iavg = temp;
  for(typename std::vector<std::vector<FLOAT> >::iterator
        istd = SI_std.begin(); istd != SI_std.end(); ++istd)
    *istd = temp;

  typename std::vector<FLOAT>::iterator itempAvg = SI_tempAvg.begin();
  typename std::vector<FLOAT>::iterator itempStd = SI_tempStd.begin();
  typename std::vector<INT>::iterator   iN       = SI_N.begin();

  for(unsigned int i = 0; i < SI_N.size(); i++,
        ++itempAvg, ++itempStd, ++iN)
  {
    *itempAvg = 0.0F; *itempStd = 0.0F; *iN = 0;
  }
  SI_count = 0;

}

/*********************************************************************/
/*********************************************************************/
/* do image segmentation by calling private memebers to
   operate on image.
   1. Low pass image
   2. Find candidate pixels
   3. Eleminate single isolated candidates
   4. Link pixels in each blob
   5. Clean up
*/
/*********************************************************************/
/*********************************************************************/
template SI_TEMPLATE_CLASS
void segmentImageMC2<SI_TEMPLATE>::SIsegment(
                               Image<PixRGB<byte> > *image,
                               typename std::vector<Image<FLOAT> > *featureMap,
                               const bool lowPass)
{
  SI_workImage      = image;
  SIsegment(featureMap,lowPass);
}

/*********************************************************************/
template SI_TEMPLATE_CLASS
void segmentImageMC2<SI_TEMPLATE>::SIsegment(
                               typename std::vector<Image<FLOAT> > *featureMap,
                               const bool lowPass)
{
  struct timezone tz;
  struct timeval start, stop;
  tz.tz_minuteswest = 0;
  tz.tz_dsttime     = 0;
  gettimeofday(&start, &tz);

  // Low pass image if desired
  if(lowPass)
  {
    typename std::vector<Image<FLOAT> >::iterator
      iimage = featureMap->begin();
    typename std::vector<Image<FLOAT> >::iterator
      ifmap = SI_infeatureMaps.begin();

    while(iimage != featureMap->end())
    {
      *ifmap = lowPass5(*iimage);
      ++ifmap; ++iimage;
    }
    SI_featureMaps = &SI_infeatureMaps;
  }
  else
    SI_featureMaps = featureMap;


  // run doSegment which will call several consecutive methods
  SIdoSegment();

  gettimeofday(&stop,&tz);
}

/*********************************************************************/
template SI_TEMPLATE_CLASS
void segmentImageMC2<SI_TEMPLATE>::SItoggleCandidateBandPass(const bool toggle)
{
  SI_useCandidateBandPass = toggle;
}

/*********************************************************************/
template SI_TEMPLATE_CLASS
void segmentImageMC2<SI_TEMPLATE>::SItoggleRemoveSingles(const bool toggle)
{
  SI_removeSingles = toggle;
}

/*********************************************************************/
template SI_TEMPLATE_CLASS
void segmentImageMC2<SI_TEMPLATE>::SIsetKillValue(const unsigned int kv)
{
  SI_killVal = kv;
}

/*********************************************************************/
/* This method when called will take all remaning blobs from post
   processing and create a mother blob
*/
template SI_TEMPLATE_CLASS
Image<INT> segmentImageMC2<SI_TEMPLATE>::SIcreateMother(const Image<INT> &img) const
{
  Image<INT> mother;
  mother.resize(img.getWidth(),img.getHeight(),ZEROS);
  for(int x = SI_frameX1; x < SI_frameX2; x++)
  {
    for(int y = SI_frameY1; y < SI_frameY2; y++)
    {
      if(img.getVal(x,y) != 0)
        mother.setVal(x,y,1);
      else
        mother.setVal(x,y,0);
    }
  }
  return mother;
}

/*********************************************************************/
/* return blob map
 */
template SI_TEMPLATE_CLASS
Image<int> segmentImageMC2<SI_TEMPLATE>::SIreturnBlobs() const
{
  return SI_blobID;
}

/*********************************************************************/

template SI_TEMPLATE_CLASS
Image<bool> segmentImageMC2<SI_TEMPLATE>::SIreturnCandidates() const
{
  return SI_candidatePixels;
}

/*********************************************************************/

template SI_TEMPLATE_CLASS
Image<FLOAT> segmentImageMC2<SI_TEMPLATE>::SIreturnNormalizedCandidates() const
{
  Image<FLOAT> NC;
  NC.resize(SI_candidatePixels.getWidth(),SI_candidatePixels.getHeight());
  for(int x = 0; x < SI_candidatePixels.getWidth(); x++)
  {
    for(int y = 0; y < SI_candidatePixels.getHeight(); y++)
    {
      if(SI_candidatePixels.getVal(x,y))
        NC.setVal(x,y,255);
      else
        NC.setVal(x,y,0);
    }
  }
  return NC;
}

/*********************************************************************/

template SI_TEMPLATE_CLASS
Image<PixRGB<FLOAT> > segmentImageMC2<SI_TEMPLATE>::SIreturnWorkImage() const
{
  ASSERT((SI_doType == 1) || (SI_doType == 2));
  return *SI_workImage;
}

/*********************************************************************/

template SI_TEMPLATE_CLASS
INT segmentImageMC2<SI_TEMPLATE>::SInumberBlobs() const
{
  return SI_totalBlobs;
}

/*********************************************************************/

template SI_TEMPLATE_CLASS
std::vector<INT> segmentImageMC2<SI_TEMPLATE>::SIgetBlobMap() const
{
  return SI_reOrderVec;
}

/*********************************************************************/

template SI_TEMPLATE_CLASS
void segmentImageMC2<SI_TEMPLATE>::SIcalcMassCenter()
{

  Image<int>::iterator iblobID           = SI_blobID.beginw();
  Image<bool>::iterator icandidatePixels = SI_candidatePixels.beginw();

  const ushort width  = (ushort)SI_candidatePixels.getWidth();
  const ushort height = (ushort)SI_candidatePixels.getHeight();

  for(ushort y = 0 ; y < height; y++)
  {
    for(ushort x = 0 ; x < width; x++)
    {
      if((*icandidatePixels) && (*iblobID != -1))
      {
        //std::cerr << "foo " << reverseOrderVec[blobID.getVal(x,y)] << "\n";
        INT *indexBlob = &SI_reverseOrderVec[*iblobID];
        if(SI_reset[*indexBlob])
        {
          SI_reset[*indexBlob] = false;
          SI_Xsum[*indexBlob]  = x;
          SI_Ysum[*indexBlob]  = y;
          SI_mass[*indexBlob]  = 1;
          SI_xmin[*indexBlob]  = x;
          SI_ymin[*indexBlob]  = y;
          SI_xmax[*indexBlob]  = x;
          SI_ymax[*indexBlob]  = y;
        }
        else
        {
          SI_Xsum[*indexBlob] += x;
          SI_Ysum[*indexBlob] += y;
          SI_mass[*indexBlob]++;

          if(x <= SI_xmin[*indexBlob])
            SI_xmin[*indexBlob] = x;
          if(x >= SI_xmax[*indexBlob])
            SI_xmax[*indexBlob] = x;
          if(y <= SI_ymin[*indexBlob])
            SI_ymin[*indexBlob] = y;
          if(y >= SI_ymax[*indexBlob])
            SI_ymax[*indexBlob] = y;
        }
      }
      ++iblobID; ++icandidatePixels;
    }
  }

  for(INT b = 0; b < SI_totalBlobs; b++)
  {
    //std::cerr << "MASS " << mass[x] << "\n";
    if(SI_mass[b] > 0)
    {
      SI_centerX[b] = SI_Xsum[b]/SI_mass[b];
      SI_centerY[b] = SI_Ysum[b]/SI_mass[b];
    }
  }
}

/*********************************************************************/

template SI_TEMPLATE_CLASS
FLOAT segmentImageMC2<SI_TEMPLATE>::SIgetCenterX(const INT blob) const
{
  return SI_centerX[blob];
}

/*********************************************************************/

template SI_TEMPLATE_CLASS
FLOAT segmentImageMC2<SI_TEMPLATE>::SIgetCenterY(const INT blob) const
{
  return SI_centerY[blob];
}

/*********************************************************************/

template SI_TEMPLATE_CLASS
INT segmentImageMC2<SI_TEMPLATE>::SIgetMass(const INT blob) const
{
  return SI_mass[blob];
}

/*********************************************************************/

template SI_TEMPLATE_CLASS
int segmentImageMC2<SI_TEMPLATE>::SIgetXmin(const INT blob) const
{
  return SI_xmin[blob];
}

/*********************************************************************/
//! get X max for a blob
template SI_TEMPLATE_CLASS
int segmentImageMC2<SI_TEMPLATE>::SIgetXmax(const INT blob) const
{
  return SI_xmax[blob];
}

/*********************************************************************/
//! get Y min for a blob
template SI_TEMPLATE_CLASS
int segmentImageMC2<SI_TEMPLATE>::SIgetYmin(const INT blob) const
{
  return SI_ymin[blob];
}

/*********************************************************************/
//! get Y max for a blob
template SI_TEMPLATE_CLASS
int segmentImageMC2<SI_TEMPLATE>::SIgetYmax(const INT blob) const
{
  return SI_ymax[blob];
}

/*********************************************************************/
//! get the working image size in X
template SI_TEMPLATE_CLASS
int segmentImageMC2<SI_TEMPLATE>::SIgetImageSizeX() const
{
  return SI_candidatePixels.getWidth();
}

/*********************************************************************/
//! get the working image size in Y
template SI_TEMPLATE_CLASS
int segmentImageMC2<SI_TEMPLATE>::SIgetImageSizeY() const
{
  return SI_candidatePixels.getHeight();
}

/*********************************************************************/
template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIgetValue(INT *blob,
                                      typename std::vector<FLOAT> *mean,
                                      typename std::vector<FLOAT> *std,
                                      INT *in)
{
  typename std::vector<FLOAT>::iterator imean = mean->begin();
  typename std::vector<FLOAT>::iterator istd  = std->begin();
  typename std::vector<Image<FLOAT> >::iterator
    ifeatureMaps = SI_featureMaps->begin();

  bool dothis = true;
  *in = 0;
  for(INT i = 0; i < SI_featureMaps->size(); i++, ++imean, ++istd,
        ++ifeatureMaps)
  {
    Image<bool>::iterator icandidatePixels      = SI_candidatePixels.beginw();
    Image<int>::iterator iblobID               = SI_blobID.beginw();
    typename Image<FLOAT>::iterator iifeatureMaps = ifeatureMaps->beginw();

    FLOAT tot = 0;
    FLOAT ss = 0;

    while(icandidatePixels != SI_candidatePixels.endw())
    {
      if((*icandidatePixels) && (*iblobID != -1))
      {
        if(SI_reverseOrderVec[*iblobID] == *blob)
        {
          tot += *iifeatureMaps;
          ss += (pow(*iifeatureMaps,2))/SI_mass[*blob];
        }
      }
      ++icandidatePixels; ++iblobID; ++iifeatureMaps;
    }
    if(SI_mass[*blob] > 0)
    {
      *imean = tot/SI_mass[*blob];
      *istd = sqrt(fabs(ss - pow(*imean,2)));
      if(dothis == true)
        *in = SI_mass[*blob] + *in;
    }
    dothis = false;
  }
}

/*********************************************************************/

template SI_TEMPLATE_CLASS
void segmentImageMC2<SI_TEMPLATE>::SIgetValueMean(INT *blobListSize,
                                          typename std::vector<INT> *blobList,
                                          typename std::vector<FLOAT> *mean,
                                          typename std::vector<FLOAT> *stdd,
                                          FLOAT *mass)
{
  // average feature channels over N iterations for all motherblobs
  // modulus on dataset size
  if(SI_count == SI_iter)
    SI_count = 0;

  INT SI_tempN;

  typename std::vector<INT>::iterator iN                   = SI_N.begin();
  typename std::vector<std::vector<FLOAT> >::iterator iavg = SI_avg.begin();
  typename std::vector<std::vector<FLOAT> >::iterator istd = SI_std.begin();

  typename std::vector<FLOAT> *pavg = &SI_avg[SI_count];
  typename std::vector<FLOAT> *pstd = &SI_std[SI_count];
  INT *pcount = &SI_N[SI_count];
  *pcount = 0;

  // get mean values etc. for all blobs in this iteration
  for(INT i = 0; i < *blobListSize; i++)
  {
    SIgetValue(&blobList->at(i), &SI_tempAvg, &SI_tempStd, &SI_tempN);
    *pcount = SI_tempN + (*pcount);
    for(INT f = 0; f < SI_featureMaps->size(); f++)
    {
      pavg->at(f) += SI_tempAvg[f]*SI_tempN;
      pstd->at(f) += SI_tempStd[f]*SI_tempN;
    }
  }

  for(INT f = 0; f < SI_featureMaps->size(); f++)
  {
    pavg->at(f) = pavg->at(f)/(*pcount);
    pstd->at(f) = pstd->at(f)/(*pcount);
  }
  *mass = *pcount;
  SI_count++;

  INT massSum = 0;

  // combine the results over the last N iterations for everything
  for(INT c = 0; c < SI_iter; c++, ++iavg, ++istd, ++iN)
  {
    massSum += *iN;
    typename std::vector<FLOAT>::iterator iiavg = iavg->begin();
    typename std::vector<FLOAT>::iterator iistd = istd->begin();
    typename std::vector<FLOAT>::iterator imean = mean->begin();
    typename std::vector<FLOAT>::iterator istdd = stdd->begin();

    for(INT i = 0; i < SI_featureMaps->size(); i++, ++iiavg, ++iistd,
          ++imean, ++istdd)
    {
      if(*iN != 0)
      {
        *imean += *iiavg*(*iN);
        *istdd += *iistd*(*iN);
      }
    }
  }

  // find the total mean and std over all blobs and iterations
  typename std::vector<FLOAT>::iterator imean = mean->begin();
  typename std::vector<FLOAT>::iterator istdd = stdd->begin();
  for(INT i = 0; i < SI_featureMaps->size(); i++, ++imean, ++istdd)
  {
    *imean = *imean/massSum;
    *istdd = *istdd/massSum;
  }
}

/*********************************************************************/
/*********************************************************************/
/*        PRIVATE methods                                            */
/*********************************************************************/
/*********************************************************************/

/*********************************************************************/

/*! find candidate pixels by just checking image pixels to see if they
   are within threshold of desired feature values
   if so, set candidatePixels image pixel to true
*/
template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIfindCandidates()
{
  // FIX
  std::vector<int>::iterator imasterVec = SI_masterVec.begin();

  while(imasterVec != SI_masterVec.end())
    *imasterVec++ = -1;

  Image<bool>::iterator
    candidatePixelsIter    = SI_candidatePixels.beginw();

  Image<bool>::iterator
    preCandidatePixelsIter = SI_preCandidatePixels.beginw();

  // set initial candidates for each pixel

  while(candidatePixelsIter != SI_candidatePixels.endw())
  {
    if((*preCandidatePixelsIter) == true)
    {
      *candidatePixelsIter = true;
    }
    else
    {
      *candidatePixelsIter = false;
    }
    *preCandidatePixelsIter = true;
    ++candidatePixelsIter; ++preCandidatePixelsIter;
  }

  typename std::vector<Image<FLOAT> >::iterator ifeatureMaps;
  typename std::vector<FLOAT>::iterator ilowThresh  = SI_lowThresh.begin();
  typename std::vector<FLOAT>::iterator ihighThresh = SI_highThresh.begin();
  // Run over each feature channel

  for(ifeatureMaps = SI_featureMaps->begin();
      ifeatureMaps !=  SI_featureMaps->end(); ++ifeatureMaps,
        ++ilowThresh, ++ihighThresh)
  {
    typename Image<FLOAT>::iterator iifeatureMaps = ifeatureMaps->beginw();
    preCandidatePixelsIter = SI_preCandidatePixels.beginw();
    candidatePixelsIter    = SI_candidatePixels.beginw();

    // Run over each pixel

    while(iifeatureMaps != ifeatureMaps->endw())
    {
      if((*iifeatureMaps > *ihighThresh) || (*iifeatureMaps < *ilowThresh))
      {
        *preCandidatePixelsIter = false;
        *candidatePixelsIter    = false;
      }
      //else
        //LINFO("VAL %f",*iifeatureMaps);
      ++preCandidatePixelsIter; ++candidatePixelsIter; ++iifeatureMaps;
    }
  }
}

/*********************************************************************/

/*! find candidate pixels by just checking image pixels to see if they
   are within threshold of desired feature values
   if so, set candidatePixels image pixel to true
*/
template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIfindCandidatesNoBandPass()
{
  // FIX
  std::vector<int>::iterator imasterVec = SI_masterVec.begin();

  while(imasterVec != SI_masterVec.end())
    *imasterVec++ = -1;

  Image<bool>::iterator
    candidatePixelsIter    = SI_candidatePixels.beginw();

  // set initial candidates for each pixel

  while(candidatePixelsIter != SI_candidatePixels.endw())
  {
    *candidatePixelsIter = true;
    ++candidatePixelsIter;
  }


  typename std::vector<Image<FLOAT> >::iterator ifeatureMaps;
  typename std::vector<FLOAT>::iterator ilowThresh  = SI_lowThresh.begin();
  typename std::vector<FLOAT>::iterator ihighThresh = SI_highThresh.begin();
  // Run over each feature channel

  for(ifeatureMaps = SI_featureMaps->begin();
      ifeatureMaps !=  SI_featureMaps->end(); ++ifeatureMaps,
        ++ilowThresh, ++ihighThresh)
  {
    typename Image<FLOAT>::iterator iifeatureMaps = ifeatureMaps->beginw();
    candidatePixelsIter = SI_candidatePixels.beginw();

    // Run over each pixel

    while(iifeatureMaps != ifeatureMaps->endw())
    {
      //LINFO("Feature %f, high %f, low %f",*iifeatureMaps,*ihighThresh,*ilowThresh);
      if((*iifeatureMaps > *ihighThresh) || (*iifeatureMaps < *ilowThresh))
        *candidatePixelsIter = false;

      ++candidatePixelsIter; ++iifeatureMaps;
    }
  }
}

/*********************************************************************/
/*! if a pixel does not have an orthogonal neighbor then it is removed
   i.e. all pixels without 4 connectivity are removed as noise
*/
template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIremoveSingles()
{
  Image<bool>::iterator icandidatePixels = SI_candidatePixels.beginw();

  const int width  = SI_candidatePixels.getWidth();
  const int height = SI_candidatePixels.getHeight();
  for(int y = 0; y < height; y++)
  {
    for(int x = 0; x < width; x++)
    {
      if(*icandidatePixels)
      {
        unsigned int kill = 0;
        const int XLeft   = x - 1;
        const int XRight  = x + 1;
        const int YTop    = y - 1;
        const int YBottom = y + 1;
        if((XLeft >= 0) && (SI_candidatePixels.getVal(XLeft,y)))
          kill++;
        if((XRight < width)
           && (SI_candidatePixels.getVal(XRight,y)))
          kill++;
        if((YTop >= 0) && (SI_candidatePixels.getVal(x,YTop)))
          kill++;
        if((YBottom < height)
           && (SI_candidatePixels.getVal(x,YBottom)))
          kill++;
        if(kill < SI_killVal)
          *icandidatePixels = false;
      }
      ++icandidatePixels;
    }
  }
}

/*********************************************************************/
/*! if a pixel does not have an orthogonal neighbor then it is removed
   i.e. all pixels without 4 connectivity are removed as noise
*/
template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIremoveSinglesItr()
{
  Image<bool>::iterator icandidatePixels = SI_candidatePixels.beginw();

  const unsigned int width  = SI_candidatePixels.getWidth();
  const unsigned int height = SI_candidatePixels.getHeight();

  // hypothetically, this could cause a segmentation fault, but
  // we don't access this until the memory is aligned so we are OK
  Image<bool>::iterator icandidatePixelsTop    =
    SI_candidatePixels.beginw() - width;
  Image<bool>::iterator icandidatePixelsBottom =
    SI_candidatePixels.beginw() + width;
  Image<bool>::iterator icandidatePixelsLeft   =
    SI_candidatePixels.beginw() - 1;
  Image<bool>::iterator icandidatePixelsRight  =
    SI_candidatePixels.beginw() + 1;

  for(unsigned int y = 0; y < height; y++)
  {
    for(unsigned int x = 0; x < width; x++)
    {
      if(*icandidatePixels)
      {
        unsigned int kill = 0;
        if(x != 0)
          if(*icandidatePixelsLeft)
            kill++;
        if(x < width - 1)
          if(*icandidatePixelsRight)
            kill++;
        if(y != 0)
          if(*icandidatePixelsTop)
            kill++;
        if(y < height - 1)
          if(*icandidatePixelsBottom)
            kill++;
        if(kill < SI_killVal)
          *icandidatePixels = false;
      }

      ++icandidatePixels;
      ++icandidatePixelsBottom; ++icandidatePixelsRight;
      ++icandidatePixelsTop;    ++icandidatePixelsLeft;
    }
  }
}

/*********************************************************************/
/*! scan pixels from top to bottom/left to right. Link all pixels in the
    same blob. Do this using recursive back tracking to make a blob linked
    as one object and not a collection of objects
*/
template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIdiscreteLinking()
{
  const ushort width    = (ushort)SI_candidatePixels.getWidth();
  const ushort height   = (ushort)SI_candidatePixels.getHeight();
  bool trace            = false;
  INT pixID             = 0;
  SI_maxIDVal           = 1;
  SI_masters            = 1;
  SI_mastersCount       = 0;
  int lastNeighbor;

  for(ushort x = 0; x < width; x++)
  {
    trace        = false;
    lastNeighbor = -2;
    for(ushort y = 0; y < height; y++)
    {
      if(SI_candidatePixels.getVal(x,y))
      {
        if(!trace)
        {
          pixID++;
          trace        = true;
          lastNeighbor = -2;
        }
        SI_blobID.setVal(x,y,pixID);
        if(x > 0) // 1 or 0?
        {
          if(SI_candidatePixels.getVal(x-1,y))
          {
            // relink blob pixels if needed
            const INT check =  SI_blobID.getVal((x-1),y);
            if((signed)check != lastNeighbor)
            {
              SIbackwardLink(pixID,check);
              lastNeighbor = check;
            }
            else
            {
              SIbasicLink(pixID);
            }
          }
          else
          {
            lastNeighbor = -2;
            SIbasicLink(pixID);
          }
        }
        else
        {
          SIbasicLink(pixID);
        }
      }
      else
      {
        trace        = false;
        SI_blobID.setVal(x,y,0);
        lastNeighbor = -2;
      }
    }
  }
  LINFO("Candidate pixels recognized %d",pixID);
  SI_num = pixID;
}

/*********************************************************************/
/*! scan pixels from top to bottom/left to right. Link all pixels in the
    same blob. Do this using recursive back tracking to make a blob linked
    as one object and not a collection of objects
*/
template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIdiscreteLinkingOrtho()
{
  const int width       = SI_candidatePixels.getWidth();
  const int height      = SI_candidatePixels.getHeight();
  bool trace            = false;
  INT pixID             = 0;
  SI_maxIDVal           = 1;
  SI_masters            = 1;
  SI_mastersCount       = 0;
  int lastNeighbor;

  Image<bool>::iterator candidatePixelsItr = SI_candidatePixels.beginw();
  Image<int>::iterator  blobIDItr          = SI_blobID.beginw();

  for(int y = 0; y < height; y++)
  {
    trace        = false;
    lastNeighbor = -2;
    for(int x = 0; x < width; x++)
    {
      // This is: if(SI_candidatePixels.getVal(x,y))
      if(*candidatePixelsItr++)
      {
        // Start a new pixel row group
        if(!trace)
        {
          pixID++;
          trace        = true;
          lastNeighbor = -2;
        }
        // This is: SI_blobID.setVal(x,y,pixID);
        *blobIDItr++ = pixID;
        if(y > 0) // 1 or 0?
        {
          //This is: if(SI_candidatePixels.getVal(x,y-1))
          if(*(candidatePixelsItr - (1 + width)))
          {
            // relink blob pixels if needed
            //This is: const INT check =  SI_blobID.getVal(x,y-1);
            const INT check = *(blobIDItr - (1 + width));
            if((signed)check != lastNeighbor)
            {
              SIbackwardLink(pixID,check);
              lastNeighbor = check;
            }
            else
            {
              SIbasicLink(pixID);
            }
          }
          else
          {
            lastNeighbor = -2;
            SIbasicLink(pixID);
          }
        }
        else
        {
          SIbasicLink(pixID);
        }
      }
      else
      {
        trace        = false;
        //This is: SI_blobID.setVal(x,y,0);
        *blobIDItr++ = 0;
        lastNeighbor = -2;
      }
    }
  }
  //LINFO("Candidate pixels recognized %d",pixID);
  SI_num = pixID;
}

/*********************************************************************/
/*! relink pixels with new val. This allows pixels that have a new value
   to take on the old value of the blob. Slaves take on masters and
   masters can take on a slaves master.

*/
template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIbackwardLink(const INT slave,
                                                  const INT master)
{
  if(master > SI_maxIDVal)
    SI_maxIDVal = master;
  if(slave > SI_maxIDVal)
    SI_maxIDVal = slave;

  // my master has no master
  if(SI_masterVec[master] == -1)
  {
    // I have no master
    if(SI_masterVec[slave] == -1)
    {
      //LINFO("my master has no master/I have no master");
      SI_masterVec[master] = SI_masters;
      SI_masterVec[slave]  = SI_masters;
      SI_masters++;
      SI_mastersCount++;
    }
    // I already have a master
    else
    {
      //LINFO("my master has no master/I already have a master");
      SI_masterVec[master] = SI_masterVec[slave];
    }
  }
  // my master has a master
  else
  {
    // I have no master
    if(SI_masterVec[slave] == -1)
    {
      //LINFO("my master has a master/I have no master");
      SI_masterVec[slave] = SI_masterVec[master];
    }
    // I already have a master
    else
    {
      //LINFO("my master has a master/I already have a master");

      std::vector<int>::iterator masterVecItr = SI_masterVec.begin();
      for(INT i = 0; i <= SI_maxIDVal; i++, ++masterVecItr)
      {
        if(*masterVecItr == SI_masterVec[slave])
        {
          //LINFO("SET %d to %d",masterVec[i],masterVec[slave]);
          *masterVecItr = SI_masterVec[master];
        }
      }
      SI_masterVec[slave] = SI_masterVec[master];
      SI_mastersCount--;
    }
  }
}

/*********************************************************************/
template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIbasicLink(const INT node)
{
  if(node > SI_maxIDVal)
    SI_maxIDVal = node;

  // my master has no master
  if(SI_masterVec[node] == -1)
  {
    //LINFO("my master has no master/I have no master");
    SI_masterVec[node] = SI_masters;
    SI_masters++;
    SI_mastersCount++;
  }
}

/*********************************************************************/
/*! combine individual elements into a fully combined blob */
template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIcombine()
{
  SI_totalBlobs = 0;

  Image<bool>::iterator candidatePixelsItr = SI_candidatePixels.beginw();
  Image<int>::iterator  blobIDItr          = SI_blobID.beginw();

  while(candidatePixelsItr != SI_candidatePixels.endw())
  {
    if(*candidatePixelsItr++)
      *blobIDItr = SI_masterVec[*blobIDItr];
    blobIDItr++;
  }

  /*
  for(int x = 0; x < SI_blobID.getWidth(); x++)
  {
    for(int y = 0; y < SI_blobID.getHeight(); y++)
    {
      if(SI_candidatePixels.getVal(x,y))
        SI_blobID.setVal(x,y,SI_masterVec[SI_blobID.getVal(x,y)]);
    }
  }
  */

  std::vector<int>::iterator masterVecItr           = SI_masterVec.begin();
  typename std::vector<INT>::iterator reOrderVecItr = SI_reOrderVec.begin();
  std::vector<bool>::iterator resetItr              = SI_reset.begin();

  for(INT x = 0; x < SI_num; x++, ++masterVecItr)
  {
    bool add = true;
    //LINFO("Master Vec %d is %d",x,masterVec[x]);
    for(INT y = 0; y < SI_totalBlobs; y++)
    {
      if((int)SI_reOrderVec[y] == *masterVecItr)
        add = false;
    }

    if((add) && (*masterVecItr != -1))
    {
      //LINFO("DOING %d",masterVec[x]);
      *reOrderVecItr = *masterVecItr;
      SI_reverseOrderVec[*masterVecItr] = SI_totalBlobs;
      *resetItr = true;
      SI_totalBlobs++; ++reOrderVecItr; ++resetItr;
    }
  }
}

/*********************************************************************/

template SI_TEMPLATE_CLASS inline
void segmentImageMC2<SI_TEMPLATE>::SIdoSegment()
{
  ASSERT(SI_set1); ASSERT(SI_set2); ASSERT(SI_set3); ASSERT(SI_set4);
  //LINFO("FINDING CANDIDATES");
  if(SI_useCandidateBandPass == true)
    SIfindCandidates();
  else
    SIfindCandidatesNoBandPass();
  //
  if(SI_removeSingles == true)
  {
    //LINFO("REMOVING SINGLES");
    SIremoveSinglesItr();
  }
  //LINFO("LINKING PIXELS");
  //SIdiscreteLinking();
  SIdiscreteLinkingOrtho();
  //LINFO("RELABELING");
  SIcombine();
  //LINFO("DONE");
}

#undef SI_TEMPLATE_CLASS
#undef SI_TEMPLATE

template class segmentImageMC2<float, unsigned int, 1>;
template class segmentImageMC2<float, unsigned int, 2>;
template class segmentImageMC2<float, unsigned int, 3>;
template class segmentImageMC2<float, unsigned int, 4>;

// ######################################################################
/* So things look consistent in everyone's emacs... */
/* Local Variables: */
/* indent-tabs-mode: nil */
/* End: */

#endif