segment.cpp

this is image segmentation code

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "graph.h"

#define OBJ	1
#define BKG	-1

unsigned char **ReadPPM(char *file,int *width, int *height);
void WritePPM(char *file, unsigned char **buf, int width, int height);
void	LabelMap(int **map, char *file, int label, int width, int height);
void	Segment(unsigned char **buf, int **map, int width, int height);
double	**FindHist(unsigned char **buf, int **map, int width, int height, int label);
double	FindB(unsigned char **buf, int x1, int y1, int x2, int y2);
double	FindR(unsigned char **buf, int x, int y,  double **hist);

main(int argc, char **argv)
{
	if(argc != 4)
	{
		printf("Usage: %s file.ppm object.txt background.txt\n",argv[0]);
		exit(0);
	}
	unsigned char **buf;
	int	width, height;
	int	**map;
	buf = ReadPPM(argv[1],&width,&height);
	map = (int **) calloc(height, sizeof(int *));
	for(int i = 0; i < height; i++)
		map[i] = (int *) calloc(width, sizeof(int));
	LabelMap(map,argv[2],OBJ,width,height);
	LabelMap(map,argv[3],BKG,width,height);
	Segment(buf,map,width,height);
	for(int i = 0; i < height; i++)
		for(int j = 0; j < width; j++)
		{
			if(map[i][j] == BKG)
			{
				buf[i][3*j] = 0;
				buf[i][3*j+1] = 0;
				buf[i][3*j+2] = 255;
			}
			else
			{
				/*
				buf[i][3*j] = 255;
				buf[i][3*j+1] = 0;
				buf[i][3*j+2] = 0;
				*/
			}
		}
	WritePPM("out.ppm",buf,width,height);
}

void	Segment(unsigned char **buf, int **map, int width, int height)
{
	Graph::node_id node[height][width];
	Graph *g = new Graph();
	for(int i = 0; i < height; i++)
		for(int j = 0; j < width; j++)
			node[i][j] = g->add_node();
	double	K = 0;
	for(int i = 0; i < height; i++)
		for(int j = 0; j < width; j++)
		{
			for(int k = -1; k < 2; k++)
				for(int l = -1; l < 2; l++)
					if(i+k >= 0 && i+k < height && j+l >= 0 && j+l < width)
					{
						double B;
						if(k == 0 && l == 0)
							continue;
						B = FindB(buf,i,j,i+k,j+l);
						g->add_edge(node[i][j],node[i+k][j+l],short(B),short(B));
//						printf("(%d,%d):(%d,%d) %f %d\n",i,j,i+k,j+l,B,short(B));
						if(1+B > K)
							K = 1+B;
					}
		}
	double	**histFg,**histBg;
	histFg = FindHist(buf, map, width, height, OBJ);
	histBg = FindHist(buf, map, width, height, BKG);
	for(int i = 0; i < height; i++)
		for(int j = 0; j < width; j++)
		{
			if(map[i][j] == OBJ)
				g->set_tweights(node[i][j],0,int(K));
			else
			if(map[i][j] == BKG)
				g->set_tweights(node[i][j],int(K),0);
			else
			{
				double Rfg,Rbg;
				unsigned char intensity;
				float temp;
				temp = (16.0*buf[i][3*j]+11.0*buf[i][3*j+1]+5.0*buf[i][3*j+2])/32.0;
				intensity = (unsigned char)temp;
				Rfg = FindR(buf,i,j,histFg);
				Rbg = FindR(buf,i,j,histBg);
				g->set_tweights(node[i][j],short(Rfg),short(Rbg));
			}
		}
	Graph::flowtype	flow = g->maxflow();
	for(int i = 0; i < height; i++)
		for(int j = 0; j < width; j++)
		{
			if(g->what_segment(node[i][j]) == Graph::SOURCE)
				map[i][j] = BKG;
			else
				map[i][j] = OBJ;
		}
}

