[HNOI 2013] 切糕

题目描述 Description

经过千辛万苦小 A 得到了一块切糕，切糕的形状是长方体，小 A 打算拦腰将切糕切成两半分给小 B。出于美观考虑，小 A 希望切面能尽量光滑且和谐。于是她找到你，希望你能帮她找出最好的切割方案。 
出于简便考虑，我们将切糕视作一个长 P、宽 Q、高 R 的长方体点阵。我们将位于第 z层中第 x 行、第 y 列上(1≤x≤P, 1≤y≤Q, 1≤z≤R)的点称为(x,y,z)，它有一个非负的不和谐值 v(x,y,z)。一个合法的切面满足以下两个条件： 
1. 与每个纵轴(一共有 P*Q 个纵轴)有且仅有一个交点。即切面是一个函数 f(x,y)，对于所有 1≤x≤P, 1≤y≤Q,我们需指定一个切割点 f(x,y),且 1≤f(x,y)≤R。 
2. 切面需要满足一定的光滑性要求，即相邻纵轴上的切割点不能相距太远。对于所有的 1≤x,x’≤P 和 1≤y,y’ ≤Q，若|x-x’|+|y-y’|=1，则|f(x,y)-f(x’,y’)| ≤D，其中 D 是给定的一个非负整数。 
可能有许多切面f 满足上面的条件，小A 希望找出总的切割点上的不和谐值最小的那个，即 ∑v(x,y, f(x,y))最小。

输入描述 Input Description

输入文件第一行是三个正整数P,Q,R，表示切糕的长P、宽Q、高R。第二行有一个非负整数D，表示光滑性要求。接下来是R个P行Q列的矩阵，第z个矩阵的第x行第y列是v(x,y,z) (1≤x≤P, 1≤y≤Q, 1≤z≤R)。

输出描述 Output Description

输出仅包含一个整数，表示在合法基础上最小的总不和谐值。

样例输入 Sample Input

input1
2 2 2 
1 
6 1 
6 1 
2 6 
2 6 
input2 
2 2 2 
0 
5 1 
5 1 
2 5 
2 5

样例输出 Sample Output

output1
6 
output2
12

数据范围及提示 Data Size & Hint

100%的数据满足P,Q,R≤40，0≤D≤R，且给出的所有的不和谐值不超过1000。 

简化一下这道题目的意思，给定一个P*Q*R的点阵，每一个点均有一个权值，要求在每一竖列的R个点中选出一个点，使得这P*Q个点的权值最小，但是题目有一个限制条件，给定一个参数D，要求相邻两竖列选出的点的高度差不能超过D。

很明显，这一道题是要用最小割解决的，被割掉的权值对应的点就是要选出的点。但是由于D的限制，所以有些点是不允许一起被选的，如图（方框中是权值）：

对于D的限制，我们不难想到，可以从  [x][y][z]  向  [x][y][z-d]的能扩展的四个方向连接一条容量为无穷大的边。（下图中只连一条）：

建图方法应该是这样的：

我们构造一个P*Q*(R+1)的点阵，用(i,j,k)表示一个点
那么(i,j,k) -> (i,j,k+1) 的流量为原图中(i,j,k)的不和谐值。
S向所有底层连无穷边，所有顶层向T连无穷边。
那么(i,j,k) -> (i,j,k+1)被割表示f(i,j)=k。
同时我们让(i,j,k) -> (i’,j’,k-D)连一条无穷边，就能保证相邻两个不超过D了。

所以代码就很容易了，一个最大流而已：

#include <cstdio>
#include <vector>
#include <algorithm>
#include <cstring>
#define clean(x,y) memset(x,y,sizeof(x))
#define pb push_back
#define MaxN 65000
#define MaxL 45
using namespace std;
const int
	dx[]={1,-1,0,0},
	dy[]={0,0,1,-1}, 
	INF=~0U>>1;
int lab[MaxL][MaxL][MaxL];
int head[MaxN],_head[MaxN],h[MaxN],v[MaxN],tot;
int p,q,r,w,d,s,t,num;
struct edge
{
	int v,cap,next;
	edge(){}
	inline void put(int x,int y,int c)
	{
		v=y,cap=c,next=head[x];
	}
}a[MaxN*6];
inline void AddEdge(int x,int y,int c)
{
	a[tot].put(x,y,c);
	head[x]=tot++;
	a[tot].put(y,x,0);
	head[y]=tot++;
}
int sap(int x,int lim)
{
	int y,flow,sum=0;
	if(x==t) return lim;
	for(int i=head[x];~i;i=a[i].next)
	{
		y=a[i].v;
		if((h[y]+1==h[x])&&a[i].cap)
			head[x]=i,
			flow=sap(y,min(a[i].cap,lim-sum)),
			a[i].cap-=flow,
			a[i^1].cap+=flow,
			sum+=flow;
		if(sum==lim) return sum;
	}
	if(h[s]>num) return sum;
	if(!--v[h[x]]) h[s]=num;
	v[++h[x]]++;
	head[x]=_head[x];
	return sum;
}
inline void fix(int x,int y,int z)
{
	int _x,_y;
	for(int k=0;k<4;k++)
	{
		_x=x+dx[k];
		_y=y+dy[k];
		if(_x>0&&_x<=p&&_y>0&&_y<=q)
			AddEdge(lab[x][y][z],lab[_x][_y][z-d],INF);
	}
}
inline void init()
{
	clean(head,-1);
	scanf("%d%d%d%d",&p,&q,&r,&d);
	for(int k=1;k<=r+1;k++)
		for(int i=1;i<=p;i++)
			for(int j=1;j<=q;j++)
				lab[i][j][k]=++num;
	s=0,t=++num,v[0]=++num;
	for(int x=1;x<=p;x++)
		for(int y=1;y<=q;y++)
			AddEdge(s,lab[x][y][1],INF),
			AddEdge(lab[x][y][r+1],t,INF);
	for(int z=1;z<=r;z++)
		for(int x=1;x<=p;x++)
			for(int y=1;y<=q;y++)
				scanf("%d",&w),
				AddEdge(lab[x][y][z],lab[x][y][z+1],w);
	for(int z=r+1;z>d;z--)
		for(int x=1;x<=p;x++)
			for(int y=1;y<=q;y++)
				fix(x,y,z);
}
inline void work()
{
	int ans=0;
	memcpy(_head,head,sizeof(head));
	while(h[s]<num)
		ans+=sap(s,INF);
	printf("%d",ans);
}
int main()
{
	init();
	work();
	return 0;
}