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

/**
   @file
   @ingroup g_tools
   @brief load images given by the command line and compute mean value for each pixel, see biasaverage.cpp
   @author koeser, fkellner
*/

#include <Base/Image/ImageBase.hh>
#include <Base/Image/ImageIO.hh>
#include <Base/Image/ImageConvert.hh>
#include <Image/HomographyMapping.hh>
#include <Utils/Param.hh>
#include <Utils/IOUtils.hh>
#include <stdio.h>


using namespace BIAS;
using namespace std;


void usage()
{
  cout << "Usage: biasaverage [options] files"<<endl;
  cout <<" Options: "<< endl
       <<" -r         : rescale output by 2"<<endl
       << endl;
  cout <<"at each pixel position average values of all images"<<endl<<endl;
}

int DetermineFiletype(const string &fn, ImageIO::TFileFormat &thetype)
{
  string::size_type EndOfPrefix;
  string suffix;
  
  EndOfPrefix = fn.rfind('.');
  
  if (EndOfPrefix == string::npos) return -1;
  suffix = fn.substr(EndOfPrefix+1,fn.length());
  if (suffix == "pnm") thetype =  ImageIO::FF_ppm;
  else if (suffix == "ppm") thetype =  ImageIO::FF_ppm;
  else if (suffix == "pgm") thetype =  ImageIO::FF_pgm;
  else if (suffix == "mip") thetype =  ImageIO::FF_mip;
  else return -1;

  return 0;
}

