Python3, 562 bytes

Long, but fast.

E=enumerate
def f(b):
 q,s,F,W=[(b,0,0)],[],[],-1
 for b,w,c in q:
  o={j for x,r in E(b)for y,v in E(r)for X,Y in[(1,0),(-1,0),(0,1),(0,-1),(1,1),(-1,1),(1,-1),(-1,-1)]for j in[(0,min(x,x+X)),(1,min(y,y+Y))]if v and 0<=x+X<len(b)and 0<=y+Y<len(b[0])and b[x+X][y+Y]}
  if not o:F+=[(b,w,c)];W=w
  for d,x in o:
   l=[b,[*zip(*b)]][d];l=l[:x+1]+[L:=[0]*len(l[x])]+l[x+1:];l=[l,[*map(list,zip(*l))]][d]
   if(W==-1 or w+len(L)<=W)and l not in s:q+=[(l,w+len(L),c+d)];s+=[l]
 K={}
 for b,w,c in F:K[w]=K.get(w,[])+[(b,c)]
 return max(K[min(K)],key=lambda x:x[1])[0]

Try it online!

Python 2, 374 346 340 339 323 317 bytes

R=range;L=len
def f(a):
 w,h=L(a[0]),L(a);W=[]
 for i in R(2**w):
	A=zip(*a)
	for c in R(w):A[-c:-c]=[[0]*h]*(i&1<<c>0)
	for j in R(2**h):
	 B=zip(*A);x=L(B[0])
	 for r in R(h):B[-r:-r]=[(0,)*x]*(j&1<<r>0)
	 y=L(B);W+=[(w*h-x*y,x,B)]*all(sum(B[i][j:j+2]+B[i+1][j:j+2])<2for i in R(y-1)for j in R(x))
 return max(W)[2]

Try it online!

Jelly, 37 bytes

ṫƤ-S€ZƊ⁺FỊẠ
Z_,,WƲ€ŒpẎ€Ʋ⁺€ẎLÞFL$ÞṚÇÞṪ

Try it online!

Function returning a 2D array of integers. Note that naturally in Jelly singleton list displays as its value so G is used to format the output.

• Link 1: Return (validity).
• Link 2: Main program.

The program runs in exponential time, but so-far I could not think of any polynomial time algorithm. Uses Ƥ on dyadic function, that feature postdates the challenge.