MixtureOfGaussians.cpp

/*
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
Multimediale Systeme der Informationsverarbeitung
Institut fuer Informatik
Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS 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.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/



#include "Image/MixtureOfGaussians.hh"
#include <Base/Image/ImageConvert.hh>
#include <Base/ImageUtils/ImageDraw.hh>
#include <Base/Image/ROI.hh>
#include <bias_config.h>
#include <Base/Image/ImageIO.hh>
#include <Base/Common/FileHandling.hh>
#include <Filter/Erosion.hh>
#include <Base/Debug/TimeMeasure.hh>

using namespace BIAS;
using namespace std;
// public

template <class StorageType>
MixtureOfGaussians<StorageType>::
MixtureOfGaussians(const unsigned int numDistributions,
                   const unsigned int sizeWindow,
                   const float updateFactor,
                   const MoGColorModel colorModel) {
  dDistributions_ = numDistributions;
  dWindow_ = sizeWindow /2;
  bInited_ = false;
  fUpdateFactor_ = updateFactor;
  dMaxWeight_ = 50;
  dMinWeight_ = 20;
  fError_ = 2.9f * 2.9f;
  dFrameCounter_ = 0;
  bSmooth_ = false;
  colorModel_ = colorModel;
  blobDetector_.SetMinSize(1000);
  bAllowNewDistributionCreation_ = true;
  bCleanUp_ = false;
  bInverseCleanup_ = false;  
  blobDetectorRan_ = false;
  //we are working on z
  maxDepth_ = 1.0;
#ifdef BIAS_HAVE_PTHREADS
  guardimagemutex = new pthread_mutex_t;
  pthread_mutex_init(guardimagemutex, NULL);
#endif
}

template <class StorageType>
void MixtureOfGaussians<StorageType>::
SetColorModel(MoGColorModel model){
     colorModel_ = model;
}

template <class StorageType>
void MixtureOfGaussians<StorageType>::
SetSmoothing(bool smooth) {
  bSmooth_ = smooth;
  if (smooth) 
    blobDetector_.SetGaussSigma(1.0);
  else 
    blobDetector_.SetGaussSigma(0.0);
}

template <class StorageType>
int MixtureOfGaussians<StorageType>::
InitWeights(const int maxWeight,
            const int minWeight,
            const float errorMargin) {
  dMaxWeight_ = maxWeight;
  dMinWeight_ = minWeight;
  fError_ = errorMargin * errorMargin;
  return 0;
}

template <class StorageType>
MixtureOfGaussians<StorageType>::
~MixtureOfGaussians() {
#ifdef BIAS_HAVE_PTHREADS
  delete guardimagemutex;
#endif
}

template <class StorageType>
int MixtureOfGaussians<StorageType>::
Apply(BIAS::Image<StorageType>& in) {

    TimeMeasure timer;
    
    blobDetectorRan_ = false;
  if (!bInited_) {
    channels_ = in.GetChannelCount();
    width_ = in.GetWidth();
    height_ = in.GetHeight();
  }
  in.GetROI()->GetCorners(roiULx_, roiULy_, roiLRx_, roiLRy_);
  blobDetector_.SetROI(BIAS::Vector2<unsigned>(roiULx_,roiULy_),
                       BIAS::Vector2<unsigned>(roiLRx_,roiLRy_));
  if( bSmooth_ ){
    BIAS::Gauss<StorageType,StorageType> gs;   
    gs.SetSigma(2.0);
    gs.SetHalfWinSize(3);
    //gs.SetSize(5);
    gs.Filter(in,in);
  }
 
  //create a normalized image with float
  CreateNormalizedImage_(in, workImage_);  
  
#ifdef BIAS_HAVE_PTHREADS
  pthread_mutex_lock(guardimagemutex);
#endif
  if (!bInited_) {
    timer.Start();
    Init_();
    timer.Stop();
    //cout<<" MoG Init_ took:"; timer.PrintRealTime();
    timer.Reset();
  }
  else {
    timer.Start();
    Update_();
    timer.Stop();
    //cout<<" MoG Update_ took:"; timer.PrintRealTime();
    timer.Reset();
  }
  
  
    timer.Start();
  if (dFrameCounter_ > dMinWeight_) {
    //inverse cleanup, delete pixel 
    if (bInverseCleanup_) {
      Image<unsigned char> white;
      white.Init(differenceSave_.GetWidth(), differenceSave_.GetHeight(), differenceSave_.GetChannelCount());
      white.FillImageWithConstValue((unsigned char)255);
      differenceSave_ = white - differenceSave_;
      blobDetector_.SetEraseBelowMinSize(true);
      if(!blobDetectorRan_){
        blobDetector_.Detect(differenceSave_, corners_);
        blobDetectorRan_ = true;
      }
      differenceSave_ = white - differenceSave_;     
    }
    
    //clean up pixel below specified minimum size 
    if (bCleanUp_) {
      blobDetector_.SetEraseBelowMinSize(true);
      if(!blobDetectorRan_){
        blobDetector_.Detect(differenceSave_, corners_);
        blobDetectorRan_ = true;
      }
    }
    
    //erode 3 pixel if requested
    if (bErosion_) {
      Erosion<unsigned char, unsigned char> erosion;
      erosion.Erode3Fast(differenceSave_, difference_);
      erosion.Erode3Fast(difference_, differenceSave_);
    }
  }
  
  
#ifdef BIAS_HAVE_PTHREADS
  pthread_mutex_unlock(guardimagemutex);
#endif
  dFrameCounter_++;
#if 0
  if (dFrameCounter_ > 250 && dFrameCounter_ < 300) {
  if (in.GetStorageType() == ImageBase::ST_unsignedchar) {
    ImageIO::Save("mog-normalized-"+FileHandling::LeadingZeroString(dFrameCounter_,3)+".mip", workImage_);
    ImageIO::Save("mog-weight-"+FileHandling::LeadingZeroString(dFrameCounter_,3)+".mip", weight_);
    ImageIO::Save("mog-cov-"+FileHandling::LeadingZeroString(dFrameCounter_,3)+".mip", cov_);
    ImageIO::Save("mog-guess-"+FileHandling::LeadingZeroString(dFrameCounter_,3)+".mip", guess_);
  }
  }
#endif
  timer.Stop();
  //cout<<" MoG Cleanup_ took:"; timer.PrintRealTime();
  timer.Reset();
    
  return 0;
}

template <class StorageType>
int MixtureOfGaussians<StorageType>::
GetNormalizedImage(BIAS::Image<float> &normalized) {    
  normalized = workImage_;  
  return 0;
}

template <class StorageType>
int MixtureOfGaussians<StorageType>::
GetWeightImage(BIAS::Image<float> &weight) {
  weight = weight_;
  return 0;
}

template <class StorageType>
int MixtureOfGaussians<StorageType>::
GetMaxWeightImage(BIAS::Image<float> &weight) {
  if (weight.IsEmpty()) {
    weight.Init(width_,height_,channels_,ImageBase::ST_float,false);
  }
  weight.FillImageWithConstValue(0.0);
  float **idaRet = weight.GetImageDataArray();
  float **idaWeight = weight_.GetImageDataArray();
  float val;
  for (unsigned int y=0;y<height_;y++) {
    for (unsigned int x=0;x<width_*channels_;x++) {
      for (unsigned int i=0;i<dCurDistributions_;i++) {
        val = idaWeight[y][x*dDistributions_+i];
        if (val > idaRet[y][x])
          idaRet[y][x] = val; //assign biggest weight of all distributions to weight output image
      }
    }
  }
  return 0;
}

template <class StorageType>
int MixtureOfGaussians<StorageType>::
GetMatchValueImage(BIAS::Image<float> &match) {
#ifdef BIAS_HAVE_PTHREADS
  pthread_mutex_lock(guardimagemutex);
#endif
  match = matchValueImage_;
#ifdef BIAS_HAVE_PTHREADS
  pthread_mutex_unlock(guardimagemutex);
#endif

  return 0;
}

template <class StorageType>
int MixtureOfGaussians<StorageType>::
GetDifferenceImage(BIAS::Image<unsigned char> &diff) {
  // BIAS::Image<float> tmp(cov_);
  // tmp.ScaleShiftBetween(0,255);
#ifdef BIAS_HAVE_PTHREADS
  pthread_mutex_lock(guardimagemutex);
#endif
  diff = differenceSave_;
#ifdef BIAS_HAVE_PTHREADS
  pthread_mutex_unlock(guardimagemutex);
#endif
  if (dFrameCounter_ < dMinWeight_)
    return -1;
  else
    return 0;
}

template <class StorageType>
int MixtureOfGaussians<StorageType>::
GetDifferenceImageWithVisuals(BIAS::Image<unsigned char> &guiImage) {
  if (dFrameCounter_ < dMinWeight_)
    return -1;
  unsigned char **visu = differenceVisu_.GetImageDataArray();
  
#ifdef BIAS_HAVE_PTHREADS
  pthread_mutex_lock(guardimagemutex);
#endif
  unsigned char **idaD = differenceSave_.GetImageDataArray();
  for (unsigned int y=0;y<height_;y++) {
    for (unsigned int x=0;x<width_;x++) {
      for (unsigned int c=0;c<3;c++) {
        visu[y][x*3+c] = idaD[y][x];
      }
    }
  }
#ifdef BIAS_HAVE_PTHREADS
  pthread_mutex_unlock(guardimagemutex);
#endif
  if(!blobDetectorRan_){
    blobDetector_.Detect(differenceSave_, corners_);
    blobDetectorRan_ = true;
  }  
  blobDetector_.DrawInImage(differenceVisu_);
  
  guiImage = differenceVisu_;
  return 0;
}

template <class StorageType>
int MixtureOfGaussians<StorageType>::
GetBoundingBoxes(std::vector<BIAS::BIASBlob> &corners) {
  if (corners_.size()>0) {
    corners = corners_;
    return 0;
  }
  else
    return -1;
}

template <class StorageType>
int MixtureOfGaussians<StorageType>::
GetLargestBoundingBox(BIAS::BIASBlob &corners) {
    if(!blobDetectorRan_){
        blobDetector_.Detect(differenceSave_, corners_);
        blobDetectorRan_ = true;
    }
  return blobDetector_.GetLargestBoundingBox(corners);
}

template <class StorageType>
int MixtureOfGaussians<StorageType>::
GetCurrentDistributions(std::vector<Image<float> > &distributionImages) {
  distributionImages.clear();
  distributionImages.resize(dCurDistributions_);

  float **idaGuess = guess_.GetImageDataArray();

  for (unsigned int i=0; i<dCurDistributions_;i++) {
    if (!distributionImages[i].IsEmpty())
      distributionImages[i].Release();
    distributionImages[i].Init(width_, height_, channels_);
    float **idaD = distributionImages[i].GetImageDataArray();
    for (unsigned int y=0; y<height_; y++) {
      for (unsigned int x=0; x<width_; x++) {
        for (unsigned int c=0; c<channels_; c++) {
          idaD[y][x*channels_+c] = idaGuess[y][(x*channels_+c)*dDistributions_+i];
          //idaD[y][x*channels_+c] = idaCov[y][(x*channels_+c)*dDistributions_+i];
          //idaD[y][x*channels_+c] = idaWeight[y][(x*channels_+c)*dDistributions_+i];
        }
      }
    }
  }
  return 0;
}

//  Protected

template <class StorageType>
int MixtureOfGaussians<StorageType>::
CreateNormalizedImage_(const BIAS::Image<StorageType> &in, BIAS::Image<float> &out) {
    height_ = in.GetHeight();
    width_= in.GetWidth();
    
    if (in.GetStorageType() == ImageBase::ST_float) {
        out = in;
    } else {
        ImageConvert::ConvertST(in, out, ImageBase::ST_float);
    }
    
    float **idaFL = out.GetImageDataArray();
    channels_ = in.GetChannelCount();
    unsigned int depthChannel =0;
    switch(colorModel_){
    //first normal RGB
    case MOG_RGB:
    case MOG_RGBD:
        depthChannel = 3;
#ifdef BIAS_HAVE_OPENMP
#pragma omp parallel for
#endif
        for (int y=0;y<(int)height_;y++) {
            for (int x=0;x<(int)width_;x++) {
                for (unsigned c=0;c<channels_;c++) {
                    if(c==depthChannel) //depth channel
                        idaFL[y][x*channels_+c] /= maxDepth_;
                    else
                        idaFL[y][x*channels_+c] /= 255.0f;
                }
            }
        }
    break;

    case MOG_HS:    
    case MOG_HSD:{  
        depthChannel = 2;
        if(colorModel_== MOG_HS ) channels_ = 2; //HS 
        else if(colorModel_== MOG_HSD ) channels_ = 3; //HSD

#ifdef BIAS_HAVE_OPENMP
#pragma omp parallel for
#endif
        for (int y=0;y<(int)height_;y++) {
            for (int x=0;x<(int)width_;x++) {
                for (unsigned c=0;c<channels_;c++) {
                    if(c==depthChannel) //depth channel
                        idaFL[y][x*channels_+c] /= maxDepth_;
                    else
                        idaFL[y][x*channels_+c] /= 255.0f;
                }
            }
        }       
    }
    break;
    case MOG_NRGB:
    case MOG_NRGBD:{
        depthChannel = 3;
#ifdef BIAS_HAVE_OPENMP
#pragma omp parallel for 
#endif
        for (int y=0;y<(int)height_;y++) {
            for (int x=0;x<(int)width_;x++) {
                //calculate normalized 
                float norm = 0;
                for(int c =0;c<(int)channels_;c++){
                    if(c<3) norm += (idaFL[y][x*channels_+c])*(idaFL[y][x*channels_+c]);
                }
                norm = (float) sqrt(norm);
                // Safeguard, we do not want to divide by something < 5 (why not??)
                //if (norm < 5.0f) norm = 5.0f;
                if (norm <= 1.0f) norm = 1.0f;
                for(unsigned c =0;c<channels_;c++) {
                    if(c==depthChannel)
                        idaFL[y][x*channels_+c] /= maxDepth_;
                    else
                        idaFL[y][x*channels_+c] /= norm;
                }
            }
        }
    }
    break;
    default:
        break;
    }

    return 0;
}

template <class StorageType>
int MixtureOfGaussians<StorageType>::
Init_() {
  dCurDistributions_ = 1;
  guess_.ReInit(width_*channels_,height_,dDistributions_,
                ImageBase::ST_float,false, ImageBase::CM_Grey);
  cov_.ReInit(width_*channels_,height_,dDistributions_,
              ImageBase::ST_float,false, ImageBase::CM_Grey);
  weight_.ReInit(width_*channels_,height_,dDistributions_,
                 ImageBase::ST_float,false, ImageBase::CM_Grey);
  
  maxWeight_.ReInit(width_,height_,channels_,ImageBase::ST_float);
  
  matchValueImage_.ReInit(width_,height_,1,ImageBase::ST_float,false, ImageBase::CM_Grey);
  matchValueImage_.FillImageWithConstValue(0.0f);
  difference_.ReInit(width_,height_,1,ImageBase::ST_unsignedchar,false, ImageBase::CM_Grey);
  differenceSave_.ReInit(width_,height_,1,ImageBase::ST_unsignedchar,false, ImageBase::CM_Grey);
  differenceSmooth_.ReInit(width_,height_,1,ImageBase::ST_unsignedchar,false, ImageBase::CM_Grey);
  differenceVisu_.ReInit(width_,height_,3,ImageBase::ST_unsignedchar,false, ImageBase::CM_RGB);
  

  float **idaImage = workImage_.GetImageDataArray();
  float **idaWeight = weight_.GetImageDataArray();
  float **idaGuess= guess_.GetImageDataArray();
  float **idaCov = cov_.GetImageDataArray();;

    //fill guesses and covariances
  for (unsigned int c=0;c<channels_;c++) {
    for (unsigned int y=0;y<height_;y++) {
      for (unsigned int x=0;x<width_;x++) {
         float valGuess = 0;
         int cntPix = 0;
         for (unsigned int yw=y-dWindow_;yw<y+dWindow_;yw++) {
           for (unsigned int xw=x-dWindow_;xw<x+dWindow_;xw++) {
             if (yw<0 || xw<0 || yw>=height_ || xw >= width_)
               continue;
             valGuess += idaImage[yw][xw*channels_+c];
                      
             cntPix++;
           }
         }
         if (cntPix != 0){      
             valGuess /= (float)cntPix;
         }
         else{
             valGuess = 0;
         }
         idaGuess[y][(x*channels_+c)*dDistributions_] = valGuess;       
      }
    }
  }
  //fill covariances
  int xc = 0;
  for (unsigned int c=0;c<channels_;c++) {
    for (unsigned int y=0;y<height_;y++) {
      for (unsigned int x=0;x<width_;x++) {
        float val = 0;
        int cntPix = 0;
        for (unsigned int yw=y-dWindow_;yw<y+dWindow_;yw++) {
          for (unsigned int xw=x-dWindow_;xw<x+dWindow_;xw++) {
            if (yw<0 || xw<0 || yw>=height_ || xw >= width_)
              continue;
            val += (idaImage[yw][xw*channels_+c] - idaGuess[yw][(xw*channels_+c)*dDistributions_])
                * (idaImage[yw][xw*channels_+c] - idaGuess[yw][(xw*channels_+c)*dDistributions_]);
            cntPix++;
          }
        }
        if (cntPix!=0)
          val /= (float)(cntPix*cntPix);
        else
          val = 1;
        idaCov[y][(x*channels_+c)*dDistributions_] = val;
      }
    }
  }
  
  //reset Guess, Cov and Weight
  for (unsigned int y=0;y<height_;y++) {
    for (unsigned int x=0;x<width_*channels_;x++) {
      for (unsigned int i=1;i<dDistributions_;i++) {
        xc = x*dDistributions_+i;
        idaWeight[y][xc] = -1;
        idaGuess[y][xc] = -1;
        idaCov[y][xc] = -1;
      }
    }
  }
  bInited_ = true;
  return 0;
}

template <class StorageType>
int MixtureOfGaussians<StorageType>::
InitOnePixel_(unsigned int x,unsigned int y,unsigned int c,unsigned int distrib) {
  float val = 0;
  int xc = 0;
  int cntPix = 0;
  float **idaImage = workImage_.GetImageDataArray();
  float **idaGuess = guess_.GetImageDataArray();
  float **idaCov = cov_.GetImageDataArray();
  float **idaWeight = weight_.GetImageDataArray();

  for (unsigned int yw=y-dWindow_;yw<y+dWindow_;yw++) {
    for (unsigned int xw=x-dWindow_;xw<x+dWindow_;xw++) {
      if (yw<0 || xw<0 || yw>=height_ || xw >= width_)
        continue;
      val += idaImage[yw][xw*channels_+c];
      cntPix++;
    }
  }
  if (cntPix!=0)
    val /= (float)cntPix;
  else
    val = 0;
  xc = (x*channels_+c)*dDistributions_+distrib;

//  cout << idaGuess[y][xc] << " " << idaCov[y][xc];

  idaGuess[y][xc] = val;

  val = 0;
  cntPix = 0;
  for (unsigned int yw=y-dWindow_;yw<y+dWindow_;yw++) {
    for (unsigned int xw=x-dWindow_;xw<x+dWindow_;xw++) {
      if (yw<0 || xw<0 || yw>=height_ || xw >= width_)
        continue;
      xc = (xw*channels_+c)*dDistributions_+distrib;
      val += (idaImage[yw][xw*channels_+c] - idaGuess[yw][xc]) *
      (idaImage[yw][xw*channels_+c] - idaGuess[yw][xc]);
      cntPix++;
    }
  }
  if (cntPix!=0)
    val /= (float)(cntPix*cntPix);
  else
    val = 1;
  xc = (x*channels_+c)*dDistributions_+distrib;
  idaCov[y][xc] = val;
  idaWeight[y][xc] = 1;
//  cout << "  " << idaGuess[y][xc] << " " << idaCov[y][xc] << endl;
  return 0;
}

template <class StorageType>
int MixtureOfGaussians<StorageType>::
Update_() {
  // fill difference image white (matched pixel will be set to black,
  // so changes are white in the result)
  difference_.FillImageWithConstValue((unsigned char)255);

  float **idaImage = workImage_.GetImageDataArray();
  float **idaWeight = weight_.GetImageDataArray();
  float **idaGuess= guess_.GetImageDataArray();
  float **idaCov = cov_.GetImageDataArray();;
  float **idaMaxWeight = maxWeight_.GetImageDataArray();;
  float **idaMatchValue = matchValueImage_.GetImageDataArray();
  unsigned char **idaDiff  = difference_.GetImageDataArray();

  // iterate through whole picture
  //int x=0,y=0,c=0,i=0;
  for (int y=(int)roiULy_ ; y<(int)roiLRy_ ; y++) {
    int xc = 0; bool matched=false;
//#ifdef BIAS_HAVE_OPENMP
//#pragma omp parallel for //private(y,x,c,i)
//#endif
    for (int x=(int)roiULx_;x<(int)roiLRx_;x++) {
      unsigned int timesMatched_=0;
      float matchValue_=0.0f;
      unsigned distCount =0; 
      //////////////////CHANNELS/////////////////////////////
      for (int c=0;c<(int)channels_;c++) {
        // initialized with large number
        float min_weight = 10000;
        // initialized with small number
        float max_weight = -10000;
        // worst distribution (i.e. matched minimal times)
        unsigned int min_distrib = 0;
        // best distribution (i.e. matched maximal times)
        //unsigned int max_distrib = 0;

        // find min/max distribs
        for (int i=0;i<(int)dDistributions_;i++) {
          xc = (x*channels_+c)*dDistributions_+i;
          // find distribution with maximal weight
          if (idaWeight[y][xc] > max_weight) {
            max_weight = idaWeight[y][xc];
            //max_distrib = i;
          }
          // find distribution with minimal weight
          if (idaWeight[y][xc] < min_weight) {
            min_weight = idaWeight[y][xc];
            min_distrib = i;
          }
          if (idaWeight[y][xc] < 0)
            break;
        }

        matched = false;
        // value in orig image
        float val = idaImage[y][x*channels_+c];
        ///////////////////DISTRIBUTIONS///////////////////////////////
        // check for match in all distributions       
        for (int i=0;i<(int)dCurDistributions_;i++) {
          xc = (x*channels_+c)*dDistributions_+i;
          // guess and covariance in this distribution          
          float guess = idaGuess[y][xc];
          if (guess < 0) continue;
          float cov = idaCov[y][xc];          
          float absdiff = abs(val - guess);
          matchValue_ += absdiff;
          distCount++;
          float absdiffSquared = absdiff*absdiff;
          ////////////////////MATCH////////////////////////////
          if (!matched &&  absdiffSquared < fError_*cov) {
            // match! weight is increased
            if (idaWeight[y][xc] <= dMaxWeight_)
              idaWeight[y][xc]++;
            // guess and covariance get updated
            float myCovTest=0.0f,myGuessTest =0.0f;
            
            myGuessTest = (1.0f-fUpdateFactor_)*guess + fUpdateFactor_*val;
            myCovTest = (1.0f-fUpdateFactor_)*cov + fUpdateFactor_*absdiffSquared;
            
            if (myCovTest < 1e-5) {
                idaCov[y][xc] = 1e-5f;
            } else if (myCovTest > 0.05) {
                idaCov[y][xc] = 0.05f;
            }
            
            idaGuess[y][xc] = myGuessTest;
           
            timesMatched_++;
            matched = true;
          }// end if (!matched && absdiff*absdiff < fError_*cov) 
          ///////////////////END MATCH /////////////////////////////
          ///////////////////NO MATCH /////////////////////////////
          else { // no match! decrease weight                 
            if (bAllowNewDistributionCreation_) {
              if (idaWeight[y][xc] > dMinWeight_)
              idaWeight[y][xc]--;
            }
          }//end else 
        } //end going over distributions
        /////////////////END DISTRIBUTIONS////////////////////////////
        // if no match is found, overwrite least matched distribution
        //with newly initialized values
        if (!matched && bAllowNewDistributionCreation_) {
          InitOnePixel_(x,y,c,min_distrib);
          if (min_distrib >= dCurDistributions_) {
            dCurDistributions_ = min_distrib +1;
          }
        }
        // update max weight image (only for vis)
        idaMaxWeight[y][x] = max_weight;
      } //channels
      ///////////////////END CHANNELS////////////////////////
      //if matched on all channels
      if (timesMatched_ == channels_) {
        idaDiff[y][x] = 0; //diff is black
      } else {
        idaDiff[y][x] = 255; //else diff white
      }
      
      if (matchValue_ <= 0) {
        idaMatchValue[y][x] = 0;
      } else {
        idaMatchValue[y][x] = matchValue_ / (float)(distCount);
      }
    }//x
  }//y
  differenceSave_.CopyIn_NoInit(difference_.GetImageData());
  return 0;
}

namespace BIAS{
template class MixtureOfGaussians<unsigned char>;
template class MixtureOfGaussians<float>;
}