Label.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
*/

#ifdef WIN32
# pragma warning (disable: 4251) // STL MAP DLL warning noise (JW)
#endif

#include <limits.h>
#include <stdio.h>
#include "Label.hh"

#include <Base/Image/ImageIO.hh>
#include <Base/Image/ImageConvert.hh>

#include <algorithm>

#include <Base/Common/BIASpragma.hh>
#include <MathAlgo/Median1D.hh>

#define ERR_ILLT 0 /* Datenstruktur wird nicht unterstuetzt */
#define ERR_ILLI 1 /* Datentyp wird nicht unterstuetzt */
#define ERR_NOLB 0 /*  Zu viele Gebiete */

#define MAX_SEG UCHAR_MAX
#define MAX_LAB USHRT_MAX

#define SEGDEF 1

using namespace BIAS;
using namespace std;



Label::Label() 
{
  LabelInitialized_=false;
  _FilterBorderHandling=TBH_same;
  //AddDebugLevel(D_LABEL);
}

Label::~Label()
{}


void Label::LabelInit_(const ImageBase& Image)
{
  BIASDOUT(D_INIT,"LabelInit_()");
  unsigned int UpperLeftX, UpperLeftY, LowerRightX, LowerRightY;
  Image.GetROI()->GetCorners(UpperLeftX, UpperLeftY, LowerRightX, LowerRightY);
  LabelSumImage_.Init(LowerRightX - UpperLeftX, LowerRightY - UpperLeftY);
  LabelResultImage_.Init(LowerRightX - UpperLeftX, 
                         LowerRightY - UpperLeftY);
  LabelInitialized_=true;
}
/*
void Label::LabelInit_(const Image<float>& Image)
{
  BIASDOUT(D_INIT,"LabelInit_()");
  unsigned int UpperLeftX, UpperLeftY, LowerRightX, LowerRightY;
  Image.GetROI()->GetCorners(UpperLeftX, UpperLeftY, LowerRightX, LowerRightY);
  LabelSumImage_.Init(LowerRightX - UpperLeftX, LowerRightY - UpperLeftY);
  LabelResultImage_.Init(LowerRightX - UpperLeftX, 
                         LowerRightY - UpperLeftY);
  LabelInitialized_=true;
}

void Label::LabelInit_(const Image<unsigned char>& Image)
{
  BIASDOUT(D_INIT,"LabelInit_()");
  unsigned int UpperLeftX, UpperLeftY, LowerRightX, LowerRightY;
  Image.GetROI()->GetCorners(UpperLeftX, UpperLeftY, LowerRightX, LowerRightY);
  LabelSumImage_.Init(LowerRightX - UpperLeftX, LowerRightY - UpperLeftY);
  LabelResultImage_.Init(LowerRightX - UpperLeftX, 
                         LowerRightY - UpperLeftY);
  LabelInitialized_=true;
}
*/

void Label::LabelRelease_()
{
    BIASDOUT(D_INIT,"LabelShortRelease_()");
    LabelSumImage_.Release();
    LabelResultImage_.Release();
    LabelInitialized_=false;
}



