Consider the following DP (Dynamic programming): \(dp[\) \(i\) \(][\) set of columns where all squares are painted \(] = \{\) number of ways to paint the first \(i\) lines \(\}\).
Specifically, each way of painting corresponds to the keys as follows:

        1101       0111       1111
        1111       0001       0011
        1001       1110       1101
dp  [3][1001]  [3][0000]  [3][0001]

Let’s consider the transition of this DP.
For example, dp[5][1101001] can transition to:

• dp[6][sssssss] by placing ( subset of 1101001) = sssssss;
• dp[6][110ssss] by placing 111 + ( subset of 1001 ) = 111ssss;
• dp[6][11010ss] by placing 11111 + ( subset of 01 ) = 11111ss; and
• dp[6][110100s] by placing 111111 + ( subset of 1 ) = 111111s.

As you can see, the transitions can be disjointly and fully classified by fixing the first occurrence of 0 in the key and letting the latter part be a subset of the key.
The transition to subsets can be done in the manner of Fast Zeta Transform.
(*) Note that the area of 1s expands from left, so one can apply the “transition to subsets” only once to finish the transitions by performing the transitions for \(k\) 1s in ascending order of \(k\), while, before every step, updating the DP table to represent the transitions of placing “(\(k\) 1s) + (subset of the others)”. (For more details, see the sample code.)

Now, the problem has been solved if the input grid is all ?. What if some of the squares are specified?
It is still not very difficult. As we described above, we know what to place in the row, so we can do the transitions only for those we can place.
For example, when applying a transition on a subset, if 1 cannot be placed, then the 1 in the key always transitions to 0; if 0 cannot be placed, the transition is impossible if the key is 0, and 1 is always placed if the key is 1.
The (*) paragraph can be achieved by placing \(k\) consecutive 1s from the left only when it is allowed. Again, for more details, see the sample code.

This way, the time complexity of \(O(2^M NM)\) is achieved, and you can solve this problem.

Sample code (C++):

#include<bits/stdc++.h>
#define mod 998244353

using namespace std;

int main(){
  int n,m;
  cin >> n >> m;
  vector<string> s(n);
  for(auto &nx : s){cin >> nx;}
  vector<long long> dp(1<<m,0);
  dp[(1<<m)-1]=1;
  for(int i=0;i<n;i++){
    vector<long long> ndp=dp;
    for(int j=0;j<m;j++){

      vector<long long> work(1<<m,0);
      if(s[i][j]!='1'){
        // put 0
        for(int k=0;k<(1<<m);k++){
          if(k&(1<<j)){
            work[k^(1<<j)]+=ndp[k];
            work[k^(1<<j)]%=mod;
          }
          else{
            work[k]+=ndp[k];
            work[k]%=mod;
          }
        }
      }
      if(s[i][j]!='0'){
        // put 1
        for(int k=0;k<(1<<m);k++){
          if(k&(1<<j)){
            work[k]+=ndp[k];
            work[k]%=mod;
          }
        }
      }
      ndp=work;

      bool ok=true;
      for(int k=0;k<=j;k++){
        if(s[i][k]=='0'){ok=false;break;}
      }
      if(!ok){continue;}
      for(int k=0;k<(1<<m);k++){
        if(k&(1<<j)){continue;}
        ndp[k]+=dp[k];
        ndp[k]%=mod;
      }
    }
    dp=ndp;
  }

  long long res=0;
  for(auto &nx : dp){res+=nx;res%=mod;}
  cout << res << "\n";
  return 0;
}