注：下文所有的复杂度分析均认为 $n,m$ 同阶。

思路

首先可以发现，最后被淹没的形状一定是一个矩形，否则一定存在一个点有两个相邻的被淹没的点。

考虑观察这个矩形有什么性质，发现矩形一定是被一圈 1 包起来的（但是外围的四个角可以不用管），且矩形内部不能存在一行或一列全为 1。

考虑设 $rt_{i,j}$ 表示满足 $(i,j+1)\sim (i,k)$ 全部为 1 的最大的 $k$，$dn_{i,j}$ 表示满足 $(i+1,j)\sim (k,j)$ 全部为 1 的最大的 $k$。设包住矩形的左上角为 $(x_1,y_1)$，右下角为 $(x_2,y_2)$，上面的限制可以刻画为：

1. $rt_{x_1,y_1}+1\ge y_2$，加一的原因为右上角可以不为 1，后面所有加一的原因都是同理的；
2. $dn_{x_1,y_1}+1\ge x_2$；
3. $rt_{x_2,y_1}+1\ge y_2$；
4. $dn_{x_1,y_2}+1\ge x_2$；
5. 不存在一个 $i\in(x_1,x_2)$ 使得 $rt_{i,y_1}+1\ge y_2$；
6. 不存在一个 $i\in(y_1,y_2)$ 使得 $dn_{x_1,i}+1\ge x_2$。

现在固定左上角，考虑从后两个限制入手。发现若一个 $(x_2,y_2)$ 是合法的，下一个可能合法的 $y_2'$ 一定是在 $y_2$ 右边第一个满足 $dn_{x_1,y_2'}>dn_{x_1,y_2}$ 的，下一个可能合法的 $x_2'$ 同理。而这个东西可以用单调栈维护。

接下来考虑前面四条限制，对于一个 $y_2$，可能合法的 $x_2$ 一定是第一个满足 $rt_{x2,y1}+1\ge y_2$ 的，所以可以在单调栈上使用双指针找到这个 $x_2$，然后再判断其他的条件即可（此时一定满足第 3 条和第 5 条限制，前 4 条限制是好判的，第 6 条限制可以通过判断扫到的上一个 $y_2$ 是否会延伸到 $x_2$ 以下）。

分析一下每一行会扫到的 $y_2$ 的量级，发现对于 $dn_{x_1,y_2}>dn_{x_1,y_1}$ 的 $y_2$ 扫到一个后剩下的一定不会合法，$dn_{x_1,y_2}\le dn_{x_1,y_1}$ 的 $y_2$ 扫到后就会被弹出单调栈，$x_2$ 的量级也是同理的。所以合法的矩形个数是 $O(n^2)$ 的，可以被接受。

接下来考虑统计答案，对于每个点，最后形成的联通块相当于覆盖到这个点面积最小的矩形。此时做法就很多了，可以使用扫描线+标记永久化线段树，每个线段树节点维护一个可删堆即可。

时间复杂度 $O(n^2\log^2n)$，实际上根本跑不满，因为要卡满合法矩形的个数矩形的几乎不会有交，堆中元素很少，否则个数不大运算量也会变小。但是这么写有概率会被卡常，对于面积和小的跑暴力，面积和大的跑扫描线。

代码

