function [cams,P,rayf,uv_hat] = im_simulation(X_hat,X,dP,cams,camera,lato,pc,theta)
%INPUT:     X_hat=markers positions prediction
%           dP=variance of markers positions prediction
%           cams=struct
%           camera=which camera
%           lato=image plane dimension in pixel
%           pc=percentual of occluded markers
%OUTPUT:    Fictious=fictious image
%           Real=real image
%           X=real markers position
%


x_hat=X_hat(:,cams(camera).setmarkers(1,:))
x=X(:,cams(camera).setmarkers(1,:))
dp=dP(cams(camera).setmarkers(1,:))

[m,n]=size(x_hat);
K=cams(camera).intrinsic;
R=cams(camera).rotationmat;
T=cams(camera).position;
f=(lato/2)*(cos(theta/2)/sin(theta/2));
[uv_hat, X_cam, hidden, distance]=camera_simulation(R,T,x_hat,K,f,lato,lato);
X_cam=X_cam(1:3,:);
% Fictious=zeros(lato,lato);
% MPP=K*[R' -R'*T];
rayf=zeros(1,n);
for i=1:n
%     distanza=(X_cam(1:3,i)'*X_cam(1:3,i))^0.5;
%     d_max_var=5;
%     d_min_var=0.5;
%     rays(i)=ray_distanza(distanza,d_min_var,d_max_var,dp(i));
%     X_hat(:,i)
%     P=[X_cam(1,i)+dp(i) X_cam(2,i) X_cam(3,i)]';
%     P=[R T]*[P;1]
%     uvs=MPP*[P;1];
%     uv_hat
%     uvn=[uvs(1)/uvs(3); uvs(2)/uvs(3)]
%     uv=uvn-uv_hat(:,i);
%     rays(i)=norm(uv);
    d=sqrt(X_cam(:,i)'*X_cam(:,i));
    rayf(i)=(3*dp(i)^0.5)*(f/d);
end

% P=[uv_hat' distance'];
% Fictiouscntrd=genera_piano_3D(P, lato, lato, 0.5, 5, 50, 10000);


% for i=1:lato
%     for j=1:lato
%         for k=1:n
%             d=((uv_hat(:,k)-[i j]')'*(uv_hat(:,k)-[i j]'))^0.5;
%             if d<rayf(k)
%                 Fictious(j,i)=130;
%             end
%         end
%     end
% end

[uv_real, x_cam, hidden, distance1]=camera_simulation(R,T,x,K,f,lato,lato);

uv=[];
distancer=[];
set=[];
for i=1:n
    a=rand(1);
    if a>pc
        set=[set i];
    end
end

X_view=cams(camera).setmarkers(1,set);
cams(camera).setmarkers_visti=X_view;
uv=uv_real(:,set);
rayr=zeros(1,size(X_view,2));
distancer=distance1(set);

for i=1:size(X_view,2)
    rayr(i)=round((3*dp(i)^0.5)*(f/distance(i)));
end

Real=zeros(lato,lato);
for k=1:size(distancer,2)
    min_i=uv(1,k)-rayr(k);
    if min_i<1
        min_i=1;
    end
    max_i=uv(1,k)+rayr(k);
    if max_i>lato
        max_i=lato;
    end
    
    min_j=uv(1,k)-rayr(k);
    if min_j<1
        min_j=1;
    end
    max_j=uv(1,k)+rayr(k);
    if max_j>lato
        max_j=lato;
    end
    
    for i=round(min_i):round(max_i)
        for j=round(min_j):round(max_j)
            d=((uv_hat(:,k)-[i j]')'*(uv_hat(:,k)-[i j]'))^0.5;
            if d<rayr(k)
                Real(j,i)=130;
            end
        end
    end
end

    
P=[uv' distancer'];


end

