zhaodongsun / contrast-phys Goto Github PK

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch, resample
from utils_sig import hr_fft, butter_bandpass, normalize

# 读取 .npy 文件
data = np.load(r'C:\first-year graduate\contrast-phys-master\results\2\3\subject38_pred.npy', allow_pickle=True)
rppg_list = data.item().get('rppg_list')
bvp_list = data.item().get('bvp_list')

# 展平信号
rppg = np.concatenate(rppg_list)
bvp = np.concatenate(bvp_list)

# 标准化信号
rppg_norm = normalize(rppg)
bvp_norm = normalize(bvp)

# 计算信号长度和采样率
signal_length = len(rppg_norm)
fs = 30  # 假设采样率为 30Hz

# 对标准化信号进行滤波
rppg_filtered = butter_bandpass(rppg_norm, lowcut=0.6, highcut=4, fs=fs)
bvp_filtered = butter_bandpass(bvp_norm, lowcut=0.6, highcut=4, fs=fs)

# 计算预测心率和地面真值心率
hr_rppg, rppg_psd, rppg_hr = hr_fft(rppg_filtered, fs=fs)
hr_bvp, bvp_psd, bvp_hr = hr_fft(bvp_filtered, fs=fs)

# 计算评估指标
mae = np.mean(np.abs(rppg_filtered - bvp_filtered))
mse = np.mean((rppg_filtered - bvp_filtered) ** 2)
rmse = np.sqrt(mse)
r = np.corrcoef(rppg_filtered, bvp_filtered)[0, 1]

# 绘制时域图和频谱图
fig, (ax1, ax2) = plt.subplots(2, figsize=(20, 10))

# 时域图
time = np.arange(signal_length) / fs
ax1.plot(time, rppg_filtered, label='Predicted Signal')
ax1.plot(time, bvp_filtered, label='Ground Truth Signal')
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Amplitude')
ax1.set_title('Predicted Signal vs Ground Truth Signal')
ax1.legend()
ax1.grid('on')

# 频谱图
ax2.plot(rppg_hr, rppg_psd, label='Predicted Signal')
ax2.plot(bvp_hr, bvp_psd, label='Ground Truth Signal')
ax2.set_xlabel('Heart Rate (bpm)')
ax2.set_ylabel('Power Spectral Density')
ax2.set_xlim([40, 200])
ax2.set_title('Power Spectral Density of Predicted Signal vs Ground Truth Signal')
ax2.legend()
ax2.grid('on')

# 显示评估指标和预测心率
print(f"MAE: {mae:.4f}")
print(f"MSE: {mse:.4f}")
print(f"RMSE: {rmse:.4f}")
print(f"Pearson's Correlation Coefficient (r): {r:.4f}")
print(f"Predicted Heart Rate: {hr_rppg:.2f} bpm")
print(f"Ground Truth Heart Rate: {hr_bvp:.2f} bpm")

plt.show()

