multiscaleretinex.cc

torch tracking code, it is a good code

const char *help = "\
progname: multiscaleretinex.cc\n\
code2html: This program reads a pgm image (grayscale) and normalizes the lighting conditions according to the Multiscale Retinex approach.\n\
\n\
The input image is filtered (gaussian) several times, and the enhanced version\n\
is obtained by dividing the input image by the filtered ones, as described\n\
in the article by Rahmann, Woodel and Jobson:\n\
'A comparison of the MSR with other image enhancement techniques'.\n\
\n\
version: Torch3 vision2.0, 2004-2005\n\
(c) Guillaume Heusch ([e-mail])\n";

// image
#include "ImageGray.h"

// ipCore
#include "ipGaussian.h"
#include "ipHistoEqual.h"

// misc
#include "Timer.h"
#include "DiskXFile.h"
#include "CmdLine.h"

using namespace Torch;

void rescale(real *data, int width, int height);

int main(int argc, char **argv)
{
  
  Timer timer;

	char *image_filename;
	char *result_filename;
	int nb_filters = 3;
	int min_size;
	int middle_size;
	int max_size;
	bool verbose;
	bool histo;

	
	// ------------------------ COMMAND LINE ---------------------------------------------------------------------
  	CmdLine cmd;
	cmd.setBOption("write log", false);
  	cmd.info(help);

  	cmd.addText("\nArguments:");
  	cmd.addSCmdArg("image_filename", &image_filename, "input image filename");
	cmd.addSCmdArg("result_filename", &result_filename, "result filename");
	
  	cmd.addText("\nOptions:");
  	cmd.addBCmdOption("-verbose", &verbose, false, "verbose");
	cmd.addICmdOption("-min_size", &min_size, 3, "size of the smallest convolution kernel (must be odd)");
	cmd.addICmdOption("-middle_size", &middle_size, 9, "size of the 'middle' convolution kernel (must be odd)");
	cmd.addICmdOption("-max_size", &max_size, 15, "size of the biggest convolution kernel (must be odd)");
	cmd.addBCmdOption("-histo", &histo, false, "apply histogram equalization on the result");
	
	cmd.read(argc, argv);


	Allocator *allocator = new Allocator;
	// ----------------------- LOAD IMAGE --------------------------------------------------------------------------
  	DiskXFile *image_file = NULL;
	Image *image_in = NULL;

	image_in = new(allocator) ImageGray();
	image_in->setBOption("verbose", verbose);
	
	image_file = new(allocator) DiskXFile(image_filename, "r");
	image_in->loadXFile(image_file);

	if(verbose)
	{
		print("Image info:\n");
		print("   width = %d\n", image_in->width);
		print("   height = %d\n", image_in->height);
		print("   format = %s (%d)\n", image_in->coding, image_in->n_planes);
	}

	
	// ---------------- MULTI-SCALE GAUSSIAN FILTERING -------------------------------------------------------------
	Image** filtered_array = new Image*[nb_filters];
	ipCore *gaussianMachine = NULL;

	for (int filter_index = 0; filter_index < nb_filters; filter_index++)
	  filtered_array[filter_index] = NULL;
	

	for (int filter_index = 0; filter_index < nb_filters; filter_index++) {
	  
	  int size = 0;
	  
	  if (filter_index == 0)
	    size = min_size;
	  if (filter_index == 1)
	    size = middle_size;
	  if (filter_index == 2)
	    size = max_size; 

	  gaussianMachine = new(allocator) ipGaussian(size, image_in->width, image_in->height, "gray");
	  gaussianMachine->process(image_in);
	  
	  filtered_array[filter_index] = new(allocator) ImageGray();
	  filtered_array[filter_index]->copyFrom(image_in->width,image_in->height, gaussianMachine->seq_out->frames[0], "gray", 255);

	  // ************************************************
	  // IF YOU WANT TO SAVE THE FILTERED IMAGES 
	  // -----------------------------------------------
	  //char filename[50];
	  //sprintf(filename, "filtered_%ix%i", size, size);
	  //strcat(filename, ".pgm");
	  //rescale(filtered_array[filter_index]->data, filtered_array[filter_index]->width, filtered_array[filter_index]->height);
	  //filtered_array[filter_index]->updatePixmapFromData();
	  //filtered_array[filter_index]->save(filename); 
	  // ***********************************************
	}
   


	// ---------------------- BUILD AND SAVE THE RESULT --------------------------------------------
	Image *result = NULL;
	result = new(allocator) ImageGray(image_in->width, image_in->height);
	result->setBOption("verbose", verbose);

	for (int n_image = 0; n_image < nb_filters; n_image++) { 
	  for (int i=0; i<(result->width*result->height); i++)
	    result->data[i] += (log((image_in->data[i]+1)) - log((filtered_array[n_image]->data[i]+1))); // add 1 in order to avoid log(0)... 
	}

	//rescale
	rescale(result->data, result->width, result->height);

	// histogram equalization if specified
	ipCore *histoEq = NULL;
   
	if (histo) {   
	  histoEq = new(allocator) ipHistoEqual(image_in->width, image_in->height, "gray");
	  histoEq->process(result);
	}

	Image *image_out = NULL;
	image_out = new(allocator) ImageGray();
	image_out->setBOption("verbose", verbose);
 
	if (histo)
	  image_out->copyFrom(image_in->width, image_in->height, histoEq->seq_out->frames[0], "gray", 255);
	else
	  image_out->copyFrom(image_in->width, image_in->height, result->data, "gray", 255);
	
	image_out->updatePixmapFromData();
	image_out->save(result_filename);

	
	// ----------------------- CLEAN UP -----------------------------------------------------------------------
	delete[] filtered_array;
	delete allocator;
	
	print("time elapsed: %g\n", timer.getTime());
	timer.stop();
	
	return(0);
}


void rescale(real *data, int width, int height) {

  real max_val = -100000;
  real min_val = 100000;
  real range;
  int max_idx;

  // find min and max values in the image
  for (int i=0; i<(width*height); i++) {

    if (data[i] > max_val) {
      max_val = data[i];
      max_idx = i;
    }
    if (data[i] < min_val)
      min_val = data[i];  
  }
  
  range = max_val - min_val;

  // change the scale
  for (int i=0; i<(width*height); i++) {
    data[i] = (data[i] - min_val)*(1.0/range);
    data[i] = data[i]*255.0;
  } 
}