double	FindR(unsigned char **buf, int x, int y, double **hist)
{
	double lambda = 1;
	double R;
	double	factor = 100;
	double	count = 0;
	int	r,g,b;
	r = buf[x][3*y]/10;
	g = buf[x][3*y+1]/10;
	b = buf[x][3*y+2]/10;
	for(int i = 0; i < 256; i++)
		count += hist[0][i];
	if(hist[0][r] != 0 && hist[1][g] != 0 && hist[2][b] != 0)
		R = -factor*lambda*log(hist[0][r]*hist[1][g]*hist[2][b]/pow(count,3))/log(10);
	else
		R = -factor*lambda*log(0.5/pow(count,3))/log(10);
	return R;
}

double	FindB(unsigned char **buf, int x1, int y1, int x2, int y2)
{
	double e = 2.71828;
	double	rho = 5;
	double B;
	double	Ip,Iq;
	double	factor = 10;
	
	Ip = (11.0*buf[x1][3*y1]+16.0*buf[x1][3*y1+1]+5.0*buf[x1][3*y1+2])/32.0;
	Iq = (11.0*buf[x2][3*y2]+16.0*buf[x2][3*y2+1]+5.0*buf[x2][3*y2+2])/32.0;
	B = factor*pow(e,-(Ip-Iq)*(Ip-Iq)/(2*rho*rho))/hypot(double(x1-x2),double(y1-y2));
	return B;
}

double	**FindHist(unsigned char **buf, int **map, int width, int height, int label)
{
	double	**hist;
	int	count = 0;
	hist = (double **) calloc(3, sizeof(double *));
	hist[0] = (double *) calloc(256, sizeof(double));
	hist[1] = (double *) calloc(256, sizeof(double));
	hist[2] = (double *) calloc(256, sizeof(double));
	for(int i = 0; i < height; i++)
		for(int j = 0; j < width; j++)
		{
			if(map[i][j] != label)
				continue;
		//	int	intensity;
		//	intensity = int((16.0*buf[i][3*j]+11.0*buf[i][3*j+1]+5.0*buf[i][3*j+2])/32.0);
			hist[0][buf[i][3*j]/10]++;
			hist[1][buf[i][3*j+1]/10]++;
			hist[2][buf[i][3*j+2]/10]++;
			count++;
		}
	/*
	for(int i = 0; i < 256; i++)
		hist[i] = hist[i]/double(count);
	*/
	return hist;
}

void	LabelMap(int **map, char *file, int label, int width, int height)
{
	FILE	*fp = fopen(file,"r");
	char	s[80];
	while(fgets(s,80,fp))
	{
		int	startx,starty,endx,endy;
		sscanf(s,"%d%d%d%d",&startx,&starty,&endx,&endy);
		for(int i = startx; i <= endx; i++)
			for(int j = starty; j <= endy; j++)
				if(i >=0 && i < height && j >=0 && j < width)
					map[i][j] = label;
	}
}

unsigned char **ReadPPM(char *file,int *width, int *height)
{
	unsigned char **buf;
	int i,j;
	int	Width,Height;

	FILE	*fp = fopen(file,"r");
	fscanf(fp,"P6 %d %d 255\n",width,height);
	/*
	fscanf(fp,"P6\n");
	while(1)
	{
		char s[80];
		fgets(s,80,fp);
		if(s[0] >= '0' && s[0] <= '9')
		{
			sscanf(s,"%d%d",width,height);
			break;
		}
	}
	fscanf(fp,"255\n");
	*/
	Width = *width;
	Height = *height;

	buf = (unsigned char **) calloc(Height, sizeof(unsigned char *));
	for(i = 0; i < Height; i++)
		buf[i] = (unsigned char *) calloc(Width*3,sizeof(unsigned char));
	for(i = 0; i < Height; i++)
	{
		for(j = 0; j < Width*3; j++)
		{
			fscanf(fp,"%c",&buf[i][j]);
		}
	}
	fclose(fp);
	return buf;
}

void WritePPM(char *file, unsigned char **buf, int width, int height)
{
	FILE	*fp = fopen(file,"w");
	int	i,j;

	fprintf(fp,"P6\n%d %d\n255\n",width,height);
	for(i = 0; i < height; i++)
		for(j = 0; j < 3*width; j++)
			fprintf(fp,"%c",buf[i][j]);
	fclose(fp);
}