bindata2dctbindata.cc

torch tracking code, it is a good code

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 ([e-mail])\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;
	delete pf_in;

	// Free all
	delete ip;
	delete seqimage;

	return 0;
}