N=50;
I=ones(50);
I(1:5,6:10)=0;
I(25:45,30:45)=0;
U=I;

NA=0.85;lambda=193e-9;
pixel_size=5e-9;
x=linspace(-90,90,11);
y=linspace(-90,90,11);
[x,y]=meshgrid(x,y);
objective=sqrt(x.^2+y.^2)*pixel_size;
u=2*pi*objective*NA/lambda;
H=besselj(1,u)./u;
H(6,6)=0.5;

V=conv2(U,H,'same');
P=rand(50);
id = eye(N); 
S = circshift(id, [1,0]); 
D=id-S;
D(1,N)=0;
W_U=conv2(U,H,'same')+P./10;

x=binvar(50,50);
y=sdpvar(50,50,'full');
constraint=[conv2(double(x),H,'same')==double(y)];
objective=5*norm(y-W_U,'fro')^2+50*norm(D*x,'fro')^2+50*norm(x*D','fro')^2;
options.solver = 'gurobi';
options.verbose = 0;
result = optimize(constraint, objective, options);
if result.problem == 0
   U=value(x);  
else
   disp(result.info)
end
disp('U=')
disp(value(x))