for h5_path in test_list:
        h5_path = str(h5_path)

        with h5py.File(h5_path, 'r') as f:
            imgs = f['imgs']
            subject_name = os.path.basename(h5_path)[:-3]
            bvp_path = f"/share2/data/zhouwenqing/UBFC_rPPG/dataset2/{subject_name}/ground_truth.txt"
            bvp = np.loadtxt(bvp_path).reshape((-1, 1))
            # bvppeak = f['bvp_peak']
            fs = config_train['fs']
            
            # duration表示秒数，fs表示frame per seccond 
            duration = np.min([imgs.shape[0], bvp.shape[0]]) / fs
            num_blocks = int(duration // time_interval)
            # 从整个视频当中截取出num_blocks个视频片段，这些片段之间是连续的（指从原视频当中截取的方式）
            rppg_list = []
            bvp_list = []
            # bvppeak_list = []
            for b in range(num_blocks):
                rppg_clip = dl_model(imgs[b*time_interval*fs:(b+1)*time_interval*fs])
                rppg_list.append(rppg_clip)

                bvp_list.append(bvp[b*time_interval*fs:(b+1)*time_interval*fs])
                # bvppeak_list.append(bvppeak[b*time_interval*fs:(b+1)*time_interval*fs])

            rppg_list = np.array(rppg_list)
            bvp_list = np.array(bvp_list)
            results = {'rppg_list': rppg_list, 'bvp_list': bvp_list}
            np.save(pred_exp_dir+'/'+h5_path.split('/')[-1][:-3], results)
            
            bvp_list = bvp_list.reshape(num_blocks, -1)
  
            
            # 通带滤波
            hr_pred = butter_bandpass_batch(rppg_list, lowcut=0.6, highcut=4, fs=30)
            hr_gt = butter_bandpass_batch(bvp_list, lowcut=0.6, highcut=4, fs=30)
            print("通带滤波 hr_pred.shape: ", hr_pred.shape)
            print("通带滤波 hr_gt.shape: ", hr_gt.shape)
            
            # 计算预测心率和地面真值心率
            hr_pred = torch.tensor(hr_fft_batch(hr_pred, 30))
            hr_gt = torch.tensor(hr_fft_batch(hr_gt, 30))
            print("hr_fft_batch hr_pred.shape: ", hr_pred.shape)
            print("hr_fft_batch hr_gt.shape: ", hr_gt.shape)
            
            mae_all = mae_loss_func(hr_pred, hr_gt)
            mse_all = mse_loss_func(hr_pred, hr_gt)
            rmse = np.sqrt(mse_all)
            correlation_coefficients = []
            for i in range(num_blocks):
                # 计算每个维度的皮尔逊相关系数
                correlation_coefficient, _ = pearsonr(hr_pred[i].flatten(), hr_gt[i].flatten())
                correlation_coefficients.append(correlation_coefficient)
            print("Evaluation Result\n MAE: {:.4f}; RMSE: {:.4f}; R: {:.4f};".format(
                mae_all, rmse, sum(correlation_coefficients) / len(correlation_coefficients)))

File "test.py", line 132, in my_main
correlation_coefficient, _ = pearsonr(hr_pred[i].flatten(), hr_gt[i].flatten())
The File "/ share1 / home/zhouwenqing/anaconda3 / envs/RPPG/lib/python3.6 / site - packages/scipy/stats/stats. Py", line 3838, in pearsonr
raise ValueError('x and y must have length at least 2.')
ValueError: x and y must have length at least 2.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch, resample
from utils_sig import hr_fft, butter_bandpass, normalize
import os


root_dir = "./results_cnn_pure/1/1"
for file in os.listdir(root_dir):
    
    if not file.endswith(".npy"):
        continue
    print(file)
    data = np.load(os.path.join(root_dir, file), allow_pickle=True)
    rppg_list = data.item().get('rppg_list')
    bvp_list = data.item().get('bvp_list')
    
    # 展平信号
    rppg = np.concatenate(rppg_list)
    # print(rppg.shape) # (1500,)
    bvp = np.concatenate(bvp_list)

    # 标准化信号
    rppg_norm = normalize(rppg)
    bvp_norm = normalize(bvp)
    

    # 计算信号长度和采样率
    signal_length = len(rppg_norm)
    fs = 30  # 假设采样率为 30Hz

    # 对标准化信号进行滤波
    rppg_filtered = butter_bandpass(rppg_norm, lowcut=0.6, highcut=4, fs=fs)
    bvp_filtered = butter_bandpass(bvp_norm, lowcut=0.6, highcut=4, fs=fs)

    # 计算预测心率和地面真值心率
    hr_rppg, rppg_psd, rppg_hr = hr_fft(rppg_filtered, fs=fs)
    hr_bvp, bvp_psd, bvp_hr = hr_fft(bvp_filtered, fs=fs)

    # 计算评估指标
    mae = np.mean(np.abs(rppg_filtered - bvp_filtered))
    mse = np.mean((rppg_filtered - bvp_filtered) ** 2)
    rmse = np.sqrt(mse)
    # print(np.corrcoef(rppg_filtered, bvp_filtered))
    r = np.corrcoef(rppg_filtered, bvp_filtered)[0, 1]

    # 绘制时域图和频谱图
    fig, (ax1, ax2) = plt.subplots(2, figsize=(20, 10))

    # 时域图
    time = np.arange(signal_length) / fs
    ax1.plot(time, rppg_filtered, label='Predicted Signal')
    ax1.plot(time, bvp_filtered, label='Ground Truth Signal')
    ax1.set_xlabel('Time (s)')
    ax1.set_ylabel('Amplitude')
    ax1.set_title('Predicted Signal vs Ground Truth Signal')
    ax1.legend()
    ax1.grid('on')

    # 频谱图
    ax2.plot(rppg_hr, rppg_psd, label='Predicted Signal')
    ax2.plot(bvp_hr, bvp_psd, label='Ground Truth Signal')
    ax2.set_xlabel('Heart Rate (bpm)')
    ax2.set_ylabel('Power Spectral Density')
    ax2.set_xlim([40, 200])
    ax2.set_title('Power Spectral Density of Predicted Signal vs Ground Truth Signal')
    ax2.legend()
    ax2.grid('on')

    # 显示评估指标和预测心率
    print(f"MAE: {mae:.4f}")
    print(f"MSE: {mse:.4f}")
    print(f"RMSE: {rmse:.4f}")
    print(f"Pearson's Correlation Coefficient (r): {r:.4f}")
    print(f"Predicted Heart Rate: {hr_rppg:.2f} bpm")
    print(f"Ground Truth Heart Rate: {hr_bvp:.2f} bpm")

    plt.show()

import numpy as np
import h5py
from torch import nn
import torch
from PhysNetModel import PhysNet
from utils_data import *
from utils_sig import *
from sacred import Experiment
from sacred.observers import FileStorageObserver
import json
from scipy.stats import pearsonr
from ViT import VideoViT
from ViT2 import TransformerImageResizer
from UniFormer import Uniformer
from CNNandTransformer import MyNet
from scipy import interpolate
import json

ex = Experiment('model_pred', save_git_info=False)

@ex.config
def my_config():
    e = 29 # the model checkpoint at epoch e
    train_exp_num = 3 # the training experiment number
    train_exp_dir = './results4cnn_pure/%d'%train_exp_num # training experiment directory
    
    time_interval = 10 # get rppg for 30s video clips, too long clips might cause out of memory

    ex.observers.append(FileStorageObserver(train_exp_dir))

    if torch.cuda.is_available():
        device = torch.device('cuda')
        torch.backends.cudnn.enabled = True
        torch.backends.cudnn.benchmark = True
    
    else:
        device = torch.device('cpu')

@ex.automain
def my_main(_run, e, train_exp_dir, device, time_interval):

    # load test file paths
    test_list = list(np.load(train_exp_dir + '/test_list.npy'))
    pred_exp_dir = train_exp_dir + '/%d'%(int(_run._id)) # prediction experiment directory

    with open(train_exp_dir+'/config.json') as f:
        config_train = json.load(f)

    model = PhysNet(config_train['S'], config_train['in_ch']).to(device).eval()
    
    
    model.load_state_dict(torch.load(train_exp_dir+'/epoch%d.pt'%(e), map_location=device)) # load weights to the model
    # model.load_state_dict(torch.load("./demo/model_weights.pt", map_location=device))

    @torch.no_grad()
    def dl_model(imgs_clip):
        # model inference
        img_batch = imgs_clip
        img_batch = img_batch.transpose((3,0,1,2))
        # 在img_batch前面新增了一个批量大小的维度（批量大小为1）
        img_batch = img_batch[np.newaxis].astype('float32')
        img_batch = torch.tensor(img_batch).to(device)

        rppg = model(img_batch)[:,-1, :] # (1, 5, T) -> (1, T)
        rppg = rppg[0].detach().cpu().numpy()
        return rppg
    
    @torch.no_grad()
    def read_wave(info, target_length):
        bvp = []
        hr = []
        for signal in info["/FullPackage"]:
            bvp.append(signal["Value"]["waveform"])
            hr.append(signal["Value"]["pulseRate"])
        
        T = len(bvp)
        gts = list()
        x_new = np.linspace(0, T - 1, num=target_length)  
        bvp_func = interpolate.CubicSpline(range(T), bvp)
        gts.append(bvp_func(x_new))
        hr_func = interpolate.CubicSpline(range(T), hr)
        gts.append(hr_func(x_new))

        return np.asarray(gts)

    
    for h5_path in test_list:
        h5_path = str(h5_path)
        # print(h5_path)
        with h5py.File(h5_path, 'r') as f:
            imgs = f['imgs']
            subject_name = os.path.basename(h5_path)[:-3]
            json_path = f"/share2/data/PURE/{subject_name}/{subject_name}.json"
            with open(json_path) as f:
                info = json.load(f)
                gts = read_wave(info, imgs.shape[0])
                # print(gts.shape)
            bvp = gts[0].reshape(-1, 1)
            # bvppeak = f['bvp_peak']
            fs = config_train['fs']
            
            print(subject_name, imgs.shape, bvp.shape)
            
            # duration表示秒数，fs表示frame per seccond 
            duration = np.min([imgs.shape[0], bvp.shape[0]]) / fs
            num_blocks = int(duration // time_interval)
            # 从整个视频当中截取出num_blocks个视频片段，这些片段之间是连续的（指从原视频当中截取的方式）
            rppg_list = []
            bvp_list = []
            # bvppeak_list = []
            for b in range(num_blocks):
                rppg_clip = dl_model(imgs[b*time_interval*fs:(b+1)*time_interval*fs])
                rppg_list.append(rppg_clip)

                bvp_list.append(bvp[b*time_interval*fs:(b+1)*time_interval*fs])
                # bvppeak_list.append(bvppeak[b*time_interval*fs:(b+1)*time_interval*fs])

            rppg_list = np.array(rppg_list)
            bvp_list = np.array(bvp_list)
            bvp_list = bvp_list.reshape(num_blocks, -1)
            print(num_blocks, rppg_list.shape, bvp_list.shape)
            # bvppeak_list = np.array(bvppeak_list)
            # results = {'rppg_list': rppg_list, 'bvp_list': bvp_list, 'bvppeak_list':bvppeak_list}
            results = {'rppg_list': rppg_list, 'bvp_list': bvp_list}
            np.save(pred_exp_dir+'/'+h5_path.split('/')[-1][:-3], results)