An Automatic Square Jigsaw Solver

jigsaw_solver.cpp

#include <opencv.hpp>
#include <vector>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <string>
#include <iostream>

using namespace cv;
using namespace std;

int frame_no = 0;
/* 分块信息 */
const int square_len = 28;
const int jigsaw_rows = 20;     // 分块以后jigsaw的大小信息 横向27块，纵向20块
const int jigsaw_cols = 27;
const int total_parts = 540;    // 总共有多少块

vector<Rect> part_pos;  // 打乱以后各块的位置
vector<Rect> orig_pos;  // 分块以后各块的位置信息
vector<Rect> dst_pos;  // 结果图像中的各块位置信息 是源图像的4倍大小 免去了超界的麻烦


bool part_used[total_parts]; // 这个part使用过没有
int cur_part_count;  // 还剩多少个可用的part

/* 相似性度量信息 */
double confidence[total_parts][total_parts][4];		// 0 1 2 3分别代表第二块在第一块的上 左 下 右方
int best_neighbors[total_parts][4];			// 每一块(part)它上下左右和它最佳匹配的part的编号
double max_confidence = -99999999, min_confidence = 99999999;

/* 邻接与网格占用信息 */
struct neighbor
{
	Point pos;		// 位置
	int adjacent_count;  // 邻接的已填充的block的数量
};
vector<neighbor> neighbors; // 当前的邻接空隙节点集
int block_part_idx[jigsaw_rows*2][jigsaw_cols*2]; // block被哪一个part占用 
int old_block_part_idx[jigsaw_rows*2][jigsaw_cols*2]; 
bool is_visited[jigsaw_rows*2][jigsaw_cols*2];  // segmentation阶段是否被访问到
int xmin, xmax, ymin, ymax; // 当前填充区域的包围盒


/* 聚类信息 */
int seg_labels[jigsaw_rows*2][jigsaw_cols*2];   // 分割后各块的编号
struct segment
{
	vector<Point> blocks;   // 属于这一类的块的集合
	Point upleft, downright;   // 包围盒
};

vector<segment> segments;
int biggest_seg_idx;   // 最大的一类在segments中的下标
int old_biggest_seg_size;  // 上一次最大的一类的大小

// 平方
double sqr(double a)
{
	return a*a;
}

// 比较函数
bool cmp(neighbor a, neighbor b)
{
	return a.adjacent_count > b.adjacent_count;
}

// 数字转字符串 a>0
string int2str(int a)
{
	string ans = "";
	if (a == 0)
		return "0";
	while (a>0)
	{
		ans += a % 10 + '0';
		a /= 10;
	}
	reverse(ans.begin(), ans.end());
	return ans;
}

// 拷贝part 将第一幅图的一部分拷入第二幅图的指定位置
void copy_part(Mat img1, Rect roi1, Mat &img2, Rect roi2)
{
	Mat buf_img1(img1, roi1);
	Mat buf_img2(img2, roi2);

	buf_img1.copyTo(buf_img2);	
}

// 生成jigsaw
void generate_jigsaw(Mat src, Mat &jigsaw)
{
	int i,j;
		
	for (i = 0; i < src.rows; i+=square_len)
		for (j = 0; j < src.cols; j+=square_len)
		{
			part_pos.push_back(Rect(j, i, square_len, square_len));
			orig_pos.push_back(Rect(j, i, square_len, square_len));			
		}
	random_shuffle(part_pos.begin(), part_pos.end());
	
	for (i = 0; i < part_pos.size(); i++)
	{
		Mat buf_src(src, part_pos[i]);
		Mat buf_jigsaw(jigsaw, orig_pos[i]);
		buf_src.copyTo(buf_jigsaw);
	}			

	printf("%d rows and %d cols\n", src.rows / square_len, src.cols / square_len);
}