int main (int argc, char *argv[])
{
  
  Param params(true);
  string* imageList = params.AddParamString("images", "list of images", "", 'i');
  string* output = params.AddParamString("output", "output filename (defaults to average.mip", "average.mip", 'o');
  string* varoutput = params.AddParamString("varoutput", "variance output filename (defaults to variance.mip", "variance.mip", 'O');
  string* outputUC = params.AddParamString("outputUC", "output filename of uchar image (defaults to average.pgm", "average.pgm", 'u');
  double* rescaleFactor = params.AddParamDouble("rescale", "rescale factor", 1.0, 0.001, 1000.0, 'r');
  bool* convertToRGB = params.AddParamBool("convertToRGB", "Convert to RGB before calculating", false,'c');
  bool* computeVariance = params.AddParamBool("computeVariance", "Also compute a variance image",false,'v');
  bool* analyzeSamples = params.AddParamBool("analyzeSamples", "Analyze 5 sample points",false,'a');
  
  if(!IOUtils::ParseCommandLineEvalHelp(params, argc, (char**)argv))
    return 0;

  if (imageList->length() == 0) {
    params.Usage();
    return 0; // for Dart testing 
  }
  
  vector<string> imageNames;
  if(Param::ParseListFile(*imageList, imageNames)!=0) {
    BIASERR("error parsing image names!" << *imageList);
    return -1;
  }
  
  bool rescale = false;
  if (*rescaleFactor != 1.0) {
    std::cout << "rescale"<<endl;
    rescale = true;
  }
  
  BIAS::ImageBase CurImage;
 
  BIAS::Image<float> averageImage;
  BIAS::Image<float> varianceImage;
  bool initialized = false;
  bool ucharout = true;

  int countImages = 0;
  for (unsigned int i=0;i<imageNames.size();i++) {
    if (ImageIO::Load(imageNames[i], CurImage) != 0) {
      BIASERR("Error reading "<< imageNames[i]);
      continue;
    }
    if(*convertToRGB){
      if(CurImage.GetColorModel() != ImageBase::CM_RGB){
        if (CurImage.GetStorageType() == ImageBase::ST_float) {
          BIAS::Image<float> in(CurImage),out;
          ImageConvert::ToRGB(in,out);    
          CurImage=out;
        }
        else if (CurImage.GetStorageType() == ImageBase::ST_unsignedchar) {
          BIAS::Image<unsigned char> in(CurImage),out;
          ImageConvert::ToRGB(in,out);    
          CurImage=out;
        }
      }
    }
    if (!initialized) {
      averageImage.Init(CurImage.GetWidth(), CurImage.GetHeight(), 
                        CurImage.GetChannelCount());
      averageImage.SetZero();
      varianceImage.Init(CurImage.GetWidth(), CurImage.GetHeight(),
                              CurImage.GetChannelCount());
      varianceImage.SetZero();

      initialized = true;
    }

    const unsigned int channels = CurImage.GetChannelCount();
    BIASASSERT(channels==1 || CurImage.IsInterleaved()); 
    BIASASSERT(averageImage.SamePixelAndChannelCount(CurImage));
    
    //cout << i << "/" << imageNames.size() << endl;
    
    void **ida = CurImage.GetImageDataArray();
    for (unsigned int y=0; y<CurImage.GetHeight(); y++) {

      for (unsigned int x=0; x<CurImage.GetWidth(); x++) {
        for (unsigned int c=0; c<channels; c++) {
                    //check storage type
          if (CurImage.GetStorageType() == ImageBase::ST_float) {
            averageImage.GetImageDataArray()[y][x*channels+c] +=
              ((float**)ida)[y][x*channels+c];
            if (ucharout) {
              ucharout = false;
              //cout << "-- Float input image found. disabling uchar image output" << endl;
            }
          } else if (CurImage.GetStorageType() == ImageBase::ST_unsignedchar) {
            averageImage.GetImageDataArray()[y][x*channels+c] +=
              float( ((unsigned char**)ida)[y][x*channels+c]) ;
          } else {
            BIASERR("Unknown storage type" << imageNames[i]);
          }
        }
      }
    }
    countImages++;
  }
  averageImage.ScaleShift(1.0/double(countImages),0);

  if(*computeVariance)
  {
    //now calculate variance image
    float **idaVar = varianceImage.GetImageDataArray();
    float **idaAv = averageImage.GetImageDataArray();
    for (unsigned int i=0;i<imageNames.size();i++) {

        if (ImageIO::Load(imageNames[i], CurImage) != 0) {
            BIASERR("Error reading "<< imageNames[i]);
            continue;
        }
        if(*convertToRGB){
            if(CurImage.GetColorModel() != ImageBase::CM_RGB){
                if (CurImage.GetStorageType() == ImageBase::ST_float) {
                    BIAS::Image<float> in(CurImage),out;
                    ImageConvert::ToRGB(in,out);
                    CurImage=out;
                }
                else if (CurImage.GetStorageType() == ImageBase::ST_unsignedchar) {
                    BIAS::Image<unsigned char> in(CurImage),out;
                    ImageConvert::ToRGB(in,out);
                    CurImage=out;
                }
            }
        }

        void **ida = CurImage.GetImageDataArray();
        const unsigned int channels = CurImage.GetChannelCount();
        for (unsigned int y=0; y<averageImage.GetHeight(); y++) {
            for (unsigned int x=0; x<averageImage.GetWidth(); x++) {
                for (unsigned int c=0; c<channels; c++) {
                    idaVar[y][x*channels+c] += pow( (idaAv[y][x*channels+c]-((float**)ida)[y][x*channels+c]),2);
                }
            }
        }
    }
    varianceImage.ScaleShift(1.0/double(countImages),0);
  }// end if(computeVariance)

  if(*analyzeSamples){
    unsigned width = averageImage.GetWidth();
    unsigned height = averageImage.GetHeight();
    float **idaVar = varianceImage.GetImageDataArray();
    float **idaAv = averageImage.GetImageDataArray();

    vector<Vector2<unsigned> > points;
    Vector2<unsigned> p(10,10);     points.push_back(p);
    p.Set(width-10,height-10);  points.push_back(p);
    p.Set(width/2,height/2);    points.push_back(p);
    p.Set(width-10,10);     points.push_back(p);
    p.Set(10,height-10);    points.push_back(p);
    for(unsigned i=0;i<points.size();i++){
        cout<<"Mean at "<<points[i]<<": "<<idaAv[points[i][1]][points[i][0]]<<endl;
        cout<<"Variance at "<<points[i]<<": "<<idaVar[points[i][1]][points[i][0]]<<endl;
    }

    unsigned channels = averageImage.GetChannelCount();
    float *mean = new float[channels];
    float *var = new float[channels];
    averageImage.GetMeanPixelValue(mean);

    varianceImage.GetMeanPixelValue(var);
    cout<<"Image mean and variance variance : ";
    for(unsigned i=0;i<channels;i++) cout<<mean[i]<<" "<<var[i]<<endl;
    delete[] mean;
    delete[] var;
  }

  Image<float> rescaled,varRescaled;
  HomographyMapping<float, float> converter;
  HMatrix H(MatrixIdentity);
  if (rescale) {
    H[0][0] = H[1][1] = 1.0 / (*rescaleFactor);
    rescaled.Init((unsigned int)ceil( (*rescaleFactor)*double(averageImage.GetWidth()) ),
            (unsigned int)ceil( (*rescaleFactor)*double(averageImage.GetHeight()) ), averageImage.GetChannelCount());
    if(*computeVariance){
        varRescaled.Init((unsigned int)ceil( (*rescaleFactor)*double(averageImage.GetWidth()) ),
                (unsigned int)ceil( (*rescaleFactor)*double(averageImage.GetHeight()) ), averageImage.GetChannelCount());
    }
    converter.SetHomography(H);
    converter.Map(averageImage, rescaled, MapBicubic);
    if(*computeVariance){
        converter.Map(varianceImage, varRescaled, MapBicubic);
    }
  } else {
    rescaled = averageImage;
    if(*computeVariance)
        varRescaled = varianceImage;
  }
  if(*convertToRGB == false) {
    rescaled.SetColorModel(CurImage.GetColorModel());
  }
  ImageIO::Save(*output, rescaled);
  if(*computeVariance)
    ImageIO::Save(*varoutput, varRescaled);

  if (ucharout) {
    Image<unsigned char> ucImg;
    ImageConvert::ConvertST(rescaled, ucImg, ImageBase::ST_unsignedchar);
    if(*convertToRGB == false) {
      ucImg.SetColorModel(CurImage.GetColorModel());
    }
    ImageIO::Save(*outputUC, ucImg);
  }
  return 0;
    
}