#include<bits/stdc++.h>
using namespace std;
const int N = 5005;
int n,m,dn[N][N],rt[N][N],stk[N][N],top[N];
int l1[N*N],r1[N*N],l2[N*N],r2[N*N],s[N*N],cnt;
long long ans = 0;
char ch[N][N];
namespace sub1{
	int id[N*N],a[N][N];
	inline void main()
	{
		for(int i = 1;i<=cnt;i++) id[i] = i;
		sort(id+1,id+cnt+1,[&](int x,int y){return s[x]>s[y];});
		for(int i = 1;i<=cnt;i++)
		{
			int now = id[i];
			for(int x = l1[now];x<=r1[now];x++)
				for(int y = l2[now];y<=r2[now];y++)
					a[x][y] = s[now];
		}
		for(int i = 1;i<=n;i++)
			for(int j = 1;j<=m;j++)
				if(ch[i][j]=='0') ans+=a[i][j];
		cout<<ans;
	}
}
namespace sub2{
	vector<pair<pair<int,int>,int>> vec[N];
	priority_queue<int> q[N<<2],del[N<<2];
	#define ls (k<<1)
	#define rs (k<<1|1)
	void change(int k,int l,int r,int x,int y,int v)
	{
		if(l>y||r<x) return;
		if(l>=x&&r<=y)
		{
			if(v<0)
			{
				del[k].push(v);
				while(!q[k].empty()&&!del[k].empty()&&q[k].top()==del[k].top()) q[k].pop(),del[k].pop();
			}
			else q[k].push(-v);
			return;
		}
		int mid = (l+r)/2;
		change(ls,l,mid,x,y,v),change(rs,mid+1,r,x,y,v); 
	}
	void ask(int k,int l,int r,int x,int mn)
	{
		if(!q[k].empty()) mn = min(mn,-q[k].top());
		if(l==r)
		{
			if(ch[x][l]=='0') ans+=mn;
			return;
		}
		int mid = (l+r)/2;
		ask(ls,l,mid,x,mn),ask(rs,mid+1,r,x,mn);
	}
	inline void main()
	{
		for(int i = 1;i<=cnt;i++)
		{
			vec[l1[i]].push_back({{l2[i],r2[i]},s[i]}); 
			vec[r1[i]+1].push_back({{l2[i],r2[i]},-s[i]}); 
		}
		for(int i = 2;i<n;i++)
		{
			for(auto _:vec[i])
			{
				int l = _.first.first,r = _.first.second,v = _.second;
				change(1,2,m-1,l,r,v);
			}
			ask(1,2,m-1,i,2e9);
		}
		cout<<ans;
	}
}
signed main()
{
//	freopen(".in","r",stdin);
//	freopen(".out","w",stdout);
	ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
	cin>>n>>m;
	n+=2,m+=2;
	for(int i = 2;i<n;i++)
		for(int j = 2;j<m;j++)
			cin>>ch[i][j];
	for(int i = 1;i<=n;i++) ch[i][1] = ch[i][m] = '1';
	for(int i = 1;i<=m;i++) ch[1][i] = ch[n][i] = '1';
	for(int i = n;i;i--)
		for(int j = m;j;j--)
		{
			if(ch[i][j]=='1')
			{
				if(ch[i+1][j]=='1') dn[i][j] = dn[i+1][j];
				else dn[i][j] = i;
				if(ch[i][j+1]=='1') rt[i][j] = rt[i][j+1];
				else rt[i][j] = j;
			}
			else
			{
				if(ch[i+1][j]=='1') dn[i][j] = dn[i+1][j];
				if(ch[i][j+1]=='1') rt[i][j] = rt[i][j+1];
			}
		}
	long long sum = 0;
	for(int i = n;i;i--)
	{
		top[0] = 0;
		for(int j = m;j;j--)
		{
			int las1 = 0;
			while(top[0])
			{
				int now = stk[0][top[0]];
				while(top[j])
				{
					int k = stk[j][top[j]];
					if(rt[k][j]>rt[i][j]||rt[k][j]+1>=now) break;
					top[j]--;
				}
				int k = stk[j][top[j]];
				if(top[j]&&now<=min(rt[i][j],rt[k][j])+1&&k<=min(dn[i][j],dn[i][now])+1&&k>dn[i][las1])
				{
					if(i+1<k&&j+1<now)
					{
						int w = (k-i-1)*(now-j-1);
						sum+=w;
						cnt++;
						s[cnt] = w,l1[cnt] = i+1,r1[cnt] = k-1,l2[cnt] = j+1,r2[cnt] = now-1;
					}
				}
				if(dn[i][now]>dn[i][j]) break;
				las1 = stk[0][top[0]];
				top[0]--;
			}
			stk[0][++top[0]] = j;
			while(top[j]&&rt[stk[j][top[j]]][j]<=rt[i][j]) top[j]--;
			stk[j][++top[j]] = i;
		}
	}
	if(sum<=3e8) sub1::main();
	else sub2::main();
	return 0;
}