// 计算两个parts的一个相对位置的confidence, t=1 square error, t=2 Lpq, t=3 prediction 
// 左表示part2在part1的左边 以此类推
double calc_confidence(Mat jigsaw, int idx1, int idx2, int pos, int t)
{
	Mat part1(jigsaw, orig_pos[idx1]);
	Mat part2(jigsaw, orig_pos[idx2]);

	double conf = 0;
	int i,j;
	switch (pos)
	{
	// 上
	case 0:
		{
			for (i = 0; i < square_len; i++)
			{
				if (t == 1)
				{
					conf += sqr(part2.at<Vec3b>(square_len-1, i)[0] - part1.at<Vec3b>(0, i)[0]);
					conf += sqr(part2.at<Vec3b>(square_len-1, i)[1] - part1.at<Vec3b>(0, i)[1]);
					conf += sqr(part2.at<Vec3b>(square_len-1, i)[2] - part1.at<Vec3b>(0, i)[2]);
				}
				else if (t == 2)
				{
					conf += pow(abs(part2.at<Vec3b>(square_len-1, i)[0] - part1.at<Vec3b>(0, i)[0]), 3/10.0);
					conf += pow(abs(part2.at<Vec3b>(square_len-1, i)[1] - part1.at<Vec3b>(0, i)[1]), 3/10.0);
					conf += pow(abs(part2.at<Vec3b>(square_len-1, i)[2] - part1.at<Vec3b>(0, i)[2]), 3/10.0);				
				}
				else
				{
					conf += pow(pow(abs(2.0*part1.at<Vec3b>(0, i)[0] - part1.at<Vec3b>(1, i)[0] - part2.at<Vec3b>(square_len-1, i)[0]), 3/10.0) +
							pow(abs(2.0*part2.at<Vec3b>(square_len-1, i)[0] - part2.at<Vec3b>(square_len-2, i)[0] - part1.at<Vec3b>(0, i)[0]), 3/10.0), 5/24.0);

					conf += pow(pow(abs(2.0*part1.at<Vec3b>(0, i)[1] - part1.at<Vec3b>(1, i)[1] - part2.at<Vec3b>(square_len-1, i)[1]), 3/10.0) +
						pow(abs(2.0*part2.at<Vec3b>(square_len-1, i)[1] - part2.at<Vec3b>(square_len-2, i)[1] - part1.at<Vec3b>(0, i)[1]), 3/10.0), 5/24.0);

					conf += pow(pow(abs(2.0*part1.at<Vec3b>(0, i)[2] - part1.at<Vec3b>(1, i)[2] - part2.at<Vec3b>(square_len-1, i)[2]), 3/10.0) +
						pow(abs(2.0*part2.at<Vec3b>(square_len-1, i)[2] - part2.at<Vec3b>(square_len-2, i)[2] - part1.at<Vec3b>(0, i)[2]), 3/10.0), 5/24.0);
				}
			}
			break;
		}
	// 左
	case 1:
		{
			for (i = 0; i < square_len; i++)
			{
				if (t == 1)
				{
					conf += sqr(part2.at<Vec3b>(i, square_len-1)[0] - part1.at<Vec3b>(i, 0)[0]);
					conf += sqr(part2.at<Vec3b>(i, square_len-1)[1] - part1.at<Vec3b>(i, 0)[1]);
					conf += sqr(part2.at<Vec3b>(i, square_len-1)[2] - part1.at<Vec3b>(i, 0)[2]);
				}
				else if (t == 2)
				{
					conf += pow(abs(part2.at<Vec3b>(i, square_len-1)[0] - part1.at<Vec3b>(i, 0)[0]), 3/10.0);
					conf += pow(abs(part2.at<Vec3b>(i, square_len-1)[1] - part1.at<Vec3b>(i, 0)[1]), 3/10.0);
					conf += pow(abs(part2.at<Vec3b>(i, square_len-1)[2] - part1.at<Vec3b>(i, 0)[2]), 3/10.0);
				}
				else
				{
					conf += pow(pow(abs(2.0*part1.at<Vec3b>(i, 0)[0] - part1.at<Vec3b>(i, 1)[0] - part2.at<Vec3b>(i, square_len-1)[0]), 3/10.0) +
						pow(abs(2.0*part2.at<Vec3b>(i, square_len-1)[0] - part2.at<Vec3b>(i, square_len-2)[0] - part1.at<Vec3b>(i, 0)[0]), 3/10.0), 5/24.0);

					conf += pow(pow(abs(2.0*part1.at<Vec3b>(i, 0)[1] - part1.at<Vec3b>(i, 1)[1] - part2.at<Vec3b>(i, square_len-1)[1]), 3/10.0) +
						pow(abs(2.0*part2.at<Vec3b>(i, square_len-1)[1] - part2.at<Vec3b>(i, square_len-2)[1] - part1.at<Vec3b>(i, 0)[1]), 3/10.0), 5/24.0);

					conf += pow(pow(abs(2.0*part1.at<Vec3b>(i, 0)[2] - part1.at<Vec3b>(i, 1)[2] - part2.at<Vec3b>(i, square_len-1)[2]), 3/10.0) +
						pow(abs(2.0*part2.at<Vec3b>(i, square_len-1)[2] - part2.at<Vec3b>(i, square_len-2)[2] - part1.at<Vec3b>(i, 0)[2]), 3/10.0), 5/24.0);
				}
			}
			break;
		}
	// 下
	case 2:
		{
			for (i = 0; i < square_len; i++)
			{
				if (t == 1)
				{
					conf += sqr(part2.at<Vec3b>(0, i)[0] - part1.at<Vec3b>(square_len-1, i)[0]);
					conf += sqr(part2.at<Vec3b>(0, i)[1] - part1.at<Vec3b>(square_len-1, i)[1]);
					conf += sqr(part2.at<Vec3b>(0, i)[2] - part1.at<Vec3b>(square_len-1, i)[2]);
				}
				else if (t == 2)
				{
					conf += pow(abs(part2.at<Vec3b>(0, i)[0] - part1.at<Vec3b>(square_len-1, i)[0]), 3/10.0);
					conf += pow(abs(part2.at<Vec3b>(0, i)[1] - part1.at<Vec3b>(square_len-1, i)[1]), 3/10.0);
					conf += pow(abs(part2.at<Vec3b>(0, i)[2] - part1.at<Vec3b>(square_len-1, i)[2]), 3/10.0);
				}
				else
				{
					conf += pow(pow(abs(2.0*part1.at<Vec3b>(square_len-1, i)[0] - part1.at<Vec3b>(square_len-2, i)[0] - part2.at<Vec3b>(0, i)[0]), 3/10.0) +
						pow(abs(2.0*part2.at<Vec3b>(0, i)[0] - part2.at<Vec3b>(1, i)[0] - part1.at<Vec3b>(square_len-1, i)[0]), 3/10.0), 5/24.0);

					conf += pow(pow(abs(2.0*part1.at<Vec3b>(square_len-1, i)[1] - part1.at<Vec3b>(square_len-2, i)[1] - part2.at<Vec3b>(0, i)[1]), 3/10.0) +
						pow(abs(2.0*part2.at<Vec3b>(0, i)[1] - part2.at<Vec3b>(1, i)[1] - part1.at<Vec3b>(square_len-1, i)[1]), 3/10.0), 5/24.0);

					conf += pow(pow(abs(2.0*part1.at<Vec3b>(square_len-1, i)[2] - part1.at<Vec3b>(square_len-2, i)[2] - part2.at<Vec3b>(0, i)[2]), 3/10.0) +
						pow(abs(2.0*part2.at<Vec3b>(0, i)[2] - part2.at<Vec3b>(1, i)[2] - part1.at<Vec3b>(square_len-1, i)[2]), 3/10.0), 5/24.0);
				}
			}
			break;
		}
	// 右
	case 3:
		{
			for (i = 0; i < square_len; i++)
			{
				if (t == 1)
				{
					conf += sqr(part2.at<Vec3b>(i, 0)[0] - part1.at<Vec3b>(i, square_len-1)[0]);
					conf += sqr(part2.at<Vec3b>(i, 0)[1] - part1.at<Vec3b>(i, square_len-1)[1]);
					conf += sqr(part2.at<Vec3b>(i, 0)[2] - part1.at<Vec3b>(i, square_len-1)[2]);
				}
				else if (t == 2)
				{
					conf += pow(abs(part2.at<Vec3b>(i, 0)[0] - part1.at<Vec3b>(i, square_len-1)[0]), 3/10.0);
					conf += pow(abs(part2.at<Vec3b>(i, 0)[1] - part1.at<Vec3b>(i, square_len-1)[1]), 3/10.0);
					conf += pow(abs(part2.at<Vec3b>(i, 0)[2] - part1.at<Vec3b>(i, square_len-1)[2]), 3/10.0);
				}
				else
				{
					double tmp1 = pow(double(2.0*part1.at<Vec3b>(i, square_len-1)[0] - part1.at<Vec3b>(i, square_len-2)[0] - part2.at<Vec3b>(i, 0)[0]), 3/10.0);
					double tmp2 = pow(double(2.0*part2.at<Vec3b>(i, 0)[0] - part2.at<Vec3b>(i, 1)[0] - part1.at<Vec3b>(i, square_len-1)[0]), 3/10.0);
					double tmp3 = pow(tmp1 + tmp2, 5/24.0 );
					double tmp4 = double(2.0*part2.at<Vec3b>(i, 0)[0] - part2.at<Vec3b>(i, 1)[0] - part1.at<Vec3b>(i, square_len-1)[0]);


					conf += pow(pow(abs(2.0*part1.at<Vec3b>(i, square_len-1)[0] - part1.at<Vec3b>(i, square_len-2)[0] - part2.at<Vec3b>(i, 0)[0]), 3/10.0) +
						pow(abs(2.0*part2.at<Vec3b>(i, 0)[0] - part2.at<Vec3b>(i, 1)[0] - part1.at<Vec3b>(i, square_len-1)[0]), 3/10.0), 5/24.0);

					conf += pow(pow(abs(2.0*part1.at<Vec3b>(i, square_len-1)[1] - part1.at<Vec3b>(i, square_len-2)[1] - part2.at<Vec3b>(i, 0)[1]), 3/10.0) +
						pow(abs(2.0*part2.at<Vec3b>(i, 0)[1] - part2.at<Vec3b>(i, 1)[1] - part1.at<Vec3b>(i, square_len-1)[1]), 3/10.0), 5/24.0);

					conf += pow(pow(abs(2.0*part1.at<Vec3b>(i, square_len-1)[2] - part1.at<Vec3b>(i, square_len-2)[2] - part2.at<Vec3b>(i, 0)[2]), 3/10.0) +
						pow(abs(2.0*part2.at<Vec3b>(i, 0)[2] - part2.at<Vec3b>(i, 1)[2] - part1.at<Vec3b>(i, square_len-1)[2]), 3/10.0), 5/24.0);
				}
			}
			break;
		}
	default:break;
	}
	
	if (t == 2)
		conf =  pow(conf, 5/24.0);

	return conf;
}