int Label::Label4Neighbour_()
{
  BIASCDOUT(D_FILTERBASE_CALLSTACK, "Label::Label4Neighbour_\n");
  int nx, ny, i, k, no_seg, lab_seg2, lab_segv, no_all, *seg_list = NULL;
  register int value, lab_seg1;
  register LABEL_CALC_TYPE *adr, *adr_out;
  // initialize
  try { seg_list = new int[MAX_LAB]; }
  catch (bad_alloc) { BIASERR("new failed "); }
  nx = LabelSumImage_.GetWidth();
  ny = LabelSumImage_.GetHeight();
  adr = LabelSumImage_.GetImageData();
  adr_out = LabelResultImage_.GetImageData();
  no_seg = 0;
    
  /*****
   *
   *     FIRST LINE
   *
   *****/
  adr++;
  *adr_out++ = 0;
  seg_list[0] = SEGDEF;
  lab_seg1 = 0;
  for ( k = nx-1 ; k > 0 ; k--)
    {
      if (*adr != *(adr-1))
        {
          no_seg++;
          seg_list[no_seg] = SEGDEF;
          lab_seg1 = no_seg;
        };
      adr++;
      *adr_out++ = lab_seg1;
    };

  /*****
   *
   *     NEXT LINES
   *
   *****/
  for ( i = ny-1 ; i > 0 ; i--)
    {

      /*****
       *
       *     FIRST PIXEL
       *
       *****/
      value = *adr;
      lab_seg1 = -1;
      if ( *(adr-nx) == value ) lab_seg1 = *(adr_out-nx);

      if ( lab_seg1 == -1 )
        {
          no_seg++;
          seg_list[no_seg] = SEGDEF;
          lab_seg1 = no_seg;
        };
      adr++;
      adr_out++;

      /*****
       *
       *     CHECK LABEL OF NEXT PIXELS
       *
       *****/
      lab_segv = lab_seg1;
      for ( k = nx-1 ; k > 0; k--)
        {
          lab_seg1 = -1;
          lab_seg2 = *(adr_out-nx); 
          value = *adr;
          if ( *(adr-1) == value)
            {
              lab_seg1 = lab_segv;
            }
          else if ( *(adr-nx) == value)
            {
              lab_seg1 = *(adr_out-nx);
            };
          /*****
           *
           *     NEW REGION ?
           *
           *****/
          if ( lab_seg1 == -1)
            {
              no_seg++;
              seg_list[no_seg] = SEGDEF;
              lab_seg1 = no_seg;
              if (no_seg > MAX_LAB) {
                BIASERR("error: too many blobs " << ERR_NOLB);
              }
            };

          /*****
           *
           *     ACTUALIZE REGION LINE
           *
           *****/
          *(adr_out-1) = lab_segv;
          lab_segv = lab_seg1;

          /*****
           *
           *     CHECK REGION MERGING
           *
           *****/
          if ( k != 1 && value == *(adr-nx) && lab_seg1 != lab_seg2)
            {
              int m;
              seg_list[lab_seg2] = -lab_seg1-1;
              for ( m = nx ; m > 0; m--)
                {
                  if ( *(adr_out-m) == lab_seg2)
                    *(adr_out-m) = lab_seg1;
                };
              lab_segv = lab_seg1;
            };

          adr_out++;
          adr++;
        };
      *(adr_out-1) = lab_segv;
    };

  /*****
   *
   *     GET NEW LABEL
   *
   *****/
  no_all = 0;
  for ( i = 0; i <= no_seg ; i++)
    {
      if (seg_list[i] == SEGDEF)
        {
          seg_list[i] = no_all;
          no_all++;
        };
    };

  /*****
   *
   *     LABEL MERGED REGIONS
   *
   *****/
  for ( i = 0; i <= no_seg ; i++)
    {
      int ind;
      ind = i;
      while (seg_list[ind] < 0 )
        {
          ind = -seg_list[ind]-1;
        };
      seg_list[i] = seg_list[ind];
    };

  /*****
   *
   *     NEW IMAGE-LABEL
   *
   *****/
  adr_out = LabelResultImage_.GetImageData();
  for ( i = nx*ny ; i > 0 ; i--)
    {
      *adr_out = seg_list[*adr_out];
      adr_out++;
    };

  delete [] seg_list;
  return (no_all);
}



int Label::Label4Neighbour(const BIAS::Image<LABEL_CALC_TYPE>& Image, 
                           BIAS::Image<LABEL_CALC_TYPE>& DestinationImage)
{
  BIASCDOUT(D_FILTERBASE_CALLSTACK, "Label::Label4Neighbour\n");
  int CornerCount=0;
  if (!LabelInitialized_){ 
    LabelInit_(Image);
  }
  if ( ! Image.SamePixelAndChannelCount(LabelSumImage_) ){
    LabelRelease_();
    LabelInit_(Image); 
  }
  Image.GetCopyOfROI(LabelSumImage_); 
    
  CornerCount = Label4Neighbour_();
  
  DestinationImage.Paste2ROI(LabelResultImage_);

  return CornerCount;
}


