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

/**
 * @example ExampleJointHistrogram.cpp
   @relates Image, HistogramImage, JointHistogram
   @brief Usage for the JointHistogram class and creation of
          Normalized Mutual Information
   @ingroup g_examples
   @author  ischiller
*/


#include <Base/Image/Image.hh>
#include <Base/Image/ImageIO.hh>
#include <Base/Image/ImageConvert.hh>
#include <Base/ImageUtils/ImageDraw.hh>
#include <Base/Math/Random.hh>
#include <Image/JointHistogram.hh>
#include <Image/HistogramImage.hh>
#include <iostream>
#include <Filter/GradientGauss.hh>
#include <Filter/GradientSimple.hh>

using namespace BIAS;
using namespace std;

int main(int argc, char *argv[])
{
  int res = 0;
  Image<unsigned char> image1, image2;
  Image<float> gauss1, gauss2, gauss1X, gauss1Y,gauss2X,gauss2Y;

  JointHistogram<float> jhist;
  BIAS::GradientSimple<unsigned char, float> gradientGenerator;

  if (argc <5) {
    cerr<< argv[0] <<" <reference image> <output> <usegradient> filenames"<<endl;
    exit(1);
  }

  //if ((res=ImageIO::Load(argv[1], image1)) != 0){
  if ((res=ImageIO::Load(argv[1], image1)) != 0){
    cerr << "error loading file " << argv[1] << "   " << res << endl;
    exit (-1);
  }

  if (image1.GetChannelCount() != 1){
    cerr << "wrong channel count in image"<< endl;
    BIAS::ImageConvert::IP_ToGrey(image1);

  }
  int numberOfImages = argc-1;
  std::vector<string> images;
  float nmi=0.0;
  float maxNMI=-100.0;
  float minJointEntropy=100000.0;
  string minJointEntropyImageName;
  string maxNMIImageName;
  //go over all images and calculate nmi
  for(int i=4;i<=numberOfImages;i++){

    //if ((res=ImageIO::Load(argv[i], image2)) != 0){
    if ((res=ImageIO::Load(argv[i], image2)) != 0){
      cerr << "error loading file " << argv[1] << "   " << res << endl;
      exit (-1);
    }
    else
      cout<<"\nUsing image: "<<argv[i]<<" for comparison"<<endl;

    if (image2.GetChannelCount() != 1){
      cerr << "wrong channel count in image"<< endl;
      BIAS::ImageConvert::IP_ToGrey(image2);
    }

    if(atoi(argv[3]) > 0){
      gradientGenerator.Filter(image1, gauss1X, gauss1Y,  gauss1);
      gradientGenerator.Filter(image2, gauss2X, gauss2Y, gauss2);
      if(atoi(argv[3]) == 2){
        cout<<"Using angle directions"<<endl;
        float** idaX1 = gauss1X.GetImageDataArray();
        float** idaY1 = gauss1Y.GetImageDataArray();
        float** dst1 = gauss1.GetImageDataArray();
        float** idaX2 = gauss2X.GetImageDataArray();
        float** idaY2 = gauss2Y.GetImageDataArray();
        float** dst2 = gauss2.GetImageDataArray();
        double angle =0.0;
        double gradientThreshold = 30;
        for(unsigned int x=0;x<image1.GetWidth();x++){
          for(unsigned int y=0;y<image1.GetHeight();y++){
            double YGradValue = idaY1[y][x];
            double XGradValue = idaX1[y][x];
            //calulate the angle of the gradient
            angle = atan2(YGradValue,XGradValue);

            // // bring angle to [0..2pi]
            while (angle<0.0f) angle += 2.0f * (float)M_PI;
            if (angle>float(2.0f*M_PI)) angle = 2.0f * (float)M_PI;
            angle = (angle*180)/M_PI;

            if(dst1[y][x] >= gradientThreshold)
              dst1[y][x] = (float)angle;
            else
              dst1[y][x] = (float)0;

            //second image
            YGradValue = idaY2[y][x];
            XGradValue = idaX2[y][x];
            angle = atan2(YGradValue,XGradValue);

            // bring angle to [0..2pi]
            while (angle<0.0f) angle += 2.0f * (float)M_PI;
            if (angle>float(2.0f*M_PI)) angle = 2.0f * (float)M_PI;

            angle = (angle*180)/M_PI;
            if(dst2[y][x] >= gradientThreshold)
              dst2[y][x] = (float)angle;
            else
              dst2[y][x] = (float)0;
          }//end for
        }// end for
      }
    }
    else
      {
        BIAS::ImageConvert::ConvertST(image1,gauss1, gauss1.GetStorageType());
        BIAS::ImageConvert::ConvertST(image2,gauss2, gauss2.GetStorageType());
      }

    jhist.Compute(gauss1,gauss2);
    std::string filename;
    filename = argv[2];
    filename +="-";
    filename +=argv[i];

    BIAS::Image<float> tmpImage;
    jhist.Draw(tmpImage);

    //BIAS::ImageIO::Save(filename.c_str(),tmpImage);
    BIAS::ImageIO::Save(filename.c_str(),tmpImage);

    float entropy = jhist.CalcShannonEntropy();
    cout<<"Shannon Joint Entropy :"<<entropy<<endl;
    if(entropy < minJointEntropy)
      {
        minJointEntropy = entropy;
        minJointEntropyImageName = argv[i];
        BIAS::Image<float> tmpImage;
        jhist.Draw(tmpImage);
        //BIAS::ImageIO::Save(argv[2],tmpImage);
        BIAS::ImageIO::Save(argv[2],tmpImage);
      }

    HistogramImage hist;
    hist.InitHist(256,2);

    hist.AddHist(gauss1,0);
    hist.AddHist(gauss2,1);

    float entropy0 = hist.CalcShannonEntropy(0);
    float entropy1 = hist.CalcShannonEntropy(1);
    cout<<"Shannon Entropies  :"<<entropy0<<","<<entropy1<<endl;

    hist.Draw();
    //  BIAS::ImageIO::Save("hist.mip",hist);

    nmi = (entropy0 + entropy1)/entropy;
    if(nmi > maxNMI)
      {
        maxNMI = nmi;
        maxNMIImageName = argv[i];
      }

    cout<<"Normalized Mutual Information: "<<nmi<<endl;

  }// end for

  cout<<"\n\n***************************************************************\n";
  cout<<"*******************  Final results: *************************** \n";
  cout<<"ImageName with maximum NMI: "<<maxNMIImageName<<endl;
  cout<<"MaxNMI :                    "<<maxNMI<<endl;
  cout<<"ImageName with min JointEntropy: "<<minJointEntropyImageName<<endl;
  cout<<"Min Joint Entropy :              "<<minJointEntropy<<endl;
  cout<<"***************************************************************\n";
  return 0;
}