// 计算所有的confidence和bestneighbors
void init_confidence_bestneighbors(Mat jigsaw)
{
	int i, j, k;
	// 利用生成好的jigsaw来计算 每块的位置在orig_pos中存储
	for (i = 0; i < total_parts; i++)
		for (j = 0; j < total_parts; j++)
		{
			if (i == j)
				confidence[i][j][0] = confidence[i][j][1] = confidence[i][j][2] = confidence[i][j][3] = 0;

			for (k = 0; k < 4; k++)
			{
				confidence[i][j][k] = calc_confidence(jigsaw, i, j, k, 2);
				if (max_confidence < confidence[i][j][k])
					max_confidence = confidence[i][j][k];
				if (min_confidence > confidence[i][j][k])
					min_confidence = confidence[i][j][k];
			}
		}
	// 归一化confidence
	for (i = 0; i < total_parts; i++)
		for (j = 0; j < total_parts; j++)
			for (k = 0; k < 4; k++)
			{
				confidence[i][j][k] = (max_confidence - confidence[i][j][k]) / max_confidence;
				// confidence[i][j][k] = exp(- confidence[i][j][k]);
			}

	// 初始化best_neighbors
	for (i = 0; i < total_parts; i++)
	{
		for (k = 0; k < 4; k++)
		{
			double best_conf = -9999;
			double best_neighbor_idx = -1;
			for (j = 0; j < total_parts; j++)
				if (confidence[i][j][k] > best_conf)
				{
					best_conf = confidence[i][j][k];
					best_neighbor_idx = j;
				}
			best_neighbors[i][k] = best_neighbor_idx;
		}
	}
}


