bindata2dctbindata.cc

 const char *help = "\
 progname: bindata2dctbindata.cc\n\
 code2html: This program reads a bindata and decomposes it in DCT/DCTmod2 blocks.\n\
 version: Torch3 vision2.0, 2004-2005\n\
 (c) Sebastien Marcel (marcel@idiap.ch)\n";
 
 #include "ImageGray.h"
 #include "ipDCT2D.h"
 #include "ipBlock.h"
 #include "ipDCTmod2.h"
 #include "DiskXFile.h"
 #include "CmdLine.h"
 
 using namespace Torch;
 
 // Allowed dimension for DCT mod2
 const int allowed_dim[] = {6, 10, 15, 21, 28, 36, 43, 49, 54, 58, 61, 63, 64};
 
 int main(int argc, char *argv[])
 {
 	char *filename_in;
 	int width;
 	int height;
 	char *filename_out;
 	int dct_dim;
 	int delta_neighbour;
 	int delta_keep;
 	int block_size;
 	int block_overlap;
 	bool block_to_frame;
 	bool dctmod2;
 	bool force;
 	bool verbose;
 	bool display;
 	bool xy;
 	bool xyIsEmbedded;
 	bool blocknorm;
 	bool unnorm;
 	
 	CmdLine cmd;
 	cmd.setBOption("write log", false);
 	cmd.info(help);
 	cmd.addText("\nArguments:");
 	cmd.addSCmdArg("imagefile in", &filename_in, "image file in");
 	cmd.addICmdArg("width", &width, "width");
 	cmd.addICmdArg("height", &height, "height");
 	cmd.addSCmdArg("output filename", &filename_out, "output bindata filename");
 	cmd.addText("\nBlock decomposition Options:");
 	cmd.addBCmdOption("-blocktoframe", &block_to_frame, false, "a block becomes a frame");
 	cmd.addICmdOption("-block", &block_size, 8, "block size");
 	cmd.addICmdOption("-overlap", &block_overlap, 4, "block overlap");
 	cmd.addText("\nDCT Options:");
 	cmd.addICmdOption("-dctdim", &dct_dim, 15, "DCT 2D dimension");
 	cmd.addICmdOption("-delta", &delta_neighbour, 1, "delta neighbour");
 	cmd.addICmdOption("-deltakeep", &delta_keep, 3, "delta keep");
 	cmd.addBCmdOption("-dctmod2", &dctmod2, false, "apply dctmod2");
 	cmd.addBCmdOption("-force", &force, false, "force dctmod2 for blocks larger that 8x8");
 	cmd.addText("\nExtra Options:");
 	cmd.addBCmdOption("-xy", &xy, false, "add xy coordinates of blocks");
 	cmd.addBCmdOption("-xyprovided", &xyIsEmbedded, false, "xy coordinates are provided");
 	cmd.addBCmdOption("-blocknorm", &blocknorm, false, "mean/stdv normalize blocks");
 	cmd.addBCmdOption("-unnorm", &unnorm, false, "unnormalize input image");
 	cmd.addText("\nMisc Options:");
 	cmd.addBCmdOption("-display", &display, false, "display dct features");
 	cmd.addBCmdOption("-verbose", &verbose, false, "verbose");
 	cmd.read(argc, argv);
 	
 	//
 	if(verbose) print("DCT 2D feature extraction.\n");
 
 	if((dctmod2 == true) && (block_size != 8))
 	{
 		warning("DCTmod2 is usually designed only for 8x8 blocks.");
 		
 		if(force == false) 
 		{
 			print("Disabling dctmod2.\n");
 
 		   	dctmod2 = false;
 		}
 	}
 
 	if(dctmod2 == true)
 	{
 		int dim_ok = 0;
 		int max_ = sizeof(allowed_dim)/sizeof(int);
 		
 		for(int i = 0 ; i < max_ ; i++)
 			if(dct_dim == allowed_dim[i]) dim_ok = 1;
 	
 		if(!dim_ok)
 		{
 		   	error("DCT dim (%d) is incorrect.", dct_dim);
 		   	print("   Correct values are: ");
 			for(int i = 0 ; i < max_ ; i++)
 			   	print("%d ", allowed_dim[i]);
 	  		print("\n");
 	
 			exit(0);
 		}
 	}
   	
 	//
 	int n_images;
 	int image_size;
 	DiskXFile *pf_in;
 	int n_inputs;
 	int n_patterns;
 	DiskXFile *pf_out;
 
 	// Load the image
 	pf_in = new DiskXFile(filename_in, "r");
 	pf_in->read(&n_images, sizeof(int), 1);
 	pf_in->read(&image_size, sizeof(int), 1);
 
 	if(verbose)
 	{
 		print("Input file ...\n");
 		print("   n_inputs = %d\n", image_size);
 		print("   n_patterns = %d\n", n_images);  
 	}
 
 	int offset = 0;
 	
 	if(xyIsEmbedded) offset = 2;
 
 	if(image_size-offset != width * height)
 		error("Incorrect image size %d != %dx%d.", image_size, width, height);
 	
 	Sequence *seqimage = new Sequence(1, image_size-offset);
 	
 	// feature image machine to apply for each block
 	ipCore *ip = NULL;
 		
 	int n_block, n_output;
 	int rowblocks;
 	int colblocks;
 	
 	// Choose the block by block convolution
 	if(dctmod2 == true)
 	{
 		ip = new ipDCTmod2(width, height, "gray", dct_dim, block_size, block_overlap, delta_neighbour, delta_keep);
 
 		ipDCTmod2 *ip_ = (ipDCTmod2 *)ip;
 		n_block = ip_->n_block;
 		n_output = ip_->output_size;
 		rowblocks = ip_->getrowblocks();
 		colblocks = ip_->getcolblocks();
 	}
 	else
 	{
 		ipCore *ip_dct = NULL;
 		ip_dct = new ipDCT2D(block_size, "gray", dct_dim);
 		ip_dct->setBOption("verbose", verbose);
 
 		ip = new ipBlock(width, height, "gray", ip_dct, dct_dim, block_size, block_overlap);
 
 		ipBlock *ip_ = (ipBlock *)ip;
 		n_block = ip_->n_block;
                 n_output = dct_dim;
 		rowblocks = colblocks = block_size;
 	}
 	ip->setBOption("verbose", verbose);
 
 	if(verbose)
 	{
 		print("-> n_block = %d\n", n_block);
 		print("-> n_output = %d\n", n_output);
 	}
 
 	//
 	pf_out = new DiskXFile(filename_out, "w");
 	if(pf_out == NULL)
 	{
   		error("Opening bindata file %s", filename_out);      
     		return 0;
 	}
 	
 	//
 	int frame_size;
 	
 	if(block_to_frame)
 	{
 	  	if(verbose) print("Block to frame: each block becomes a frame.\n");
 		
 		n_inputs = n_output + offset;
 		if(xy) frame_size = n_inputs + 2;
 		else frame_size = n_inputs;
 
 		n_patterns = n_images * n_block;
 	}
 	else
 	{
 		if(xy) warning("Please use blocktoframe to add xy coordinates.");
 		
 		n_inputs = n_block * (n_output + offset);
 		frame_size = n_inputs;
 		n_patterns = n_images;
 	}
 
 	if(verbose)
 	{
 		print("Output file ...\n");
 		print("   n_inputs = %d\n", frame_size);
 		print("   n_patterns = %d\n", n_patterns);  
 	}
 
 	//
 	pf_out->write(&n_patterns, sizeof(int), 1);
 	pf_out->write(&frame_size, sizeof(int), 1);
 	
 	//
 	real *blocks_mean = NULL;
 	real *blocks_stdv = NULL;
 
 	real *data_in = new real [image_size];
 	real *data_out = new real [n_inputs];
 
 	for(int p = 0 ; p < n_images ; p++)
 	{
 	   	pf_in->read(data_in, sizeof(real), image_size);
 		
 		if(unnorm)
 		{
 			for(int i = offset ; i < image_size ; i++)
 			{
 				int z_ = (int) (data_in[i] * 255.0);
 
 				seqimage->frames[0][i-offset] = z_;
 			}
 		}
 		else
 		{
 			for(int i = offset ; i < image_size ; i++)
 				seqimage->frames[0][i-offset] = data_in[i];
 		}
 
 		// Apply the pattern machine (DCT) to each block of the image
 		ip->process(seqimage);
 
 		if(blocknorm)
 		{
 			blocks_mean = new real [n_output];
 			blocks_stdv = new real [n_output];
 			
 			for(int i = 0 ; i < n_output ; i++)
 			{
 				blocks_mean[i] = 0;
 				blocks_stdv[i] = 0;
 			}
 			
 			for(int b = 0 ; b < n_block ; b++)
 			{
 				real *src_ = ip->seq_out->frames[b];
 				
 				for(int i = 0 ; i < n_output ; i++)
 				{
 					real z = src_[i];
 					blocks_mean[i] += z;
 					blocks_stdv[i] += z*z;
 				}
 			}
 
 			for(int i = 0 ; i < n_output ; i++)
 			{
 				blocks_mean[i] /= (real) n_block;
 				blocks_stdv[i] /= (real) n_block;
 				blocks_stdv[i] -= blocks_mean[i]*blocks_mean[i];
 				if(blocks_stdv[i] <= 0)
 				{
 					warning("MeanVarNorm: input column %d has a null stdv. Replaced by 1.", i);
 					blocks_stdv[i] = 1.;
 				}
 				else blocks_stdv[i] = sqrt(blocks_stdv[i]);
 			}
 			
 			for(int b = 0 ; b < n_block ; b++)
 			{
 				real *src_ = ip->seq_out->frames[b];
 				
 				for(int i = 0 ; i < n_output ; i++)
 					src_[i] = (src_[i] - blocks_mean[i]) / blocks_stdv[i];
 			}
 		}
 
 		if(block_to_frame)
 		{
 			for(int b = 0 ; b < n_block ; b++)
 			{
 			   	// (July, 22 2005) There might be a BUG there: TO CHECK
 			   	int x = (b % rowblocks);
 				int y = (b / rowblocks);
 				// the correct thing could be
 			   	//int x = (b % colblocks);
 				//int y = (b / colblocks);
 				
 		   		if(display == true) 
 				{
 					if(xy) print("Block %d (%dx%d): ", b, x, y);
 					else print("Block %d: ", b);
 				}
 			
 				real *output_ = ip->seq_out->frames[b];
 	
 				// if xy is embedded
 				for(int i = 0 ; i < offset ; i++)
 				{
 		   			if(display == true) print("%g ", data_in[i]);
 					
 					data_out[i] = data_in[i];
 				}
 				
 				for(int i = 0 ; i < n_output ; i++)
 				{
 		   			if(display == true) print("%g ", output_[i]);
 
 		   			data_out[i+offset] = output_[i];
 				}
 
 		   		if(display == true) print("\n");
 	
 				if(xy)
 				{
 					real x_ = (real) x;
 					real y_ = (real) y;
 
 					pf_out->write(&x_, sizeof(real), 1);
 					pf_out->write(&y_, sizeof(real), 1);
 				}
 				
 				for(int i = 0 ; i <  n_inputs ; i++)
 					pf_out->write(&data_out[i], sizeof(real), 1);   
 			}
 		}
 		else
 		{
 			int j = 0;
 		
 			for(int b = 0 ; b < n_block ; b++)
 			{
 		   		if(display == true) print("Block %d: ", b);
 			
 				real *output_ = ip->seq_out->frames[b];
 			
 				// if xy is embedded
 				for(int i = 0 ; i < offset ; i++)
 				{
 		   			if(display == true) print("%g ", data_in[i]);
 					
 					data_out[i] = data_in[i];
 				}
 				
 				for(int i = offset ; i < n_output ; i++)
 				{
 		   			if(display == true) print("%g ", output_[i]);
 
 		   			data_out[j] = output_[i];
 				
 					j++;
 				}
 		   		if(display == true) print("\n");
 			}
 		
 			for(int i = 0 ; i <  n_inputs ; i++)
 				pf_out->write(&data_out[i], sizeof(real), 1);   
 		}
 	}
 		
 	if(blocknorm)
 	{
 		delete [] blocks_mean;
 		delete [] blocks_stdv;
 	}
 	
 	delete [] data_out;
 	delete [] data_in;
 	delete pf_out;
 	
 
 	// Free all
 	delete ip;
 	delete seqimage;
 
 	return 0;
 }