二维前缀DP求解矩阵区域和

矩阵区域和是一道经典的二维前缀和求解问题：

给定一个m*n的矩阵mat和一个整数k，对于矩阵中的每个位置(i,j)，需要计算以该位置为中心、边长为2k+1的正方形区域内所有元素的和，当边界超出矩阵时，会被截断。简单来说，就是给矩阵的每个位置画一个框，框的大小由k决定，然后把框里的所有数字加起来。

以示例1的中心位置(1,1)为例，k=1覆盖了周围3*3的所有元素，和为1+2+3+4+5+6+7+8+9=45。

求任意一个子矩阵的和，可以通过二维前缀和来进行求解，构建一个前缀和矩阵dp，其中dp[i][j]表示原矩阵从(0,0)到(i-1,j-1)这个子矩阵的所有元素之和。

int m=mat.size(),n=mat[0].size();
vector<vector<int>> dp(m+1,vector<int>(n+1,0));

vector<vector<int>> dp(m+1，vector<int>(n+1,0))，构建二维前缀和数组dp，多开一行一列方便处理边界问题，可通过容斥原理推出dp[i][j]的递推表达式。

dp[i][j]表示从原矩阵左上角(0,0)到(i-1,j-1)这个矩形内所有数的和：

可以把dp[i][j]看成这三部分组成：dp[i-1][j]、dp[i][j-1]、mat[i-1][j-1]，如下图所示：

1、上面区域：dp[i-1][j]，即前i-1列，前j列。

2、左边区域：dp[i][j-1]，即前i行，前j-1列。

右下角单独格子：mat[i-1][j-1]

如果直接将这三部分加起来，会发现左上角区域，即dp[i-1][j-1]被多加了一次，需要将其减掉，故最终dp[i][j]的递推表达式为：dp[i][j]=dp[i-1][j]+dp[i][j-1]-dp[i-1][j-1]+mat[i-1][j-1]。因此二维dp需多开一行一列，这样在计算第一行和第一列时，能安全地用到dp[0][j]或dp[i][0]这些为0的空区域。

class Solution {
public:
    vector<vector<int>> matrixBlockSum(vector<vector<int>>& mat, int k) {
      int m=mat.size(),n=mat[0].size();
      vector<vector<int>> dp(m+1,vector<int>(n+1));
      for(int i=1;i<=m;i++)
      {
          for(int j=1;j<=n;j++)
          {
              dp[i][j]=dp[i-1][j]+dp[i][j-1]-dp[i-1][j-1]+mat[i-1][j-1];
          }
      }
    }
};

有了递推公式dp[i][j]=dp[i-1][j]+dp[i][j-1]-dp[i-1][j-1]+mat[i-1][j-1]，通过两层for循环就可以对dp数组进行填充，使得dp[i][j]表示原矩阵从左上角到第i-1行第j-1列的所有元素之和。

有了前缀和数组dp之后，就可以通过dp来快速推出某个位置(i,j)的区域和，如下图所示：

首先分别确定以(i,j)为中心，半径为k的矩阵区域左上角、右下角的坐标，需要注意坐标是否越界，即左上角坐标为(max(i-k,0),max(j-k,0))，右下角坐标为(min(i+k,m-1),min(j+k,n-1))。

映射到dp，由于dp多增了一行一列，故左上角坐标(x1,y1)、右下角坐标(x2,y2)需相应+1。

问题就转化为：已知一个子矩阵，上边界为x1-1，下边界为x2-1，左边界为y1-1，右边界为y2-1，求这个子矩阵内所有数的和。

映射到dp，坐标相应加1，如下图所示：

要计算(x1,y1)到(x2,y2)的区域和，可以转化为先求(0,0)到(x2,y2)这块大矩形的区域和，该区域包含了四部分，如下图所示：

该区域包含了目标区域，但多了几块多余的部分：

1、上半部分：dp[x1-1][y2]

2、左半部分：dp[x2][y1-1]

根据容斥原理，用大矩形dp[x2][y2]减掉上半部分dp[x1-1][y2]、左半部分dp[x2][y1-1]，此时二者重叠部分dp[x1-1][y1-1]被多减了一次，要加回来，就求出了目标区域的和，即：ret=dp[x2][y2]-dp[x1-1][y2]-dp[x2][y1-1]+dp[x1-1][y1-1]。

代码实现：

class Solution {
public:
    vector<vector<int>> matrixBlockSum(vector<vector<int>>& mat, int k) {
      int m=mat.size(),n=mat[0].size();
      vector<vector<int>> dp(m+1,vector<int>(n+1));
      for(int i=1;i<=m;i++)
      {
          for(int j=1;j<=n;j++)
          {
              dp[i][j]=dp[i-1][j]+dp[i][j-1]-dp[i-1][j-1]+mat[i-1][j-1];
          }
      }
      vector<vector<int>> res(m,vector<int>(n));
      for(int i=0;i<m;i++)
      {
          for(int j=0;j<n;j++)
          {
              int x1=max(0,i-k)+1,x2=min(m-1,i+k)+1;
              int y1=max(0,j-k)+1,y2=min(n-1,j+k)+1;
              res[i][j]=dp[x2][y2]-dp[x1-1][y2]-dp[x2][y1-1]+dp[x1-1][y1-1];
          }
      }
      return res;  
    }
};

在原来代码的基础上，vector<vector<int>> res(m,vector<int>(n))，创建一个与原矩阵相同大小的矩阵res，通过for循环再次遍历原矩阵的每个位置(i,j)，int x1=max(0,i-k)+1,x2=min(m-1,i+k)+1;int y1=max(0,j-k)+1,y2=min(n-1,j+k)+1;根据参数k计算出以该位置(i,j)为中心的矩阵边界，+1转化为在dp数组中的索引，ret[i][j]=dp[x2][y2]-dp[x1-1][y2]-dp[x2][y1-1]+dp[x1-1][y1-1]，最后通过容斥原理利用四个前缀和O(1)求得该区域的和，并存入res[i][j]中，return res即可。整个算法的时间复杂度为O(m*n)，空间复杂度也为O(m*n)。

矩阵区域和的本质是通过二维前缀和来求解，而二维前缀和的核心就是容斥原理，分为两步，构建前缀和数组dp时，dp[i][j]表示从原矩阵左上角(0,0)到(i-1,j-1)的矩形和，要得到dp[i][j]，可以通过容斥原理把上方区域dp[i-1][j]、左方区域dp[i][j-1]和当前格mat[i-1][j-1]加起来，但此时左上角dp[i-1][j-1]被重复加了一次，因此要减掉，即dp[i][j]=dp[i-1][j]+dp[i][j-1]-dp[i-1][j-1]+mat[i-1][j-1]，多开一行一列，是为了让边界情况也能使用该公式。查询任意子矩形的区域和，也是通过容斥原理，目标区域=大矩形dp[x2][y2]减去上方部分dp[x1-1][y2]减去左方部分dp[x2][y1-1]，但此时左上角dp[x1-1][y1-1]被多减了一次，要加回来，故res=dp[x2][y2]-dp[x1-1][y2]-dp[x2][y1-1]+dp[x1-1][y1-1]，两个公式本质都是容斥，构建时从小到大累加，加两块减一块，查询时从大到小裁剪，减两块加一块，一正一反，互为逆运算。这样就能在O(m*n)预处理后，O(1)求解任意矩阵区域和。