// 两个part i和j在 k这个相对位置关系下是否为bestbuddy
bool is_bestbuddy(int i, int j, int k)
{
	int opposite = (k + 2) % 4;  // j和i的相对位置关系

	return best_neighbors[i][k] == j && best_neighbors[j][opposite] == i;
}

void add_a_neighbor(int x, int y)
{
	int i, idx = -1;
	for (i = 0; i < neighbors.size(); i++)
		if (neighbors[i].pos.x == x && neighbors[i].pos.y == y)
			idx = i;

	if (idx >= 0)
	{
		neighbors[idx].adjacent_count++;
	}
	else
	{
		struct neighbor tmp;
		tmp.pos.x = x;
		tmp.pos.y = y;
		tmp.adjacent_count = 1;
		neighbors.push_back(tmp);
	}
}

void delete_a_neighbor(int neighbor_idx)
{
	neighbors.erase(neighbors.begin() + neighbor_idx);
}

bool check_bounds(int x, int y)
{
	if (abs(x - xmin) >= jigsaw_cols || abs(x - xmax) >= jigsaw_cols || abs(y - ymin) >= jigsaw_rows || abs(y - ymax) >= jigsaw_rows)
		return false;
	else
		return true;
}

void update_bounds(int x, int y)
{
	if (x < xmin)
		xmin = x;
	if (x > xmax)
		xmax = x;
	if (y < ymin)
		ymin = y;
	if (y > ymax)
		ymax = y;
}