int Label::Label4Neighbour(const Image<unsigned char>& image, 
                           Image<unsigned char>& DestinationImage)
{
  BIASCDOUT(D_FILTERBASE_CALLSTACK, "Label::Label4Neighbour\n");
  int CornerCount=0;

  if (!LabelInitialized_) 
    LabelInit_(image);

 
  if ( ! image.SamePixelAndChannelCount(LabelSumImage_) ){
    LabelRelease_();
    LabelInit_(image); 
  }

  // explicit cast
  Image<unsigned char> roi;
  image.GetCopyOfROI(roi);
  ImageConvert::ConvertST(roi, LabelSumImage_, LABEL_STORAGE_TYPE);

  //      ImageIO::Save("labelsumimage.",LabelSumImage_);
    

  CornerCount = Label4Neighbour_();
    
  if (CornerCount > MAX_SEG){
    BIASERR("Label::Label4Neighbour() Segmentzahl "
            << CornerCount <<" der Label zu gross\n");
    CornerCount = MAX_SEG;
  }

  DestinationImage.Paste2ROI(LabelResultImage_);

  return CornerCount;
}

int Label::Process(const Image<unsigned char>& im)
{
  BIASCDOUT(D_FILTERBASE_CALLSTACK, "Label::Process\n");
  if  (!LabelInitialized_) 
    LabelInit_(im);
  
  Image<unsigned char> roiImg;
  im.GetCopyOfROI(roiImg);
  unsigned minx, miny, maxx, maxy;
  im.GetROI()->GetCorners(minx, miny, maxx, maxy);
  /*static int coun=0;
    ostringstream os("");
    os << setw(4)<<setfill('0')<<coun;
    ImageIO::Save("roi-"+os.str()+".mip", roiImg);*/
  if (!LabelResultImage_.SamePixelAndChannelCount(roiImg) ) {
    LabelRelease_();
    LabelInit_(roiImg);
  }
  ImageConvert::ConvertST( roiImg, LabelSumImage_ ,LABEL_STORAGE_TYPE);
  //  ImageIO::Save("lsi-"+os.str()+".mip", LabelSumImage_);
  
  int width, height;
  width=LabelResultImage_.GetWidth();
  height=LabelResultImage_.GetHeight();
  BIASCDOUT(D_LABEL, "image size: "<<width<<"x"<<height<<endl);
  //LabelResultImage_.FillImageWithConstValue(0);
  //coun++;
  int count = Label4Neighbour_();
  //ImageIO::Save("lri-"+os.str()+".mip", LabelResultImage_);
  
  BIASCDOUT(D_LABEL, "found "<<count<<" blobs\n");
  BlobCount_.clear();
  BlobMap_.clear();
  TopLefts_.clear(); 
  BotRights_.clear();
  Centers_.clear();
  MedianCenters_.clear();

  if (count<=0) { 
    BlobCount_.push_back(0);
  } else {
    //ImageIO::Save("sim.mip", LabelResultImage_);
    // initialize data structs
    BlobCount_.resize(count-1);
    TopLefts_.resize(count-1);
    BotRights_.resize(count-1);
    Centers_.resize(count-1);
    MedianCenters_.resize(count-1);
    SortVec_.resize(count);
    Num_.clear();
    for (int n=0; n<count; n++){
      SortVec_[n].bc=0;
      SortVec_[n].tl.Set(width, height);
      SortVec_[n].br.Set(0, 0);
      SortVec_[n].center.Set(0.0, 0.0);
      SortVec_[n].x.clear();
      SortVec_[n].y.clear();
    }
    // now get values
    LABEL_CALC_TYPE **ida=LabelResultImage_.GetImageDataArray();
    LABEL_CALC_TYPE c;
    unsigned int i;
    for (int y=0; y<height; y++){
      for (int x=0; x<width; x++){
        c=ida[y][x];
        // find the correct SortVec index
        i=0;
        while (i<Num_.size() && Num_[i]!=c) i++;
        // create SortVec it if not found
        if (i==Num_.size()) Num_.push_back(c);
        // fill SortVec_[i] 
        /*if (SortVec_[i].bc>10000){
          cout << SortVec_[i].bc<<" : c: "<<c<<"\ti="<<i<<"  "<<Num_.size()
          <<" "<<SortVec_.size()<<endl;
          }*/
        SortVec_[i].bc++; // BlobCount_[i]++;
        SortVec_[i].center += Vector2<double>((double)x, (double)y);
        SortVec_[i].x.push_back(x);
        SortVec_[i].y.push_back(y);
        if (SortVec_[i].tl[0]>x) SortVec_[i].tl[0]=x;
        if (SortVec_[i].tl[1]>y) SortVec_[i].tl[1]=y;
        if (SortVec_[i].br[0]<x) SortVec_[i].br[0]=x;
        if (SortVec_[i].br[1]<y) SortVec_[i].br[1]=y;
      }
    }
    
    sort(SortVec_.begin(), SortVec_.end());
    const Vector2<double> offset((double)minx, (double)miny);
    BIASCDOUT(D_LABEL, "neglecting bg label: "<<SortVec_[0].center <<"\t"
              << (double)SortVec_[0].bc << "\t"<<SortVec_[0].tl<<"\t"
              << SortVec_[0].br
              << (SortVec_[0].center / (double)SortVec_[0].bc) <<"\t"
              <<endl);
    static Median1D<int> med;
    double mx, my;
    for (int i=1, j=0; i<count; i++, j++){
      BlobCount_[j]=SortVec_[i].bc;
      BlobMap_[i]=SortVec_[i].bc;
      TopLefts_[j]=SortVec_[i].tl;
      BotRights_[j]=SortVec_[i].br;
      Centers_[j] = (SortVec_[i].center / (double)SortVec_[i].bc) + offset;
      med.Compute(SortVec_[i].x);
      mx = (double)med.GetMedian() + offset[0];
      med.Compute(SortVec_[i].y);
      my = (double)med.GetMedian() + offset[1];
      MedianCenters_[j] = Vector2<double>(mx, my);
      //cout << "c: "<<Centers_[j]<<"\tmc: "<<MedianCenters_[j]<<endl;
      BIASCDOUT(D_LABEL, SortVec_[i].center <<"\t"
                << (double)SortVec_[i].bc << "\t"<<Centers_[j] <<"\t"
                << SortVec_[i].tl <<"\t"<< SortVec_[i].br<< "\t"
                << (SortVec_[i].center / (double)SortVec_[i].bc) <<"\n");
    }
  }
  return 0;
}

