Rectangling Panoramic Images via Warping论文的matlab代码

#include 
#include 
#include 
#include 
#include 
#include 
using namespace std;

unsigned char* img;
int w h;
int* mask;
int* outMask;
unsigned char* outImg;
int* xDispMap;
int* yDispMap;
int* outXDispMap;
int* outYDispMap;

void readPPM（const char* name unsigned char** img int* w int* h） {
	FILE* file = fopen（name “rb“）;
	fscanf（file “P6
%d %d
255
“ w h）;
	int size = （*w） * （*h） * 3;
	*img = new unsigned char[size];
	fread（*img sizeof（unsigned char） size file）;
	fclose（file）;
}
void writePPM（const char* name const unsigned char* img int w int h） {
	FILE* file = fopen（name “wb“）;
	fprintf（file “P6
%d %d
255
“ w h）;
	int size = w * h * 3;
	fwrite（img sizeof（unsigned char） size file）;
	fclose（file）;
}
void writePGM（const char* name const unsigned char* img int w int h） {
	FILE* file = fopen（name “wb“）;
	fprintf（file “P5
%d %d
255
“ w h）;
	int size = w * h;
	fwrite（img sizeof（unsigned char） size file）;
	fclose（file）;
}
void energy（const unsigned char* image int W int H int* energy int* mask）{
	memset（energy 0 sizeof（energy[0]）*W*H）;
	//int addValue = max（WH）*3*10;
	int addValue = 1000000;
	#pragma omp parallel for
	for（int i=0; i		for（int j=0; j			for（int k=0; k<3; ++k）
				energy[i*W+j] += abs（（int）image[（i*W+（j+1））*3+k] - （int）image[（i*W+j）*3+k]）;
	#pragma omp parallel for
	for（int i=0; i		for（int j=0; j			for（int k=0; k<3; ++k）
				energy[i*W+j] += abs（（int）image[（（i+1）*W+j）*3+k] - （int）image[（i*W+j）*3+k]）;
	#pragma omp parallel for
	for（int i=0; i		for（int j=0; j			energy[i*W+j] += mask[i*W+j] * addValue;
}
///////////////////////Vertical Seam////////////////////////////
void accumulateVerEn（int * energy int W int H int* verticalEnergy int* verticalRoute）{
	/////////Initialize////////////
	memset（verticalEnergy 0 sizeof（verticalEnergy[0]）*W*H）;
	memset（verticalRoute 0 sizeof（verticalRoute[0]）*W*H）;
	#pragma omp parallel for
	for（int i=0; i		verticalEnergy[i] = energy[i];

	////////Accumulte the min energy and find the route//////////////
	for（int i=1; i		for（int j=0; j			if（j-1<0）{
				verticalEnergy[i*W+j] = energy[i*W+j] + min（verticalEnergy[（i-1）*W+j] verticalEnergy[（i-1）*W+j+1]）;//e（ij） + min（ M（i-1j）  M（i-1 j+1））
				verticalRoute[i*W+j] = （verticalEnergy[（i-1）*W+j]<=verticalEnergy[（i-1）*W+j+1]） ? 0 : 1;
			}
			else if（j+1>=W）{
				verticalEnergy[i*W+j] = energy[i*W+j] + min（verticalEnergy[（i-1）*W+j-1] verticalEnergy[（i-1）*W+j]）;//e（ij） + min（ M（i-1j-1）  M（i-1 j））
				verticalRoute[i*W+j] = （verticalEnergy[（i-1）*W+j-1]			}
			else{
				verticalEnergy[i*W+j] = energy[i*W+j] + min（min（verticalEnergy[（i-1）*W+j-1] verticalEnergy[（i-1）*W+j]） verticalEnergy[（i-1）*W+j+1]）;//e（ij） + min（ M（i-1j-1）  M（i-1 j） M（i-1 j+1））
				verticalRoute[i*W+j] = （verticalEnergy[（i-1）*W+j-1]