// 添加neighbor blocks 需要考虑边界条件
void add_neighbors(int x, int y)
{
	if (check_bounds(x-1, y) && block_part_idx[y][x-1] < 0)
	{
		add_a_neighbor(x-1, y);		
		update_bounds(x, y);
	}

	if (check_bounds(x+1, y) && block_part_idx[y][x+1] < 0)
	{
		add_a_neighbor(x+1, y);		
		update_bounds(x, y);
	}
	if (check_bounds(x, y-1) && block_part_idx[y-1][x] < 0)
	{
		add_a_neighbor(x, y-1);
		update_bounds(x, y);
	}		
	if (check_bounds(x, y+1) && block_part_idx[y+1][x] < 0)
	{
		add_a_neighbor(x, y+1);	
		update_bounds(x, y);
	}
}


// 计算对于一个neighbor block，与它配对的最佳part是哪个，并返回其confidence
double calc_max_confidence(int neighbor_idx, int &part_idx)
{
	int i, j;
	Point pos = neighbors[neighbor_idx].pos;
	int count = neighbors[neighbor_idx].adjacent_count;
	double max_conf = -9999;
	int dir[4][2] = {{0, -1}, {-1, 0}, {0, 1}, {1, 0}}; // 第二块在当前块的上 左 下 右

	for (i = 0; i < total_parts; i++)  // i为当前块的编号
	{
		if (part_used[i])
			continue;

		double conf = 0;
		for (j = 0; j < 4; j++)
		{
			int idx = block_part_idx[pos.y + dir[j][1]][pos.x + dir[j][0]]; // idx为第二块的编号
			if (idx >= 0) // 被占用才计算
				conf += confidence[i][idx][j];   // j为相对位置信息
		}

		conf /= count;
		if (conf > max_conf)
		{
			max_conf = conf;
			part_idx = i;
		}
	}

	return max_conf;
}