int Label::Filter(const Image<unsigned char>& src, Image<unsigned char>& dst)
{
  BIASCDOUT(D_FILTERBASE_CALLSTACK, "Label::Filter\n");
 
  Process(src);
  if (BlobCount_.size()>255) {
    BIASERR("Label.Filter(): resulting image invalid, cause more than 255 "
            "regions");
  }
  bool UseLargeImage = false;
  Image<LABEL_CALC_TYPE> LargeImage;
  if (!LabelResultImage_.SamePixelAndChannelCount(src)) {
    LargeImage.Init(src.GetWidth(), src.GetHeight(), 1);
    unsigned int UpperLeftX, UpperLeftY, LowerRightX, LowerRightY;
    src.GetROI()->GetCorners(UpperLeftX, UpperLeftY, LowerRightX, LowerRightY);
    LargeImage.GetROI()->SetCorners(UpperLeftX, UpperLeftY, LowerRightX, LowerRightY);
    LargeImage.Paste2ROI(LabelResultImage_);
    UseLargeImage = true;
  }

  if (!dst.SamePixelAndChannelCount(src)) {
    if (!dst.IsEmpty()) dst.Release();
    dst.Init(src.GetWidth(), src.GetHeight(), 1);
  }
  
  if (UseLargeImage)  ImageConvert::ConvertST(LargeImage, dst,
                                              ImageBase::ST_unsignedchar);
  else ImageConvert::ConvertST(LabelResultImage_, dst,
                               ImageBase::ST_unsignedchar);
  return 0;
}

void Label::
GetBordersValid_(int& border_x, int& border_y) const
{
  border_x = border_y = 0;
}