blockdct2dpgm.cc

torch tracking code, it is a good code

const char *help = "\
progname: blockDCT2Dpgm.cc\n\
code2html: This program reads a pgm image, decomposes it in blocks and applies DCT to each block.\n\
version: Torch3 vision2.0, 2004-2005\n\
(c) Sebastien Marcel ([e-mail])\n";

#include "ImageGray.h"
#include "ipDCT2D.h"
#include "ipDCT2Dinverse.h"
#include "ipBlock.h"
#include "DiskXFile.h"
#include "CmdLine.h"

using namespace Torch;

int main(int argc, char **argv)
{
	char *image_filename;
	int block_size;
	int block_overlap;
	int dimDCT;
	bool verbose;
	bool inverse;
	bool printblocks;
  
	
  	CmdLine cmd;
	cmd.setBOption("write log", false);
  	cmd.info(help);
  	cmd.addText("\nArguments:");
  	cmd.addSCmdArg("image filename", &image_filename, "image filename");
  	cmd.addText("\nOptions:");
  	cmd.addICmdOption("-blocksize", &block_size, 8, "size of blocks");
  	cmd.addICmdOption("-dim", &dimDCT, 15, "DCT dimensionality");
  	cmd.addICmdOption("-blockoverlap", &block_overlap, 4, "overlap between blocks");
  	cmd.addBCmdOption("-verbose", &verbose, false, "verbose");
  	cmd.addBCmdOption("-inverse", &inverse, false, "inverse");
  	cmd.addBCmdOption("-pb", &printblocks, false, "print blocks");
	cmd.read(argc, argv);


	Image *image_in = NULL;

	image_in = new ImageGray();
	image_in->setBOption("verbose", verbose);
	image_in->load(image_filename);

	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);
	}

	ipCore *dct2D = NULL;

	dct2D = new ipDCT2D(block_size, "gray", dimDCT);
	dct2D->setBOption("verbose", verbose);

	ipBlock *convolution = NULL;

	convolution = new ipBlock(image_in->width, image_in->height, "gray", dct2D, dimDCT, block_size, block_overlap);
	convolution->setBOption("verbose", verbose);
	convolution->process(image_in);

	int convo_output_size = dimDCT;	
	for(int i = 0 ; i < convolution->n_block ; i++)
	{
	   	real *blocks_ = convolution->seq_out->frames[i];
	
	   	if(printblocks)
		{
		   	real energy = 0.0;
			real log_energy;

			print("%d : [ ", i);
			for(int j = 0 ; j < convo_output_size ; j++)
			{
				print("%g ", blocks_[j]);
				energy += blocks_[j] * blocks_[j];
			}
			energy = sqrt(energy);
			log_energy = energy < 1.0 ? 0.0 : (real) log(energy); 
			print("] (log E=%g)\n", log_energy);
		}
	}

	if(inverse)
	{
		ipDCT2Dinverse *idct = new ipDCT2Dinverse(dimDCT, block_size);
		
		print("ipDCT (dctdim = %d, output_size=%d) ...\n", dimDCT, convo_output_size);
	
		real *image_ = new real [image_in->width * image_in->height];
		int *contributions_ = new int [image_in->width * image_in->height];

		for(int i = 0 ; i < image_in->width * image_in->height ; i++)
		{
			image_[i] = 0.0;
			contributions_[i] = 0;
		}


		Sequence *seq_tmp = new Sequence(1, convolution->seq_out->frame_size);
		
		for(int rowblock=0; rowblock<convolution->rowblocks; rowblock++)
		{
	   		for(int colblock=0; colblock<convolution->colblocks; colblock++) 
			{
				int row = convolution->row_offset + rowblock*convolution->pixeladv;
				int col = convolution->col_offset + colblock*convolution->pixeladv;

				// Storing blocks as [colblock][rowblock] (i.e. [x][y])
				int index_block = colblock + rowblock * convolution->colblocks;
				seq_tmp->copyFrom(convolution->seq_out->frames[index_block]);
				
				idct->process(seq_tmp);
   		
				real *idct_output = idct->seq_out->frames[0];

				int idx_idct_ = 0;
				for(int y = row ; y < row + block_size ; y++)
				{
				   	int bias_ = y * image_in->width;
					
					for(int x = col ; x < col + block_size ; x++)
					{
					   	int idx_ = bias_ + x;

						image_[idx_] += idct_output[idx_idct_++];
						contributions_[idx_]++;
					}
				}
			}
		}	

		delete seq_tmp;

		for(int i = 0 ; i < image_in->width * image_in->height ; i++)
			if(contributions_[i] == 0.0) 
			{
			   	warning("No contributions for that pixel.");

			   	image_[i] = 0.0;
			}
			else 
			{
			   	if(image_[i] > 255.0)
				{
			   		image_[i] = 255.0;
				}
				else if(image_[i] < 0.0)
				{
			   		image_[i] = 0.0;
				}
					
			   	image_[i] /= (float) contributions_[i];
			}
			
		ImageGray *grayimage_= new ImageGray(image_in->width, image_in->height);
	
		grayimage_->copyFrom(image_in->width, image_in->height, image_, "gray");
		grayimage_->save("inverse.pgm");
		
		delete grayimage_;
		delete [] contributions_;
		delete [] image_;
		delete idct;
	}

	delete convolution;
	delete dct2D;
	delete image_in;

	return(0);
}