// 选择part 与 block之间最佳的一个组合(平均的confidence最大,但需要优先选择邻接已占用块最多的)
// (part_idx, x, y)是part的编号和block坐标的组合，返回block对应的neighbor编号，便于删除neighbor
int find_best_partblock(int &part_idx, int &x, int &y)
{
	int i;
	// 排序，最大邻接占用块的靠前
	sort(neighbors, cmp);
	int max_adjacent_count = neighbors[0].adjacent_count;
	double max_conf = -9999;
	int neighbor_idx;

	for (i = 0; i < neighbors.size() && neighbors[i].adjacent_count == max_adjacent_count; i++)
	{
		int idx;
		double conf = calc_max_confidence(i, idx);
		if (check_bounds(neighbors[i].pos.x, neighbors[i].pos.y) && conf > max_conf)
		{
			max_conf = conf;
			neighbor_idx = i;
			x = neighbors[i].pos.x;
			y = neighbors[i].pos.y;
			part_idx = idx;
		}
	}

	return neighbor_idx;
}

int find_best_seed()
{
	int i, j;
	int max_best_buddies = -1;
	int best_seed = -1;

	for (i = 0; i < total_parts; i++)
	{
		int best_buddies = 0;
		for (j = 0; j < 4; j++)
		{
			int k = best_neighbors[i][j];
			if (is_bestbuddy(i,k,j))
				best_buddies++;			
		}

		if (best_buddies > max_best_buddies)
		{
			max_best_buddies = best_buddies;
			best_seed = i;
		}
	}

	return best_seed;
}

// 放置阶段
void placing(Mat jigsaw, Mat &dst)
{
	int x, y, idx;
	Rect roi;
	roi.x = dst.cols / 4;
	roi.y = dst.rows / 4;
	roi.width = dst.cols / 2;
	roi.height = dst.rows / 2;

	// 放置其它的
	while(cur_part_count > 0)
	{
		// char ch = waitKey(0);
		int part_idx, neighbor_idx;
		neighbor_idx = find_best_partblock(part_idx, x, y);
		idx = y * jigsaw_cols * 2 + x;
		copy_part(jigsaw, orig_pos[part_idx], dst, dst_pos[idx]);

		block_part_idx[y][x] = part_idx;
		part_used[part_idx] = true;
		cur_part_count--;
		
		delete_a_neighbor(neighbor_idx);		
		add_neighbors(x, y);	

		imshow("dst", dst);
		frame_no++;	
		/*if (frame_no < 540)
			imwrite(("res/"+int2str(frame_no)+".jpg").c_str(), dst);*/

		waitKey(1);                 
	}
}

// 深度优先搜索
void dfs(int x, int y, int label, segment &tmp)
{
	int part_idx = block_part_idx[y][x];
	int idx_up = block_part_idx[y-1][x], idx_left = block_part_idx[y][x-1], idx_down = block_part_idx[y+1][x], idx_right = block_part_idx[y][x+1];
	
	is_visited[y][x] = true;			// 这块已经被访问过了
	seg_labels[y][x] = label;		    // 聚类的编号
	tmp.blocks.push_back(Point(x, y));   // 保存当前的block到相应的类别中
	// 更新当前的包围盒
	if (y > tmp.downright.y)
		tmp.downright.y = y;
	if (y < tmp.upleft.y)
		tmp.upleft.y = y;
	if (x > tmp.downright.x)
		tmp.downright.x = x;
	if (x < tmp.upleft.x)
		tmp.upleft.x = x;

	if (!is_visited[y-1][x] && idx_up >= 0 && is_bestbuddy(part_idx, idx_up, 0))
		dfs(x, y-1, label, tmp);
	
	if (!is_visited[y][x-1] && idx_left >= 0 && is_bestbuddy(part_idx, idx_left, 1))
		dfs(x-1, y, label, tmp);
	
	if (!is_visited[y+1][x] && idx_down >= 0 && is_bestbuddy(part_idx, idx_down, 2))
		dfs(x, y+1, label, tmp);
	
	if (!is_visited[y][x+1] && idx_right >= 0 && is_bestbuddy(part_idx, idx_right, 3))
		dfs(x+1, y, label, tmp);
	
}

