%% Predictive Cultural Heritage Digital Twin – Enhanced Visualization

clear; clc; close all;

%% 1. Baseline Artifact
[X,Y,Z] = sphere(50);
X = X(:) + 0.01*randn(numel(X),1);
Y = Y(:) + 0.01*randn(numel(Y),1);
Z = Z(:) + 0.01*randn(numel(Z),1);
PC_baseline = pointCloud([X,Y,Z]);

figure('Name','Baseline Artifact');
pcshow(PC_baseline);
title('Baseline Artifact Point Cloud');
xlabel('X'); ylabel('Y'); zlabel('Z');
axis equal; view(30,20);

%% 2. Environmental Conditions
dt = 60; Tf = 6*3600; N = Tf/dt;
Tout = 18 + 2*sin((1:N)*2*pi/1440);   % temperature °C
RHout = 45 + 5*sin((1:N)*2*pi/1440); % relative humidity %

% Plot environmental conditions
figure('Name','Environmental Conditions');
subplot(2,1,1);
plot(1:N, Tout, 'r'); hold on; plot(1:N, RHout, 'b');
legend('Outdoor Temp [°C]','Outdoor RH [%]');
xlabel('Time steps'); ylabel('Environmental Values');
title('Simulated Outdoor Conditions');

subplot(2,1,2);
plot(1:N, Tout, 'r'); hold on; plot(1:N, RHout, 'b');
xlabel('Time steps'); ylabel('Values');
xlim([1 600]); % Zoom in first 10 minutes for clarity
title('Environmental Conditions (Zoomed)');

%% 3. Artifact Microclimate
C_T = 1e4; R_T = 20;  
C_RH = 50; R_RH = 5;  
T_art = zeros(1,N); RH_art = zeros(1,N);  
T_art(1) = 20; RH_art(1) = 50;

for k = 2:N
    T_art(k) = T_art(k-1) + dt/C_T * ((Tout(k-1)-T_art(k-1))/R_T);
    RH_art(k) = RH_art(k-1) + dt/C_RH * ((RHout(k-1)-RH_art(k-1))/R_RH);
end

% Plot microclimate
figure('Name','Artifact Microclimate');
subplot(2,1,1); plot(1:N,T_art,'r'); hold on; plot(1:N,T_art,'ro','MarkerSize',2);
xlabel('Time steps'); ylabel('Artifact Temp [°C]');
title('Artifact Temperature Over Time');

subplot(2,1,2); plot(1:N,RH_art,'b'); hold on; plot(1:N,RH_art,'bo','MarkerSize',2);
xlabel('Time steps'); ylabel('Artifact RH [%]');
title('Artifact Relative Humidity Over Time');

%% 4. Add Measurement Noise
sigma_T = 0.3; sigma_RH = 1.0;
T_meas = T_art + sigma_T*randn(size(T_art));
RH_meas = RH_art + sigma_RH*randn(size(RH_art));

% Plot noisy measurements
figure('Name','Noisy Measurements');
subplot(2,1,1); plot(T_art,'r'); hold on; plot(T_meas,'k.');
legend('True Temp','Measured Temp'); title('Temperature Measurements');
xlabel('Time steps'); ylabel('Temp [°C]');

subplot(2,1,2); plot(RH_art,'b'); hold on; plot(RH_meas,'k.');
legend('True RH','Measured RH'); title('Humidity Measurements');
xlabel('Time steps'); ylabel('RH [%]');

%% 5. Kalman Filter Estimation
a_T = 1 - dt/(R_T*C_T); b1_T = dt/(R_T*C_T);
a_RH = 1 - dt/(R_RH*C_RH); b1_RH = dt/(R_RH*C_RH);

Qw_T = 1e-3; Qw_RH = 1e-2;
Rv_T = sigma_T^2; Rv_RH = sigma_RH^2;

xhat_T = zeros(1,N); P_T = 1; xhat_T(1) = T_meas(1);
xhat_RH = zeros(1,N); P_RH = 1; xhat_RH(1) = RH_meas(1);

for k = 2:N
    % Temperature
    x_pred = a_T*xhat_T(k-1) + b1_T*Tout(k-1);
    P_pred = a_T*P_T*a_T + Qw_T;
    K = P_pred/(P_pred + Rv_T);
    xhat_T(k) = x_pred + K*(T_meas(k)-x_pred);
    P_T = (1-K)*P_pred;
    
    % Humidity
    x_pred = a_RH*xhat_RH(k-1) + b1_RH*RHout(k-1);
    P_pred = a_RH*P_RH*a_RH + Qw_RH;
    K = P_pred/(P_pred + Rv_RH);
    xhat_RH(k) = x_pred + K*(RH_meas(k)-x_pred);
    P_RH = (1-K)*P_pred;
end

% Plot estimated microclimate
figure('Name','Kalman Filter Estimates');
subplot(2,1,1);
plot(T_art,'r'); hold on; plot(T_meas,'k.'); plot(xhat_T,'b-');
legend('True','Measured','Estimated'); title('Temperature Estimation');

subplot(2,1,2);
plot(RH_art,'b'); hold on; plot(RH_meas,'k.'); plot(xhat_RH,'r-');
legend('True','Measured','Estimated'); title('Humidity Estimation');

%% 6. Time-Lapse Artifact Scans with Environmental Degradation
num_scans = 5;
degradation_factor = 0.001;
PC_scans = cell(1,num_scans);
PC_scans{1} = PC_baseline;

figure('Name','Time-Lapse Scans');
for k = 2:num_scans
    theta = deg2rad(5*randn);
    Ry = [cos(theta) 0 sin(theta);0 1 0;-sin(theta) 0 cos(theta)];
    t = 0.005*randn(1,3);
    
    pts = (Ry*PC_scans{k-1}.Location')' + t;
    env_factor = 1 + 0.5*max(0,xhat_T(k-1)-22)/5 + 0.5*max(0,xhat_RH(k-1)-55)/10;
    pts = pts + degradation_factor*env_factor*randn(size(pts));
    PC_scans{k} = pointCloud(pts);
    
    subplot(2,3,k-1);
    pcshow(PC_scans{k}); title(['Scan ',num2str(k)]);
    axis equal; view(30,20);
end

%% 7. ICP Registration and Merge
PC_registered = PC_scans{1};
for k = 2:num_scans
    [~, PC_reg, ~] = pcregistericp(PC_scans{k}, PC_registered, ...
        'Metric','pointToPlane','MaxIterations',60);
    PC_registered = pcmerge(PC_registered, PC_reg, 0.005);
end

figure('Name','Merged Artifact');
pcshow(PC_registered);
title('Merged Artifact After ICP Registration'); axis equal;

%% 8. Compute and Visualize Damage Heatmap
baseline_pts = PC_baseline.Location;
merged_pts = PC_registered.Location;
[idx,dists] = knnsearch(baseline_pts, merged_pts);
damage_map = zeros(size(baseline_pts,1),1);
for i = 1:size(baseline_pts,1)
    matches = find(idx==i);
    if ~isempty(matches)
        damage_map(i) = mean(dists(matches));
    end
end

figure('Name','Damage Heatmap');
pcshow(PC_baseline.Location, damage_map);
title('Artifact Surface Damage Heatmap');
colorbar; axis equal;

%% 9. Predictive Alerts
threshold = 0.003;
num_alerts = sum(damage_map>threshold);
fprintf('Points exceeding threshold: %d\n', num_alerts);
if num_alerts>0
    disp('ALERT: Areas of artifact are at risk!');
else
    disp('Artifact condition within safe limits.');
end
