multiscaleRetinex.cc

  1 const char *help = "\
  2 progname: multiscaleretinex.cc\n\
  3 code2html: This program reads a pgm image (grayscale) and normalizes the lighting conditions according to the Multiscale Retinex approach.\n\
  4 \n\
  5 The input image is filtered (gaussian) several times, and the enhanced version\n\
  6 is obtained by dividing the input image by the filtered ones, as described\n\
  7 in the article by Rahmann, Woodel and Jobson:\n\
  8 'A comparison of the MSR with other image enhancement techniques'.\n\
  9 \n\
 10 version: Torch3 vision2.0, 2004-2005\n\
 11 (c) Guillaume Heusch (heusch@idiap.ch)\n";
 12 
 13 // image
 14 #include "ImageGray.h"
 15 
 16 // ipCore
 17 #include "ipGaussian.h"
 18 #include "ipHistoEqual.h"
 19 
 20 // misc
 21 #include "Timer.h"
 22 #include "DiskXFile.h"
 23 #include "CmdLine.h"
 24 
 25 using namespace Torch;
 26 
 27 void rescale(real *data, int width, int height);
 28 
 29 int main(int argc, char **argv)
 30 {
 31   
 32   Timer timer;
 33 
 34 	char *image_filename;
 35 	char *result_filename;
 36 	int nb_filters = 3;
 37 	int min_size;
 38 	int middle_size;
 39 	int max_size;
 40 	bool verbose;
 41 	bool histo;
 42 
 43 	
 44 	// ------------------------ COMMAND LINE ---------------------------------------------------------------------
 45   	CmdLine cmd;
 46 	cmd.setBOption("write log", false);
 47   	cmd.info(help);
 48 
 49   	cmd.addText("\nArguments:");
 50   	cmd.addSCmdArg("image_filename", &image_filename, "input image filename");
 51 	cmd.addSCmdArg("result_filename", &result_filename, "result filename");
 52 	
 53   	cmd.addText("\nOptions:");
 54   	cmd.addBCmdOption("-verbose", &verbose, false, "verbose");
 55 	cmd.addICmdOption("-min_size", &min_size, 3, "size of the smallest convolution kernel (must be odd)");
 56 	cmd.addICmdOption("-middle_size", &middle_size, 9, "size of the 'middle' convolution kernel (must be odd)");
 57 	cmd.addICmdOption("-max_size", &max_size, 15, "size of the biggest convolution kernel (must be odd)");
 58 	cmd.addBCmdOption("-histo", &histo, false, "apply histogram equalization on the result");
 59 	
 60 	cmd.read(argc, argv);
 61 
 62 
 63 	Allocator *allocator = new Allocator;
 64 	// ----------------------- LOAD IMAGE --------------------------------------------------------------------------
 65   	DiskXFile *image_file = NULL;
 66 	Image *image_in = NULL;
 67 
 68 	image_in = new(allocator) ImageGray();
 69 	image_in->setBOption("verbose", verbose);
 70 	
 71 	image_file = new(allocator) DiskXFile(image_filename, "r");
 72 	image_in->loadXFile(image_file);
 73 
 74 	if(verbose)
 75 	{
 76 		print("Image info:\n");
 77 		print("   width = %d\n", image_in->width);
 78 		print("   height = %d\n", image_in->height);
 79 		print("   format = %s (%d)\n", image_in->coding, image_in->n_planes);
 80 	}
 81 
 82 	
 83 	// ---------------- MULTI-SCALE GAUSSIAN FILTERING -------------------------------------------------------------
 84 	Image *filtered_array [nb_filters];
 85 	ipCore *gaussianMachine = NULL;
 86 
 87 	for (int filter_index = 0; filter_index < nb_filters; filter_index++)
 88 	  filtered_array[filter_index] = NULL;
 89 	
 90 
 91 	for (int filter_index = 0; filter_index < nb_filters; filter_index++) {
 92 	  
 93 	  int size = 0;
 94 	  
 95 	  if (filter_index == 0)
 96 	    size = min_size;
 97 	  if (filter_index == 1)
 98 	    size = middle_size;
 99 	  if (filter_index == 2)
100 	    size = max_size; 
101 
102 	  gaussianMachine = new(allocator) ipGaussian(size, image_in->width, image_in->height, "gray");
103 	  gaussianMachine->process(image_in);
104 	  
105 	  filtered_array[filter_index] = new(allocator) ImageGray();
106 	  filtered_array[filter_index]->copyFrom(image_in->width,image_in->height, gaussianMachine->seq_out->frames[0], "gray", 255);
107 
108 	  // ************************************************
109 	  // IF YOU WANT TO SAVE THE FILTERED IMAGES 
110 	  // -----------------------------------------------
111 	  //char filename[50];
112 	  //sprintf(filename, "filtered_%ix%i", size, size);
113 	  //strcat(filename, ".pgm");
114 	  //rescale(filtered_array[filter_index]->data, filtered_array[filter_index]->width, filtered_array[filter_index]->height);
115 	  //filtered_array[filter_index]->updatePixmapFromData();
116 	  //filtered_array[filter_index]->save(filename); 
117 	  // ***********************************************
118 	}
119    
120 
121 
122 	// ---------------------- BUILD AND SAVE THE RESULT --------------------------------------------
123 	Image *result = NULL;
124 	result = new(allocator) ImageGray(image_in->width, image_in->height);
125 	result->setBOption("verbose", verbose);
126 
127 	for (int n_image = 0; n_image < nb_filters; n_image++) { 
128 	  for (int i=0; i<(result->width*result->height); i++)
129 	    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)... 
130 	}
131 
132 	//rescale
133 	rescale(result->data, result->width, result->height);
134 
135 	// histogram equalization if specified
136 	ipCore *histoEq = NULL;
137    
138 	if (histo) {   
139 	  histoEq = new(allocator) ipHistoEqual(image_in->width, image_in->height, "gray");
140 	  histoEq->process(result);
141 	}
142 
143 	Image *image_out = NULL;
144 	image_out = new(allocator) ImageGray();
145 	image_out->setBOption("verbose", verbose);
146  
147 	if (histo)
148 	  image_out->copyFrom(image_in->width, image_in->height, histoEq->seq_out->frames[0], "gray", 255);
149 	else
150 	  image_out->copyFrom(image_in->width, image_in->height, result->data, "gray", 255);
151 	
152 	image_out->updatePixmapFromData();
153 	image_out->save(result_filename);
154 
155 	
156 	// ----------------------- CLEAN UP -----------------------------------------------------------------------
157 	delete allocator;
158 	
159 	print("time elapsed: %g\n", timer.getTime());
160 	timer.stop();
161 	
162 	return(0);
163 }
164 
165 
166 void rescale(real *data, int width, int height) {
167 
168   real max_val = -100000;
169   real min_val = 100000;
170   real range;
171   int max_idx;
172 
173   // find min and max values in the image
174   for (int i=0; i<(width*height); i++) {
175 
176     if (data[i] > max_val) {
177       max_val = data[i];
178       max_idx = i;
179     }
180     if (data[i] < min_val)
181       min_val = data[i];  
182   }
183   
184   range = max_val - min_val;
185 
186   // change the scale
187   for (int i=0; i<(width*height); i++) {
188     data[i] = (data[i] - min_val)*(1.0/range);
189     data[i] = data[i]*255.0;
190   } 
191 }