// 聚类阶段
int segmentation()
{
	// 聚类前的初始化
	memset(is_visited, false, sizeof(is_visited));
	memset(seg_labels, 0, sizeof(seg_labels));
	segments.clear();

	int y, x;
	int label = 1;
	int max_seg_size = -1;
	biggest_seg_idx = -1;


	for (y = ymin; y <= ymax; y++)
		for (x = xmin; x <= xmax; x++)
			if (!is_visited[y][x])
			{
				segment tmp;				 // 一个类的列表
				tmp.upleft.x = tmp.upleft.y = 9999;
				tmp.downright.x = tmp.downright.y = -9999;
				dfs(x, y, label, tmp);

				segments.push_back(tmp);    // 存入所有聚类的信息列表中
				if (int(tmp.blocks.size()) > max_seg_size)  // 这里unsigned int和int比较，要注意不要犯错
				{
					max_seg_size = tmp.blocks.size();
					biggest_seg_idx = label - 1;
				}

				label++;
			}
	
	FILE *fp = fopen("labels.txt","w+");
	for (y = ymin; y <= ymax; y++)
	{
		for (x = xmin; x <= xmax; x++)
			fprintf(fp, "%-3d ", seg_labels[y][x]);
		fprintf(fp, "\n");
	}

	fclose(fp);
	return label - 1;
}

void shifting(int iter_level, Mat jigsaw, Mat &dst)
{
	int i,j;

	// 初始化
	memset(part_used, false, sizeof(part_used));
	cur_part_count = total_parts;

	memcpy(old_block_part_idx, block_part_idx, sizeof(block_part_idx));
	for (i = 0; i < jigsaw_rows * 2; i++)			// 每一块block设置为未被占用
		for (j = 0; j < jigsaw_cols * 2; j++)
			block_part_idx[i][j] = -1;	

	dst = Mat::zeros(jigsaw.rows * 2, jigsaw.cols * 2, CV_8UC3);  // 全黑
	neighbors.clear();			// 上一次计算的neighbors要清空
	xmin = ymin = 9999;         // 重置包围盒
	xmax = ymax = -9999;       

	// 放置阶段
	if (iter_level == 0)    // 第一次迭代放置第一个
	{		
		int x, y, dst_idx, idx;
		y = jigsaw_rows; x = jigsaw_cols;
		dst_idx = y * jigsaw_cols * 2 + x;
		// idx = int(rand() * 1.0 / RAND_MAX * total_parts);   // 随机选择
		idx = 0;											    // 选第一个
		// idx = find_best_seed();									// 选择best buddy最多的那个
		copy_part(jigsaw, orig_pos[idx], dst, dst_pos[dst_idx]);

		block_part_idx[y][x] = idx;  //占用掉了这个block
		part_used[idx] = true;        //使用了这个part
		cur_part_count--;
		
		update_bounds(x, y);
		add_neighbors(x, y);

		imshow("dst", dst);			
	}
	else					// 找上次迭代最大的那一个segment放入
	{
		Point seg_center;
		Point center(jigsaw_cols, jigsaw_rows);
		Point offset;

		seg_center.x = (segments[biggest_seg_idx].downright.x + segments[biggest_seg_idx].upleft.x) / 2 ;
		seg_center.y = (segments[biggest_seg_idx].downright.y + segments[biggest_seg_idx].upleft.y) / 2;
		offset = center - seg_center;		

		int x, y, idx, dst_idx;
		for (i = 0; i < segments[biggest_seg_idx].blocks.size(); i++)
		{
			y = segments[biggest_seg_idx].blocks[i].y;
			x = segments[biggest_seg_idx].blocks[i].x;
			idx = old_block_part_idx[y][x];   // 它对应的是那个part

			y = y + offset.y;				// 把这一块移到屏幕中央的位置
			x = x + offset.x;
			dst_idx = y * jigsaw_cols * 2 + x;        // 要拷贝入的dst的位置

			copy_part(jigsaw, orig_pos[idx], dst, dst_pos[dst_idx]);

			block_part_idx[y][x] = idx;
			part_used[idx] = true; 
			cur_part_count--;
			update_bounds(x, y);			
			imshow("dst", dst);	
			// waitKey(0);
		}

		// 添加neighbors
		for (i = 0; i < segments[biggest_seg_idx].blocks.size(); i++)
		{
			y = segments[biggest_seg_idx].blocks[i].y;
			x = segments[biggest_seg_idx].blocks[i].x;
			y = y + offset.y;
			x = x + offset.x;
			add_neighbors(x, y);
		}

		waitKey(500);
	}

 	placing(jigsaw, dst);
	// printf("Bounding box: (%d, %d) (%d, %d)\n", xmin, ymin, xmax, ymax);

	int total_labels = segmentation();
	printf("iteration #%d: totally %d segments, biggest segment with %d blocks\n", iter_level+1, total_labels, segments[biggest_seg_idx].blocks.size());
	// waitKey(0);
	if (segments[biggest_seg_idx].blocks.size() > old_biggest_seg_size)   // 还有增长空间 继续
	{
		old_biggest_seg_size = segments[biggest_seg_idx].blocks.size();
		shifting(iter_level+1, jigsaw, dst);		
	}
	else																// 收敛了
		return;
}

void solve_jigsaw(Mat jigsaw, Mat &dst)
{
	int i, j;
	init_confidence_bestneighbors(jigsaw);

	 // 结果图像中各block的位置
	for (i = 0; i < dst.rows; i+=square_len)  
		for (j = 0; j < dst.cols; j+=square_len)
			dst_pos.push_back(Rect(j, i, square_len, square_len));
	
	shifting(0, jigsaw, dst);	
}

int main()
{
	string jigsaw_path = "jigsaw_pics/";
	string file_name;

	srand(unsigned(time(NULL)));

	cout << "please enter image file name:";
	cin >> file_name;
	Mat src;
	src = imread(jigsaw_path+file_name);
	if (src.cols != 756 || src.rows != 560)
	{
		cerr << "image size must be 756x560" << endl;
		return 0;
	}

	Mat jigsaw(src.rows, src.cols, CV_8UC3, Scalar(0,0,0));
	Mat dst(src.rows * 2, src.cols * 2, CV_8UC3, Scalar(0,0,0));		

	generate_jigsaw(src, jigsaw);
	imshow("src",src);
	imshow("jigsaw",jigsaw);
	solve_jigsaw(jigsaw, dst);

	// 显示最终结果
	Point ul, dr;
	ul.x = dst_pos[ymin * jigsaw_cols * 2 + xmin].x; ul.y = dst_pos[ymin * jigsaw_cols * 2 + xmin].y;
	dr.x = dst_pos[ymax * jigsaw_cols * 2 + xmax].x + dst_pos[ymax * jigsaw_cols * 2 + xmax].width; dr.y = dst_pos[ymax * jigsaw_cols * 2 + xmax].y + dst_pos[ymax * jigsaw_cols * 2 + xmax].height;
	Rect roi(ul, dr);
	imshow("dst", dst(roi));	
	
	waitKey(0);
	return 0;
}