diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..4dac285 --- /dev/null +++ b/.gitignore @@ -0,0 +1,15 @@ +code/__pycache__/* + +rPPG-checkpoints +rPPG-checkpoints/* +.vscode +.vscode/* +.vs +.vs/* +0.png +1.png +log.txt +.gitignore +picture1.jpg +picture.jpg +picture2.jpg diff --git a/README.md b/README.md index c46b40d..1df33b0 100644 --- a/README.md +++ b/README.md @@ -1,86 +1,110 @@ -## MTTS-CAN: Multi-Task Temporal Shift Attention Networks for On-Device Contactless Vitals Measurement +# Analysis and optimization of photoplethysmography imaging methods for non-contact measurement of heart variability parameters [![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT) [![made-with-python](https://img.shields.io/badge/Made%20with-Python-1f425f.svg)](https://www.python.org/) +Deep learning neuronal networks based on remote photoplethysmography. Extracting the pulse signal from video using machine learning with a view to heart rate variability parameters. +Source code of the master thesis titles: "Analysis and optimization of photoplethysmography imaging methods for non-contact measurement of heart variability parameters" + +## Cite as + +Sarah Quehl. (2022, April 11). Analysis and Optimization of Photoplethysmography Imaging Methods for Non-Contact Measurement of Heart Variability Parameters + +## Abstract +Heart rate variability is an important physiological parameter for health and refers to the natural +variation of the time between two heartbeats. Heart rate variability describes the adaptability of an +organism to external and internal factors and can be measured with common measuring devices, like +electrocardiogram or photoplethysmogram. Today, this is even possible with the smartphone via apps. +Photoplethysmography Imaging as a non-contact method is a further development of state-of-the-art +photoplethysmography for recording cardiac activity by detecting minimal pulse-induced fluctuations +on the skin with a RGB camera. Most Photoplethysmography Imaging methods focus on heart rate +measurement and do not consider heart rate variability. In recent years, many new approaches based +on signal filtering or neural networks have been presented. However, the accuracy required for medical +purposes, especially with regard to heart rate variability, has not yet been achieved and represents a +major challenge. +This thesis compares current Photoplethysmography Imaging methods based on neural networks. For +this purpose, four basis methods are implemented and tested for functionality. Based on these findings, +two new networks were developed, the PTS-CAN and the PPTS-CAN. These are based on multi-objective +optimization and add one and two additional outputs to the neural network, respectively. The additional +output of the PTS-CAN outputs a binary signal that has a value of one at peaks. For this output two new +loss functions were developed, which have the goal to reduce the temporal error of the peaks. For this +purpose, two new loss functions named ownGauss and the TE were developed, the last one allows an +interpretation of the error in seconds. Both manipulate the ground truth to generate a loss, to reward +the peaks that are close to the real peak and to punish peaks that are further away. A further output +was added to the first model, which outputs various variable parameters and is evaluated by the mean +absolute percentage error loss function. All used models are trained on the same database and are +compared. In addition, there is a comparison with the first developed methods on the subject of vital +parameters extraction from video. A final comparison shows an improvement in HR and HRV parameter +calculation with the new methods. The heart rate calculation can be improved by about 20%. In the field +of HRV parameters, an improvement of 5,7% can be achieved for the parameter SDNN, for example. +In a cross-validation, improvements are achieved over the baseline methods and there is also a slight +improvement over the basis models. For the parameters in the frequency domain, the improvements are +a bit less clear than in the time domain, since the frequency analysis is more challenging here. +A project was generated, which can be used as a basis for further experiments with further approaches +and loss functions. The integration of further network architectures as well as loss functions is easily +possible. + +## Preprocessing +It is recommended to save the important information of each video into a hdf5-file using the `prepare_databases.py` script. Here pixel data, ground truth and various parameters are integrated. + +## Training -## Paper - -#### [Xin Liu](https://homes.cs.washington.edu/~xliu0/), [Josh Fromm](https://www.linkedin.com/in/josh-fromm-2a4a2258/), [Shwetak Patel](https://ubicomplab.cs.washington.edu/members/), [Daniel McDuff](https://www.microsoft.com/en-us/research/people/damcduff/), “Multi-Task Temporal Shift Attention Networks for On-Device Contactless Vitals Measurement”, NeurIPS 2020, Oral Presentation (105 out of 9454 submissions) - -#### Link: - - -## New Pre-Trained Model (Updated Nov 2021) - -Working in pregress. Check back later! - -#### Abstract - -Telehealth and remote health monitoring have become increasingly important during the SARS-CoV-2 pandemic and it is widely expected that this will have a lasting impact on healthcare practices. These tools can help reduce the risk of exposing patients and medical staff to infection, make healthcare services more accessible, and allow providers to see more patients. However, objective measurement of vital signs is challenging without direct contact with a patient. We present a video-based and on-device optical cardiopulmonary vital sign measurement approach. It leverages a novel multi-task temporal shift convolutional attention network (MTTS-CAN) and enables real-time cardiovascular and respiratory measurements on mobile platforms. We evaluate our system on an ARM CPU and achieve state-of-the-art accuracy while running at over 150 frames per second which enables real-time applications. Systematic experimentation on large benchmark datasets reveals that our approach leads to substantial (20\%-50\%) reductions in error and generalizes well across datasets. - - +`python code/train.py --exp_name test --exp_name [e.g., test] --data_dir [DATASET_PATH] --temporal [e.g., MMTS_CAN]` -## Waveform Samples +examples: -### Pulse +python code/train.py --exp_name test1 --data_dir /mnt/share/StudiShare/sarah/Databases/ --temporal TS_CAN --database_name MIX2 -![pulse_waveform](./pulse_waveform.png) +#### Issues: -### Respiration +In PPTS_CAN, the frame rate used is derived from the video length used (which results from the data sets). This must still be passed in generalized form in the layers. -![resp_waveform](./resp_waveform.png) +## Inference +`python code/predict_vitals_oneVideo.py --video_path [VIDEO_PATH] --save_dir [SAVE_PATH] --trained_model [CHECKPOINT_PATH] + --model_name [e.g., TS_CAN, PTS_CAN, PPTS_CAN] --parameter [e.g., "bpm, sdnn, pnn50, lfhf"]` -## Citation +## Path dependencies in the following scripts +final_evaluation.py -``` bash -@article{liu2020multi, - title={Multi-Task Temporal Shift Attention Networks for On-Device Contactless Vitals Measurement}, - author={Liu, Xin and Fromm, Josh and Patel, Shwetak and McDuff, Daniel}, - journal={arXiv preprint arXiv:2006.03790}, - year={2020} -} -``` +model_evaluation.py -## Demo +pre_process.py -**Try out our live demo via link [here](https://vitals.cs.washington.edu/).** +predict_vitals_comparison.py -Our demo code: https://github.com/ubicomplab/rppg-web +predict_vitals_new.py +predict_vitals_oneVideo.py -## TVM +predict.vitals.py -If you want to use TVM, pleaea follow [this tutorial](https://tvm.apache.org/docs/) to set it up. Then, you will need to replace the code in `incubator-tvm/python/tvm/relay/frontend/keras.py` with our `code/tvm-ops-mtts-can.py`. We implemented required tensor operations for attention, tensor shift module used in our models. +layer_output.py -## Training -`python code/train.py --exp_name test --exp_name [e.g., test] --data_dir [DATASET_PATH] --temporal [e.g., MMTS_CAN]` +In the current scripts, the data has been divided into the folders 1)Training and 2)Validation. -## Inference +## evaluation_iPhys.py +Script for evaluating the prediction of the iPhys models (GreenChannel, POH, CHROM) with the same procedure and products as in the finalEvaluation.py script. -`python code/predict_vitals.py --video_path [VIDEO_PATH]` +### Requirements: +Predictions of the models, saved as a .txt file with the names: `*GC.txt`, `*ICA_POH.txt`, `*CHROM.txt` -The default video sampling rate is 30Hz. - -#### Note - -During the inference, the program will generate a sample pre-processed frame. Please ensure it is in portrait orientation. If not, you can comment out line 30 (rotation) in the `inference_preprocess.py`. +They are located in the same folder as the ground truth files. ## Requirements Tensorflow 2.0+ +tested with Tensorflow-gpu=2.3 +`conda create -n tf-gpu tensorflow-gpu cudatoolkit=10.1` -- this command takes care of both CUDA and TF environments. -`conda create -n tf-gpu tensorflow-gpu cudatoolkit=10.1` -- this command takes care of both CUDA and TF environments. - -`pip install opencv-python scipy numpy matplotlib` +`pip install opencv-python scipy numpy matplotlib heartpy scikit-learn` -If`pip install opencv-python` does not work, I found these commands always work on my mac. +If`pip install opencv-python` does not work, I found these commands always work on my mac. ``` conda install -c menpo opencv -y @@ -88,15 +112,10 @@ pip install opencv-python ``` - +## Basis Paper +The code is based on the following paper: +#### [Xin Liu](https://homes.cs.washington.edu/~xliu0/), [Josh Fromm](https://www.linkedin.com/in/josh-fromm-2a4a2258/), [Shwetak Patel](https://ubicomplab.cs.washington.edu/members/), [Daniel McDuff](https://www.microsoft.com/en-us/research/people/damcduff/), “Multi-Task Temporal Shift Attention Networks for On-Device Contactless Vitals Measurement”, NeurIPS 2020, Oral Presentation (105 out of 9454 submissions)´ ## Contact -Please post your technical questions regarding this repo via Github Issues. - - - - - - - +Please post your technical questions regarding this repo via Github Issues. diff --git a/code/custom_fit.py b/code/custom_fit.py new file mode 100644 index 0000000..2817d70 --- /dev/null +++ b/code/custom_fit.py @@ -0,0 +1,184 @@ +import tensorflow as tf +from tensorflow import keras +from tensorflow.python.keras.engine import data_adapter +from tensorflow.python.eager import backprop +from tensorflow.python.keras.mixed_precision.experimental import loss_scale_optimizer as lso +from tensorflow.python.distribute import parameter_server_strategy + + +class CustomModel(keras.Model): + def train_step(self, data): + # Unpack the data. Its structure depends on your model and + # on what you pass to `fit()`. + data = data_adapter.expand_1d(data) + x, y, sample_weight = data_adapter.unpack_x_y_sample_weight(data) + + # with tf.GradientTape() as tape: + # y_pred = self(x, training=True) # Forward pass + # # Compute the loss value + # # (the loss function is configured in `compile()`) + # loss = self.compiled_loss(y, y_pred, regularization_losses=self.losses) + + # # Compute gradients + # trainable_vars = self.trainable_variables + # gradients = tape.gradient(loss, trainable_vars) + # # Update weights + # self.optimizer.apply_gradients(zip(gradients, trainable_vars)) + # # Update metrics (includes the metric that tracks the loss) + # self.compiled_metrics.update_state(y, y_pred) + # # Return a dict mapping metric names to current value + # return {m.name: m.result() for m in self.metrics} + + with backprop.GradientTape() as tape: + y_pred = self(x, training=True) + + # y_pred = get_peaks(y_pred) + # y = get_peaks(y) + # y, y_pred = filt_peaks(y, y_pred) + # y = tf.cast(y, tf.float32) + # y_pred = tf.cast(y_pred, tf.float32) + + loss = self.compiled_loss( + y, y_pred, sample_weight, regularization_losses=self.losses) + # For custom training steps, users can just write: + # trainable_variables = self.trainable_variables + # gradients = tape.gradient(loss, trainable_variables) + # self.optimizer.apply_gradients(zip(gradients, trainable_variables)) + # The _minimize call does a few extra steps unnecessary in most cases, + # such as loss scaling and gradient clipping. + _minimize(self.distribute_strategy, tape, self.optimizer, loss, + self.trainable_variables) + + self.compiled_metrics.update_state(y, y_pred, sample_weight) + return {m.name: m.result() for m in self.metrics} + + +def _minimize(strategy, tape, optimizer, loss, trainable_variables): + """Minimizes loss for one step by updating `trainable_variables`. + + This is roughly equivalent to + + ```python + gradients = tape.gradient(loss, trainable_variables) + self.optimizer.apply_gradients(zip(gradients, trainable_variables)) + ``` + + However, this function also applies gradient clipping and loss scaling if the + optimizer is a LossScaleOptimizer. + + Args: + strategy: `tf.distribute.Strategy`. + tape: A gradient tape. The loss must have been computed under this tape. + optimizer: The optimizer used to minimize the loss. + loss: The loss tensor. + trainable_variables: The variables that will be updated in order to minimize + the loss. + """ + + with tape: + if isinstance(optimizer, lso.LossScaleOptimizer): + loss = optimizer.get_scaled_loss(loss) + + gradients = tape.gradient(loss, trainable_variables) + + # Whether to aggregate gradients outside of optimizer. This requires support + # of the optimizer and doesn't work with ParameterServerStrategy and + # CentralStroageStrategy. + aggregate_grads_outside_optimizer = ( + optimizer._HAS_AGGREGATE_GRAD and # pylint: disable=protected-access + not isinstance(strategy.extended, + parameter_server_strategy.ParameterServerStrategyExtended)) + + if aggregate_grads_outside_optimizer: + # We aggregate gradients before unscaling them, in case a subclass of + # LossScaleOptimizer all-reduces in fp16. All-reducing in fp16 can only be + # done on scaled gradients, not unscaled gradients, for numeric stability. + gradients = optimizer._aggregate_gradients(zip(gradients, # pylint: disable=protected-access + trainable_variables)) + if isinstance(optimizer, lso.LossScaleOptimizer): + gradients = optimizer.get_unscaled_gradients(gradients) + gradients = optimizer._clip_gradients(gradients) # pylint: disable=protected-access + if trainable_variables: + if aggregate_grads_outside_optimizer: + optimizer.apply_gradients( + zip(gradients, trainable_variables), + experimental_aggregate_gradients=False) + else: + optimizer.apply_gradients(zip(gradients, trainable_variables)) + +@tf.function +def get_peaks(y): + # y: (N,) + data_reshaped = tf.reshape(y, (1, -1, 1)) # (1, N, 1) + max_pooled_in_tensor = tf.nn.max_pool(data_reshaped, (20,), 1,'SAME') + maxima = tf.equal(data_reshaped,max_pooled_in_tensor) # (1, N, 1) + maxima = tf.cast(maxima, tf.float32) + maxima = tf.squeeze(maxima) # (N,1) + peaks = tf.where(maxima) # now only the Peak Indices (A, 3) + peaks = tf.reshape(peaks, (tf.size(y),)) # (A,1) + + return peaks + +# x: true y: prediction +# input: peaks of truth and prediction as tensor... +@tf.function +def filt_peaks(x,y): + def true_fn(): + return min + def false_fn(): + return tf.cast(-1, tf.int64) + max_offset = 10 + mask = tf.cast(tf.zeros(tf.size(x)),tf.bool) # tensor with size of x (truth data) + # check which peaks of truth are recognized in pred + min = 0 + min = tf.cast(min, tf.int64) + + # for item in y: # items of predicion + # diff = tf.abs(x - item) # diff of truth data and item + # min = tf.reduce_min(diff) # minimum of diff + # min = tf.cond(tf.less(min, max_offset), true_fn, false_fn) + # temp_mask = tf.equal(min, diff) + # mask = tf.logical_or(mask, temp_mask) + + # x = tf.boolean_mask(x, mask) + def fn(item): + def true_fn(): + return tf.cast(min, tf.float64) + def false_fn(): + return tf.cast(-1, tf.float64) + diff = tf.abs(x - item) # diff of truth data and item + diff = tf.cast(diff, tf.float64) + min = tf.reduce_min(diff) # minimum of diff + min = tf.cond(tf.less(min, max_offset), true_fn, false_fn) + temp_mask = tf.equal(min, diff) + return temp_mask + mask1 = tf.map_fn(fn=lambda item: fn(item), elems=y, fn_output_signature=tf.bool) + mask1 = tf.reduce_any(mask1, 0) + x = tf.boolean_mask(x,mask1) + + # check if outliners are in pred + # mask = tf.cast(tf.zeros(tf.size(y)), tf.bool) + # for item in x: + # diff = tf.abs(y - item) # diff of truth data and item + # min = tf.reduce_min(diff) # minimum of diff + # min = tf.cond(tf.less(min, max_offset), true_fn, false_fn) + # temp_mask = tf.equal(min, diff) + # mask = tf.logical_or(mask, temp_mask) + # y = tf.boolean_mask(y,mask) + + def fn2(item): + def true_fn(): + return tf.cast(min, tf.float64) + def false_fn(): + return tf.cast(-1, dtype=tf.float64) + diff = tf.abs(y - item) # diff of truth data and item + diff = tf.cast(diff, tf.float64) + min = tf.reduce_min(diff) # minimum of diff + min = tf.cond(tf.less(min, max_offset), true_fn, false_fn) + temp_mask = tf.equal(min, diff) + return temp_mask + mask2 = tf.map_fn(fn=lambda item: fn2(item), elems=x, fn_output_signature=tf.bool) + mask2 = tf.reduce_any(mask2, 0) + y = tf.boolean_mask(y,mask2) + + return x, y \ No newline at end of file diff --git a/code/data_generator.py b/code/data_generator.py index b4cbdfc..6c3fa95 100755 --- a/code/data_generator.py +++ b/code/data_generator.py @@ -1,32 +1,40 @@ ''' Data Generator for Multi-Task Temporal Shift Attention Networks for On-Device Contactless Vitals Measurement Author: Xin Liu + +Further Development: Sarah Quehl ''' import math import h5py import numpy as np -from tensorflow import keras - +from pandas.core.resample import h +import tensorflow as tf +from tensorflow.python.keras.utils import data_utils +import ast -class DataGenerator(keras.utils.Sequence): +class DataGenerator(data_utils.Sequence): 'Generates data for Keras' - def __init__(self, paths_of_videos, nframe_per_video, dim, batch_size=32, frame_depth=10, - shuffle=True, temporal=True, respiration=0): + def __init__(self, paths_of_videos, maxLen_Video, dim, batch_size=32, frame_depth=10, + shuffle=True, temporal=True, respiration=0, database_name = None, time_error_loss=False, truth_parameter=None): self.dim = dim self.batch_size = batch_size self.paths_of_videos = paths_of_videos - self.nframe_per_video = nframe_per_video + self.maxLen_Video = maxLen_Video self.shuffle = shuffle self.temporal = temporal self.frame_depth = frame_depth self.respiration = respiration + self.database_name = database_name + self.time_error_loss = time_error_loss + self.truth_parameter = truth_parameter self.on_epoch_end() - def __len__(self): + def __len__(self): 'Denotes the number of batches per epoch' - return math.ceil(len(self.paths_of_videos) / self.batch_size) + temp_var = math.ceil(len(self.paths_of_videos) / self.batch_size) + return temp_var def __getitem__(self, index): 'Generate one batch of data' @@ -44,20 +52,255 @@ def on_epoch_end(self): def __data_generation(self, list_video_temp): 'Generates data containing batch_size samples' - if self.respiration == 1: - label_key = "drsub" - else: - label_key = 'dysub' - if self.temporal == 'CAN_3D': - num_window = self.nframe_per_video - (self.frame_depth + 1) + sum_frames_batch = get_frame_sum_3D_Hybrid(list_video_temp, self.maxLen_Video) + data = np.zeros((sum_frames_batch, self.dim[0], self.dim[1],self.frame_depth, 6), dtype=np.float32) + label = np.zeros((sum_frames_batch, self.frame_depth), dtype=np.float32) + index_counter = 0 + for index, temp_path in enumerate(list_video_temp): + f1 = h5py.File(temp_path, 'r') + dXsub = np.array(f1['data']) + dysub = np.array(f1['pulse']) + # if dXsub.shape[0] > self.maxLen_Video: # only 30 sek videos + # dXsub = dXsub[0:self.maxLen_Video, :,:,:] + # dysub = dysub[0:self.maxLen_Video] + num_window = int(dXsub.shape[0]) -(self.frame_depth+1) + tempX = np.array([dXsub[f:f + self.frame_depth, :, :, :] # (491, 10, 36, 36 ,6) (169, 10, 36, 36, 6) + for f in range(num_window)]) + tempY = np.array([dysub[f:f + self.frame_depth] #(491,10,1) - (169, 10, 1) + for f in range(num_window)]) + tempX = np.swapaxes(tempX, 1, 3) # (169, 36, 36, 10, 6) + tempX = np.swapaxes(tempX, 1, 2) # (169, 36, 36, 10, 6) + tempY = np.reshape(tempY, (num_window, self.frame_depth)) # (169, 10) + data[index_counter: index_counter + num_window, :, :, :, :] = tempX + label[index_counter: index_counter + num_window, :] = tempY + index_counter += num_window + + motion_data = data[:, :, :, :, :3] + apperance_data = data[:, :, :, :, -3:] + max_data = num_window*self.frame_depth + motion_data = motion_data[0:max_data, :, :, :] + apperance_data = apperance_data[0:max_data, :, :, :] + label = label[0:max_data, :] + output = (motion_data, apperance_data) + + elif self.temporal == 'CAN': + sum_frames_batch = get_frame_sum(list_video_temp, self.maxLen_Video) + data = np.zeros((sum_frames_batch, self.dim[0], self.dim[1], 6), dtype=np.float32) + label = np.zeros((sum_frames_batch, 1), dtype=np.float32) + num_window = int(sum_frames_batch/ self.frame_depth) + index_counter = 0 + for index, temp_path in enumerate(list_video_temp): + f1 = h5py.File(temp_path, 'r') + dXsub = np.array(f1['data']) + dysub = np.array(f1['pulse']) + # if dXsub.shape[0] > self.maxLen_Video: # only 1 min videos + # current_nframe = self.maxLen_Video + # dXsub = dXsub[0:self.maxLen_Video, :,:,:] + # dysub = dysub[0:self.maxLen_Video] + # else: + current_nframe = dXsub.shape[0] + data[index_counter:index_counter+current_nframe, :, :, :] = dXsub + label[index_counter:index_counter+current_nframe, 0] = dysub # data BVP + index_counter += current_nframe + motion_data = data[:, :, :, :3] + apperance_data = data[:, :, :, -3:] + max_data = num_window*self.frame_depth + motion_data = motion_data[0:max_data, :, :, :] + apperance_data = apperance_data[0:max_data, :, :, :] + label = label[0:max_data, 0] + + output = (motion_data, apperance_data) + + elif self.temporal == 'TS_CAN': + sum_frames_batch = get_frame_sum(list_video_temp, self.maxLen_Video) + data = np.zeros((sum_frames_batch, self.dim[0], self.dim[1], 6), dtype=np.float32) + label = np.zeros((sum_frames_batch, 1), dtype=np.float32) + num_window = int(sum_frames_batch/ self.frame_depth) + index_counter = 0 + for index, temp_path in enumerate(list_video_temp): + f1 = h5py.File(temp_path, 'r') + dXsub = np.array(f1['data']) + dysub = np.array(f1['pulse']) + # if dXsub.shape[0] > self.maxLen_Video: # only 1 min videos + # current_nframe = self.maxLen_Video + # dXsub = dXsub[0:self.maxLen_Video, :,:,:] + # dysub = dysub[0:self.maxLen_Video] + # else: + current_nframe = dXsub.shape[0] + data[index_counter:index_counter+current_nframe, :, :, :] = dXsub + label[index_counter:index_counter+current_nframe, 0] = dysub # data BVP + index_counter += current_nframe + motion_data = data[:, :, :, :3] + apperance_data = data[:, :, :, -3:] + + if num_window % 2 == 1: + num_window = num_window - 1 + max_data = num_window * self.frame_depth + else: + max_data = num_window*self.frame_depth + motion_data = motion_data[0:max_data, :, :, :] + apperance_data = apperance_data[0:max_data, :, :, :] + label = label[0:max_data, 0] + apperance_data = np.reshape(apperance_data, (num_window, self.frame_depth, self.dim[0], self.dim[1], 3)) + apperance_data = np.average(apperance_data, axis=1) + apperance_data = np.repeat(apperance_data[:, np.newaxis, :, :, :], self.frame_depth, axis=1) + apperance_data = np.reshape(apperance_data, (apperance_data.shape[0] * apperance_data.shape[1], + apperance_data.shape[2], apperance_data.shape[3], + apperance_data.shape[4])) + output = (motion_data, apperance_data) + + # new Peak Temperal Shift CAN + elif self.temporal == 'PTS_CAN': + sum_frames_batch = get_frame_sum(list_video_temp, self.maxLen_Video) + data = np.zeros((sum_frames_batch, self.dim[0], self.dim[1], 6), dtype=np.float32) + label_y = np.zeros((sum_frames_batch, 1), dtype=np.float32) + label_z = np.zeros((sum_frames_batch, 1), dtype=np.float32) + num_window = int(sum_frames_batch/ self.frame_depth) + index_counter = 0 + for index, temp_path in enumerate(list_video_temp): + f1 = h5py.File(temp_path, 'r') + dXsub = np.array(f1['data']) + dysub = np.array(f1['pulse']) + dzsub = np.array(f1['peaklist']) + + if dXsub.shape[0] > self.maxLen_Video: # UBFC-PHYS + current_nframe = dXsub.shape[0] + sigma = 2.6 + fps = 35.138 + elif dXsub.shape[0] > 1300 and dXsub.shape[0] < 2200: # UBFC-rPPG + current_nframe = dXsub.shape[0] + sigma = 2.2 + fps = 29.51 + else: #COHFACE + current_nframe = dXsub.shape[0] + sigma = 1.5 + fps = 20 + data[index_counter:index_counter+current_nframe, :, :, :] = dXsub + label_y[index_counter:index_counter+current_nframe, 0] = dysub # data BVP + if(self.time_error_loss == False): + temp = gauss_loss_dataGenerator(current_nframe, dzsub, sigma) + else: + temp = time_error_loss_dataGenerator(current_nframe, dzsub, fps) + label_z[index_counter:index_counter+current_nframe, 0] = temp # data Peaks + index_counter += current_nframe + motion_data = data[:, :, :, :3] + apperance_data = data[:, :, :, -3:] + + if num_window % 2 == 1: + num_window = num_window - 1 + max_data = num_window * self.frame_depth + else: + max_data = num_window*self.frame_depth + motion_data = motion_data[0:max_data, :, :, :] + apperance_data = apperance_data[0:max_data, :, :, :] + label_y = label_y[0:max_data, 0] + label_z = label_z[0:max_data, 0] + apperance_data = np.reshape(apperance_data, (num_window, self.frame_depth, self.dim[0], self.dim[1], 3)) + apperance_data = np.average(apperance_data, axis=1) + apperance_data = np.repeat(apperance_data[:, np.newaxis, :, :, :], self.frame_depth, axis=1) + apperance_data = np.reshape(apperance_data, (apperance_data.shape[0] * apperance_data.shape[1], + apperance_data.shape[2], apperance_data.shape[3], + apperance_data.shape[4])) + label = (label_y, label_z) + output = (motion_data, apperance_data) + + # new Parameter Peak Temperal Shift CAN + elif self.temporal == 'PPTS_CAN': + sum_frames_batch = get_frame_sum(list_video_temp, self.maxLen_Video) + data = np.zeros((sum_frames_batch, self.dim[0], self.dim[1], 6), dtype=np.float32) + label_y = np.zeros((sum_frames_batch, 1), dtype=np.float32) + label_z = np.zeros((sum_frames_batch, 1), dtype=np.float32) + label_params = np.zeros(self.batch_size*len(self.truth_parameter), dtype=np.float32) + + num_window = int(sum_frames_batch/ self.frame_depth) + index_counter = 0 + param_counter = 0 + for index, temp_path in enumerate(list_video_temp): + f1 = h5py.File(temp_path, 'r') + dXsub = np.array(f1['data']) + dysub = np.array(f1['pulse']) + dzsub = np.array(f1['peaklist']) + + truthParams = np.array(f1['parameter']) + truthParams = np.array(truthParams) + truthParams = ast.literal_eval(str(truthParams)) + params = np.zeros(len(self.truth_parameter)) + for parameter_index in range(0,len(self.truth_parameter)): + if str(self.truth_parameter[parameter_index]) == "lf_hf": + nn_list = np.array(f1['nn']) + nn_list = tf.reshape(tf.convert_to_tensor(nn_list), (-1,)) + frq = tf.cast(tf.abs(tf.signal.rfft(nn_list)), tf.float32)/tf.cast(tf.size(nn_list),tf.float32) + frq = tf.multiply(tf.pow(frq,2),tf.math.sqrt(tf.cast(2, tf.float32))) + + dt = tf.math.reduce_mean(nn_list) / 1000 # in sec + t = tf.cast(tf.range(0, tf.size(frq)), tf.float32) + t = tf.cast(t, tf.float32)/(tf.cast(dt, tf.float32)*tf.cast(tf.size(frq)*2, tf.float32)) + + mask_lf = tf.cast(tf.logical_and(tf.greater_equal(t, 0.04), tf.less(t, 0.15)), tf.float32) + lf = tf.maximum(tf.reduce_sum(frq*mask_lf), 0.000001) + mask_hf = tf.cast(tf.logical_and(tf.greater_equal(t, 0,15), tf.less(t, 0.4)), tf.float32) + hf = tf.maximum(tf.reduce_sum(frq*mask_hf), 0.000001) + + lf_hf = lf/hf + params[parameter_index] = lf_hf + else: + params[parameter_index] = truthParams[str(self.truth_parameter[parameter_index])] + label_params[param_counter: param_counter+len(self.truth_parameter)] = params + + if dXsub.shape[0] > self.maxLen_Video: # UBFC-PHYS + current_nframe = dXsub.shape[0] + sigma = 2.6 + fps = 35.138 + else: #COHFACE + current_nframe = dXsub.shape[0] + sigma = 1.5 + fps = 20 + data[index_counter:index_counter+current_nframe, :, :, :] = dXsub + label_y[index_counter:index_counter+current_nframe, 0] = dysub # data BVP + if(self.time_error_loss == False): + temp = gauss_loss_dataGenerator(current_nframe, dzsub, sigma) + else: + temp = time_error_loss_dataGenerator(current_nframe, dzsub, fps) + label_z[index_counter:index_counter+current_nframe, 0] = temp # data Peaks + index_counter += current_nframe + param_counter += len(self.truth_parameter) + motion_data = data[:, :, :, :3] + apperance_data = data[:, :, :, -3:] + + if num_window % 2 == 1: + num_window = num_window - 1 + max_data = num_window * self.frame_depth + else: + max_data = num_window*self.frame_depth + motion_data = motion_data[0:max_data, :, :, :] + apperance_data = apperance_data[0:max_data, :, :, :] + label_y = label_y[0:max_data, 0] + label_z = label_z[0:max_data, 0] + apperance_data = np.reshape(apperance_data, (num_window, self.frame_depth, self.dim[0], self.dim[1], 3)) + apperance_data = np.average(apperance_data, axis=1) + apperance_data = np.repeat(apperance_data[:, np.newaxis, :, :, :], self.frame_depth, axis=1) + apperance_data = np.reshape(apperance_data, (apperance_data.shape[0] * apperance_data.shape[1], + apperance_data.shape[2], apperance_data.shape[3], + apperance_data.shape[4])) + label = (label_y, label_z, label_params) + output = (motion_data, apperance_data) + + + elif self.temporal == 'Hybrid_CAN': + sum_frames_batch = get_frame_sum_3D_Hybrid(list_video_temp, self.maxLen_Video) data = np.zeros((num_window*len(list_video_temp), self.dim[0], self.dim[1], self.frame_depth, 6), dtype=np.float32) label = np.zeros((num_window*len(list_video_temp), self.frame_depth), dtype=np.float32) + index_counter = 0 for index, temp_path in enumerate(list_video_temp): f1 = h5py.File(temp_path, 'r') - dXsub = np.transpose(np.array(f1["dXsub"])) - dysub = np.array(f1[label_key]) + dXsub = np.array(f1['data']) + dysub = np.array(f1['pulse']) + # if dXsub.shape[0] > self.maxLen_Video: # only 30 sek videos + # dXsub = dXsub[0:self.maxLen_Video, :,:,:] + # dysub = dysub[0:self.maxLen_Video] + num_window = int(dXsub.shape[0]) -(self.frame_depth+1) tempX = np.array([dXsub[f:f + self.frame_depth, :, :, :] # (169, 10, 36, 36, 6) for f in range(num_window)]) tempY = np.array([dysub[f:f + self.frame_depth] # (169, 10, 1) @@ -65,9 +308,29 @@ def __data_generation(self, list_video_temp): tempX = np.swapaxes(tempX, 1, 3) # (169, 36, 36, 10, 6) tempX = np.swapaxes(tempX, 1, 2) # (169, 36, 36, 10, 6) tempY = np.reshape(tempY, (num_window, self.frame_depth)) # (169, 10) - data[index*num_window:(index+1)*num_window, :, :, :, :] = tempX - label[index*num_window:(index+1)*num_window, :] = tempY - output = (data[:, :, :, :, :3], data[:, :, :, :, -3:]) + data[index_counter: index_counter + num_window, :, :, :, :] = tempX + label[index_counter: index_counter + num_window, :] = tempY + index_counter += num_window + motion_data = data[:, :, :, :, :3] + apperance_data = np.average(data[:, :, :, :, -3:], axis=-2) + output = (motion_data, apperance_data) + + # Multi-Task Approaches with Respiration rate + elif self.temporal == 'MT_CAN': + data = np.zeros((self.nframe_per_video * len(list_video_temp), self.dim[0], self.dim[1], 6), + dtype=np.float32) + label_y = np.zeros((self.nframe_per_video * len(list_video_temp), 1), dtype=np.float32) + label_r = np.zeros((self.nframe_per_video * len(list_video_temp), 1), dtype=np.float32) + for index, temp_path in enumerate(list_video_temp): + f1 = h5py.File(temp_path, 'r') + dXsub = np.transpose(np.array(f1["dXsub"])) # dRsub for respiration + drsub = np.array(f1['drsub']) + dysub = np.array(f1['dysub']) + data[index * self.nframe_per_video:(index + 1) * self.nframe_per_video, :, :, :] = dXsub + label_y[index*self.nframe_per_video:(index+1)*self.nframe_per_video, :] = dysub + label_r[index * self.nframe_per_video:(index + 1) * self.nframe_per_video, :] = drsub + output = (data[:, :, :, :3], data[:, :, :, -3:]) + label = (label_y, label_r) elif self.temporal == 'MT_CAN_3D': num_window = self.nframe_per_video - (self.frame_depth + 1) data = np.zeros((num_window*len(list_video_temp), self.dim[0], self.dim[1], self.frame_depth, 6), @@ -94,65 +357,30 @@ def __data_generation(self, list_video_temp): label_r[index * num_window:(index + 1) * num_window, :] = tempY_r output = (data[:, :, :, :, :3], data[:, :, :, :, -3:]) label = (label_y, label_r) - elif self.temporal == 'CAN': - data = np.zeros((self.nframe_per_video * len(list_video_temp), self.dim[0], self.dim[1], 6), dtype=np.float32) - label = np.zeros((self.nframe_per_video * len(list_video_temp), 1), dtype=np.float32) - for index, temp_path in enumerate(list_video_temp): - f1 = h5py.File(temp_path, 'r') - dXsub = np.transpose(np.array(f1["dXsub"])) #dRsub for respiration - dysub = np.array(f1[label_key]) - data[index*self.nframe_per_video:(index+1)*self.nframe_per_video, :, :, :] = dXsub - label[index*self.nframe_per_video:(index+1)*self.nframe_per_video, :] = dysub - output = (data[:, :, :, :3], data[:, :, :, -3:]) - elif self.temporal == 'MT_CAN': - data = np.zeros((self.nframe_per_video * len(list_video_temp), self.dim[0], self.dim[1], 6), - dtype=np.float32) - label_y = np.zeros((self.nframe_per_video * len(list_video_temp), 1), dtype=np.float32) - label_r = np.zeros((self.nframe_per_video * len(list_video_temp), 1), dtype=np.float32) - for index, temp_path in enumerate(list_video_temp): - f1 = h5py.File(temp_path, 'r') - dXsub = np.transpose(np.array(f1["dXsub"])) # dRsub for respiration - drsub = np.array(f1['drsub']) - dysub = np.array(f1['dysub']) - data[index * self.nframe_per_video:(index + 1) * self.nframe_per_video, :, :, :] = dXsub - label_y[index*self.nframe_per_video:(index+1)*self.nframe_per_video, :] = dysub - label_r[index * self.nframe_per_video:(index + 1) * self.nframe_per_video, :] = drsub - output = (data[:, :, :, :3], data[:, :, :, -3:]) - label = (label_y, label_r) - elif self.temporal == 'TS_CAN': - data = np.zeros((self.nframe_per_video * len(list_video_temp), self.dim[0], self.dim[1], 6), dtype=np.float32) - label = np.zeros((self.nframe_per_video * len(list_video_temp), 1), dtype=np.float32) - num_window = int(self.nframe_per_video / self.frame_depth) * len(list_video_temp) - for index, temp_path in enumerate(list_video_temp): - f1 = h5py.File(temp_path, 'r') - dXsub = np.transpose(np.array(f1["dXsub"])) #dRsub for respiration - dysub = np.array(f1[label_key]) - data[index*self.nframe_per_video:(index+1)*self.nframe_per_video, :, :, :] = dXsub - label[index*self.nframe_per_video:(index+1)*self.nframe_per_video, :] = dysub - motion_data = data[:, :, :, :3] - apperance_data = data[:, :, :, -3:] - apperance_data = np.reshape(apperance_data, (num_window, self.frame_depth, self.dim[0], self.dim[1], 3)) - apperance_data = np.average(apperance_data, axis=1) - apperance_data = np.repeat(apperance_data[:, np.newaxis, :, :, :], self.frame_depth, axis=1) - apperance_data = np.reshape(apperance_data, (apperance_data.shape[0] * apperance_data.shape[1], - apperance_data.shape[2], apperance_data.shape[3], - apperance_data.shape[4])) - output = (motion_data, apperance_data) elif self.temporal == 'MTTS_CAN': - data = np.zeros((self.nframe_per_video * len(list_video_temp), self.dim[0], self.dim[1], 6), dtype=np.float32) - label_y = np.zeros((self.nframe_per_video * len(list_video_temp), 1), dtype=np.float32) - label_r = np.zeros((self.nframe_per_video * len(list_video_temp), 1), dtype=np.float32) - num_window = int(self.nframe_per_video / self.frame_depth) * len(list_video_temp) + sum_frames_batch = get_frame_sum(list_video_temp, self.maxLen_Video) + data = np.zeros((sum_frames_batch, self.dim[0], self.dim[1], 6), dtype=np.float32) + label_y = np.zeros((sum_frames_batch, 1), dtype=np.float32) + label_r = np.zeros((sum_frames_batch, 1), dtype=np.float32) + num_window = int(sum_frames_batch/ self.frame_depth) + index_counter = 0 for index, temp_path in enumerate(list_video_temp): f1 = h5py.File(temp_path, 'r') - dXsub = np.transpose(np.array(f1["dXsub"])) #dRsub for respiration - drsub = np.array(f1['drsub']) - dysub = np.array(f1['dysub']) - data[index*self.nframe_per_video:(index+1)*self.nframe_per_video, :, :, :] = dXsub - label_y[index*self.nframe_per_video:(index+1)*self.nframe_per_video, :] = dysub - label_r[index * self.nframe_per_video:(index + 1) * self.nframe_per_video, :] = drsub + dXsub = np.array(f1['data']) + drsub = np.array(f1['respiration']) + dysub = np.array(f1['pulse']) + current_nframe = dXsub.shape[0] + data[index_counter:index_counter+current_nframe, :, :, :] = dXsub + label_y[index_counter:index_counter+current_nframe, 0] = dysub # data BVP + label_r[index_counter:index_counter+current_nframe, 0] = drsub # data Respiration + index_counter += current_nframe motion_data = data[:, :, :, :3] apperance_data = data[:, :, :, -3:] + max_data = num_window*self.frame_depth + motion_data = motion_data[0:max_data, :, :, :] + apperance_data = apperance_data[0:max_data, :, :, :] + label_y = label_y[0:max_data, 0] + label_r = label_r[0:max_data, 0] apperance_data = np.reshape(apperance_data, (num_window, self.frame_depth, self.dim[0], self.dim[1], 3)) apperance_data = np.average(apperance_data, axis=1) apperance_data = np.repeat(apperance_data[:, np.newaxis, :, :, :], self.frame_depth, axis=1) @@ -189,28 +417,80 @@ def __data_generation(self, list_video_temp): apperance_data = np.average(data[:, :, :, :, -3:], axis=-2) output = (motion_data, apperance_data) label = (label_y, label_r) - elif self.temporal == 'Hybrid_CAN': - num_window = self.nframe_per_video - (self.frame_depth + 1) - data = np.zeros((num_window*len(list_video_temp), self.dim[0], self.dim[1], self.frame_depth, 6), - dtype=np.float32) - label = np.zeros((num_window*len(list_video_temp), self.frame_depth), dtype=np.float32) - for index, temp_path in enumerate(list_video_temp): - f1 = h5py.File(temp_path, 'r') - dXsub = np.transpose(np.array(f1["dXsub"])) - dysub = np.array(f1[label_key]) - tempX = np.array([dXsub[f:f + self.frame_depth, :, :, :] # (169, 10, 36, 36, 6) - for f in range(num_window)]) - tempY = np.array([dysub[f:f + self.frame_depth] # (169, 10, 1) - for f in range(num_window)]) - tempX = np.swapaxes(tempX, 1, 3) # (169, 36, 36, 10, 6) - tempX = np.swapaxes(tempX, 1, 2) # (169, 36, 36, 10, 6) - tempY = np.reshape(tempY, (num_window, self.frame_depth)) # (169, 10) - data[index*num_window:(index+1)*num_window, :, :, :, :] = tempX - label[index*num_window:(index+1)*num_window, :] = tempY - motion_data = data[:, :, :, :, :3] - apperance_data = np.average(data[:, :, :, :, -3:], axis=-2) - output = (motion_data, apperance_data) else: raise ValueError('Unsupported Model!') return output, label + +def find_csv(video_path): + csv_path = str(video_path).replace("vid", "bvp").replace(".avi", ".csv") + return csv_path + +def get_frame_sum(list_vid, maxLen_Video): + frames_sum = 0 + counter = 0 + for vid in list_vid: + hf = h5py.File(vid, 'r') + shape = hf['data'].shape + # if shape[0] > maxLen_Video: + # frames_sum += maxLen_Video + # else: + frames_sum += shape[0] + counter += 1 + return frames_sum + +def get_frame_sum_3D_Hybrid(list_vid, maxLen_Video): + frames_sum = 0 + counter = 0 + for vid in list_vid: + hf = h5py.File(vid, 'r') + shape = hf['data'].shape + # if shape[0] > maxLen_Video: + # frames_sum += maxLen_Video - 9 + # else: + frames_sum += shape[0] - 9 + counter += 1 + return frames_sum + +def gauss_loss_dataGenerator(current_nframe, dzsub, sigma): + temp = np.zeros(current_nframe, dtype=np.float32) + for i in dzsub: + mu = i + min = int(i-sigma*3) + if min < 0: + min = 0 + max = int(i+sigma*3) + if max > len(temp): + max = len(temp)-1 + + for j in range(min, max): + temp[j] = gauss(j, sigma, mu) + return temp + +def time_error_loss_dataGenerator(current_nframe, dzsub, fps): + temp = np.zeros(current_nframe, dtype=np.float32) + m = 1/fps + for i in range(0, len(dzsub)): + peak_1 = dzsub[i] + if(i-1 >= 0): + peak_0 = dzsub[i-1] + min = int(round((peak_1 - peak_0)/2) + peak_0 + 1) + for j in range(min, peak_1+1): + temp[j] = m*(peak_1 - j) + elif(i-1 == -1): + min = 0 + for j in range(min, peak_1+1): + temp[j] = m*(peak_1 - j) + if(i+1 < len(dzsub)): + peak_2 = dzsub[i+1] + max = int(round((peak_2 - peak_1)/2) + peak_1) + for j in range(peak_1, max+1): + temp[j] = m*(j-peak_1) + elif(i+1 == len(dzsub)): + max = len(temp)-1 + for j in range(peak_1, max+1): + temp[j] = m*(j-peak_1) + return temp + +def gauss(x, sigma, mu): + return math.exp(-(x - mu)**2 / (2 * sigma**2)) / (sigma * math.sqrt(2 * math.pi)) diff --git a/code/evaluation_iPhys.py b/code/evaluation_iPhys.py new file mode 100644 index 0000000..5f91e2c --- /dev/null +++ b/code/evaluation_iPhys.py @@ -0,0 +1,246 @@ +import numpy as np +import scipy.io +import xlsxwriter +import h5py +import os +import matplotlib.pyplot as plt +from scipy.signal import butter +from inference_preprocess import preprocess_raw_frames, preprocess_raw_video, detrend +from sklearn.preprocessing import MinMaxScaler +from sklearn import metrics +import scipy.stats as sc +from glob import glob +from scipy import signal + +import heartpy as hp + + +def write_header(worksheet): + header = ['Database', 'Subj/Task', 'HR-pred', 'HR-truth', + 'p','meanNN-pred', 'meanNN-truth', 'sdnn-pred', 'sdnn-truth', + 'rmssd-pred', 'rmssd-truth', 'pNN50-pred', 'pNN50-truth', + 'LF-pred', 'LF-truth', 'HF-pred', 'HF-truth', + 'TP-pred', 'TP-truth', 'LF/HF-pred', 'LF/HF-truth', + 'sd1-pred', 'sd1_truth', 'sd2_pred', 'sd2_truth', 'MAE'] + for index in range(len(header)): + worksheet.write(0,index, header[index]) + +def predict_vitals(worksheet, video_path, save_dir): + mms = MinMaxScaler() + + ###### load video Data ####### + counter_video = 1 + old_database = "COH" + for sample_data_path in video_path: + + print("path: ",sample_data_path) + pulse_pred = open(sample_data_path, 'r').read() + pulse_pred = str(pulse_pred).split("\n") + pulse_pred = pulse_pred[:-1] + pulse_pred = np.array(list(map(float, pulse_pred))) + + mean = pulse_pred.mean() + std = np.std(pulse_pred) + upper_limit = mean + std*3 + lower_limit = mean - std*3 + for x in range(0, len(pulse_pred)): + if pulse_pred[x] > upper_limit: + pulse_pred[x] = upper_limit + elif pulse_pred[x] < lower_limit: + pulse_pred[x] = lower_limit + pulse_pred = np.array(mms.fit_transform(pulse_pred.reshape(-1,1))).flatten() # normalization + + + ##### ground truth data resampled ####### + data_path = "default" + if(str(sample_data_path).find("GC") > 0): + method = "GC" + data_path = sample_data_path.replace("_GC", "") + elif(str(sample_data_path).find("ICA_POH") > 0): + method = "ICA" + data_path = sample_data_path.replace("_ICA_POH", "") + elif(str(sample_data_path).find("CHROM") > 0): + method = "CHROM" + data_path = sample_data_path.replace("_CHROM", "") + else: + raise print("ERROR") + + if(str(sample_data_path).find("COHFACE") > 0): + database_name = "COH" + fs = 20 + truth_path = data_path.replace(".txt", "_dataFile.hdf5") # akutell für COHACE... + elif(str(sample_data_path).find("UBFC-PHYS") > 0): + database_name = "UB-Ph" + fs = 35 + truth_path = data_path.replace("vid_", "").replace(".txt","_dataFile.hdf5") + elif(str(sample_data_path).find("UBFC") > 0): + database_name = "UBFC" + fs = 30 + truth_path = data_path.replace("vid.txt", "dataFile.hdf5") + elif(str(sample_data_path).find("BP4D") > 0): + fs = 25 + database_name = "BP4D" + truth_path = data_path + "/BP_mmHg.txt" + else: + return print("Error in finding the ground truth signal...") + + if database_name != "BP4D": + gound_truth_file = h5py.File(truth_path, "r") + pulse_truth = gound_truth_file["pulse"] ### range ground truth from 0 to 1 + pulse_truth = pulse_truth[0:len(pulse_pred)] + else: + data = open(truth_path, 'r').read() + data = str(data).split("\n") + pulse_truth = np.array(list(map(float, data[0:-1]))) + mean = pulse_truth.mean() + std = np.std(pulse_truth) + upper_limit = mean + std*3 + lower_limit = mean - std*3 + for x in range(0, len(pulse_truth)): + if pulse_truth[x] > upper_limit: + pulse_truth[x] = upper_limit + elif pulse_truth[x] < lower_limit: + pulse_truth[x] = lower_limit + + pulse_truth = np.array(signal.resample(pulse_truth, len(pulse_pred))) + pulse_truth = np.array(mms.fit_transform(pulse_truth.reshape(-1,1))).flatten() # normalization + pulse_truth = detrend(np.cumsum(pulse_truth), 100) + [b_pulse_tr, a_pulse_tr] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') + pulse_truth = scipy.signal.filtfilt(b_pulse_tr, a_pulse_tr, np.double(pulse_truth)) + pulse_truth = np.array(mms.fit_transform(pulse_truth.reshape(-1,1))).flatten() + + ### same size ####### + if len(pulse_pred) > len(pulse_truth): + pulse_pred = pulse_pred[:len(pulse_truth)] + elif len(pulse_pred) < len(pulse_truth): + pulse_truth = pulse_truth[:len(pulse_pred)] + ########### Peaks ########### + try: + working_data_pred, measures_pred = hp.process(pulse_pred, fs, calc_freq=True) + working_data_truth, measures_truth = hp.process(pulse_truth, fs, calc_freq=True) + except: + continue + peaks_pred = working_data_pred['peaklist'] + peaks_truth = working_data_truth['peaklist'] + + ######## name files ############# + if(str(data_path).find("COHFACE") > 0): + nmr = str(data_path).find("COHFACE") + nameStr = str(data_path)[nmr + 7:].replace("\\", "-").replace("-data.txt", "") + elif(str(data_path).find("UBFC-PHYS") > 0): + nmr = str(data_path).find("UBFC-PHYS") + nameStr = str(data_path)[nmr + 12:].replace("\\", "-").replace("vid_", "").replace(".txt", "") + elif(str(data_path).find("UBFC") > 0): + nmr = str(data_path).find("UBFC") + nameStr = str(data_path)[nmr + 5:].replace("\\", "-").replace("vid.txt", "") + elif(str(data_path).find("BP4D") > 0): + nmr = str(data_path).find("BP4D") + nameStr = str(data_path)[nmr + 5:].replace("\\", "-") + else: + raise ValueError + ########## Plot ################## + peaks_pred_new = [] + for peak in peaks_pred: + if (peak > 400 and peak <700): + peaks_pred_new.append(peak-400) + peaks_truth_new = [] + for peak in peaks_truth: + if (peak > 400 and peak <700): + peaks_truth_new.append(peak-400) + plt.figure() #subplot(211) + plt.plot(pulse_pred[400:700], "#E6001A", label='rPPG signal') + plt.plot(peaks_truth_new, pulse_truth[400:700][peaks_truth_new], "x", color="#005AA9") + plt.plot(peaks_pred_new, pulse_pred[400:700][peaks_pred_new], "x", color ='#E6001A') + plt.title('rPPG signal with ground truth') + plt.ylabel("normalized Signal [a.u.]") + plt.xlabel("time (samples)") + plt.plot(pulse_truth[400:700], '#005AA9', linewidth=0.9, label='ground truth') + plt.legend() + plt.savefig(save_dir + database_name+ nameStr + method + "_both.svg", format="svg") + + plt.figure() + plt.subplot(211) + plt.plot(pulse_truth[400:700],"#004E8A", label='Ground truth') + plt.plot(peaks_truth_new, pulse_truth[400:700][peaks_truth_new], "x", color="#004E8A") + plt.ylabel("normalized Signal [a.u.]") + plt.title('Ground truth') + plt.subplot(212) + plt.plot(pulse_pred[400:700], "#004E8A",label='Prediction') + plt.plot(peaks_pred_new, pulse_pred[400:700][peaks_pred_new],"x", color="#004E8A") + plt.title("Predicted rPPG") + plt.ylabel("normalized Signal [a.u.]") + plt.xlabel("time (samples)") + plt.legend() + plt.savefig(save_dir + database_name+ nameStr + method +".svg", format="svg") + + ######### Metrics ############## + # MSE: + MAE = metrics.mean_absolute_error(pulse_truth, pulse_pred) + # RMSE: + RMSE = metrics.mean_squared_error(pulse_truth, pulse_pred, squared=False) + # Pearson correlation: + p = sc.pearsonr(pulse_truth, pulse_pred) + + ####### Logging ############# + if database_name != old_database: + counter_video += 1 + worksheet.write(counter_video,0, database_name) + worksheet.write(counter_video,1, nameStr) + worksheet.write(counter_video,2, str(measures_pred['bpm'])) + worksheet.write(counter_video,3, str(measures_truth['bpm'])) + worksheet.write(counter_video,4, p[0]) + worksheet.write(counter_video,5, str(measures_pred['ibi'])) + worksheet.write(counter_video,6, str(measures_truth['ibi'])) + worksheet.write(counter_video,7, str(measures_pred['sdnn'])) + worksheet.write(counter_video,8, str(measures_truth['sdnn'])) + worksheet.write(counter_video,9, str(measures_pred['rmssd'])) + worksheet.write(counter_video,10, str(measures_truth['rmssd'])) + worksheet.write(counter_video,11, str(measures_pred['pnn50'])) + worksheet.write(counter_video,12, str(measures_truth['pnn50'])) + worksheet.write(counter_video,13, str(measures_pred['lf'])) + worksheet.write(counter_video,14, str(measures_truth['lf'])) + worksheet.write(counter_video,15, str(measures_pred['hf'])) + worksheet.write(counter_video,16, str(measures_truth['hf'])) + try: + worksheet.write(counter_video,17,str(measures_pred['p_total'])) + worksheet.write(counter_video,18,str(measures_truth['p_total'])) + except: + pass + try: + worksheet.write(counter_video,19, str(measures_pred['lf/hf'])) + worksheet.write(counter_video,20, str(measures_truth['lf/hf'])) + worksheet.write(counter_video,21, str(measures_pred['sd1'])) + worksheet.write(counter_video,22, str(measures_truth['sd1'])) + worksheet.write(counter_video,23, str(measures_pred['sd2'])) + worksheet.write(counter_video,24, str(measures_truth['sd2'])) + worksheet.write(counter_video,25, MAE) + except: + pass + + counter_video += 1 + old_database = database_name + +if __name__ == "__main__": + data_dir = 'D:/Databases/3)Testing/' + + GC_path = glob(os.path.join(data_dir, "**/*", '*GC.txt'), recursive=True) + ICA_path = glob(os.path.join(data_dir, "**/*", '*ICA_POH.txt'), recursive=True) + CHROM_path = glob(os.path.join(data_dir, "**/*", '*CHROM.txt'), recursive=True) + + + save_dir = 'D:/Databases/5)Evaluation/Test/' + + workbook = xlsxwriter.Workbook(save_dir + "Result_iPhys" + ".xlsx") + worksheet_GC = workbook.add_worksheet("GC") + write_header(worksheet_GC) + predict_vitals(worksheet_GC, GC_path, save_dir) + print("Ready with this model") + worksheet_ICA = workbook.add_worksheet("ICA") + write_header(worksheet_ICA) + predict_vitals(worksheet_ICA, ICA_path, save_dir) + print("Ready with this model") + worksheet_CHROM = workbook.add_worksheet("CHROM") + write_header(worksheet_CHROM) + predict_vitals(worksheet_CHROM, CHROM_path, save_dir) + print("Ready with this model") + workbook.close() diff --git a/code/final_evaluation.py b/code/final_evaluation.py new file mode 100644 index 0000000..4f974f3 --- /dev/null +++ b/code/final_evaluation.py @@ -0,0 +1,332 @@ +from aifc import Error +import numpy as np +import scipy.io +import xlsxwriter +from model import CAN, CAN_3D, PPTS_CAN, PTS_CAN, TS_CAN, Hybrid_CAN +import h5py +import os +import matplotlib.pyplot as plt +from scipy.signal import butter +from inference_preprocess import preprocess_raw_frames, preprocess_raw_video, detrend +from sklearn.preprocessing import MinMaxScaler +from sklearn import metrics +import scipy.stats as sc +from glob import glob +from scipy import signal + +import heartpy as hp + + +def write_header(worksheet): + header = ['Database', 'Subj/Task', 'HR-pred', 'HR-truth', + 'p','meanNN-pred', 'meanNN-truth', 'sdnn-pred', 'sdnn-truth', + 'rmssd-pred', 'rmssd-truth', 'pNN50-pred', 'pNN50-truth', + 'LF-pred', 'LF-truth', 'HF-pred', 'HF-truth', + 'TP-pred', 'TP-truth', 'LF/HF-pred', 'LF/HF-truth', + 'sd1-pred', 'sd1_truth', 'sd2_pred', 'sd2_truth', 'MAE'] + for index in range(len(header)): + worksheet.write(0,index, header[index]) + +def prepare_3D_CAN(dXsub): + frame_depth = 10 + num_window = int(dXsub.shape[0]) - frame_depth + 1 + tempX = np.array([dXsub[f:f + frame_depth, :, :, :] # (491, 10, 36, 36 ,6) (169, 10, 36, 36, 6) + for f in range(num_window)]) + tempX = np.swapaxes(tempX, 1, 3) # (169, 36, 36, 10, 6) + tempX = np.swapaxes(tempX, 1, 2) # (169, 36, 36, 10, 6) + return tempX + +def prepare_Hybrid_CAN(dXsub): + frame_depth = 10 + num_window = int(dXsub.shape[0]) - frame_depth + 1 + tempX = np.array([dXsub[f:f + frame_depth, :, :, :] # (169, 10, 36, 36, 6) + for f in range(num_window)]) + tempX = np.swapaxes(tempX, 1, 3) # (169, 36, 36, 10, 6) + tempX = np.swapaxes(tempX, 1, 2) # (169, 36, 36, 10, 6) + motion_data = tempX[:, :, :, :, :3] + apperance_data = np.average(tempX[:, :, :, :, -3:], axis=-2) + return motion_data, apperance_data + +def predict_vitals(worksheet, test_name, model_name, video_path, path_results): + mms = MinMaxScaler() + img_rows = 36 + img_cols = 36 + frame_depth = 10 + batch_size = 100 + + #### initialize Model ###### + try: + model_checkpoint = os.path.join(path_results, test_name, "cv_0_epoch24_model.hdf5") + except: + model_checkpoint = os.path.join(path_results, test_name, "cv_0_epoch23_model.hdf5") + batch_size = batch_size + + if model_name == "TS_CAN": + model = TS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3)) + elif model_name == "3D_CAN": + model = CAN_3D(frame_depth, 32, 64, (img_rows, img_cols, frame_depth, 3)) + elif model_name == "CAN": + model = CAN(32, 64, (img_rows, img_cols, 3)) + elif model_name == "Hybrid_CAN": + model = Hybrid_CAN(frame_depth, 32, 64, (img_rows, img_cols, frame_depth, 3), + (img_rows, img_cols, 3)) + elif model_name == "PTS_CAN": + model = PTS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3)) + elif model_name == "PPTS_CAN": + model = PPTS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3), parameter=['bpm', 'sdnn']) + else: + raise NotImplementedError + + model.load_weights(model_checkpoint) + ###### load video Data ####### + counter_video = 1 + old_database = "COH" + for sample_data_path in video_path: + print("path: ",sample_data_path) + if sample_data_path[-4:] == ".avi": + dXsub, fs = preprocess_raw_video(sample_data_path, dim=36) + else: + dXsub, fs = preprocess_raw_frames(sample_data_path, dim=36) + print('dXsub shape', dXsub.shape, "fs: ", fs) + + if model_name == "PPTS_CAN": + dXsub_len = (dXsub.shape[0] // (frame_depth*10)) * (frame_depth*10) + dXsub = dXsub[:dXsub_len, :, :, :] + + else: + dXsub_len = (dXsub.shape[0] // frame_depth) * frame_depth + dXsub = dXsub[:dXsub_len, :, :, :] + + if model_name == "3D_CAN": + dXsub = prepare_3D_CAN(dXsub) + dXsub_len = (dXsub.shape[0] // (frame_depth)) * (frame_depth) + dXsub = dXsub[:dXsub_len, :, :, :,:] + yptest = model.predict((dXsub[:, :, :,: , :3], dXsub[:, :, :, : , -3:]), verbose=1) + elif model_name == "Hybrid_CAN": + dXsub1, dXsub2 = prepare_Hybrid_CAN(dXsub) + dXsub_len1 = (dXsub1.shape[0] // (frame_depth*10)) * (frame_depth*10) + dXsub1 = dXsub1[:dXsub_len1, :, :, :, :] + dXsub_len2 = (dXsub2.shape[0] // (frame_depth*10)) * (frame_depth*10) + dXsub2 = dXsub2[:dXsub_len2, :, :, :] + yptest = model.predict((dXsub1, dXsub2), verbose=1) + else: + yptest = model((dXsub[:, :, :, :3], dXsub[:, :, :, -3:]), training=False) + + if model_name == "3D_CAN" or model_name == "Hybrid_CAN": + pulse_pred = yptest[:,0] + elif model_name != "PTS_CAN" and model_name != "PPTS_CAN": + pulse_pred = yptest + + else: + pulse_pred = yptest[0] + + pulse_pred = detrend(np.cumsum(pulse_pred), 100) + [b_pulse_pred, a_pulse_pred] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') + pulse_pred = scipy.signal.filtfilt(b_pulse_pred, a_pulse_pred, np.double(pulse_pred)) + pulse_pred = np.array(mms.fit_transform(pulse_pred.reshape(-1,1))).flatten() + + ##### ground truth data resampled ####### + if(str(sample_data_path).find("COHFACE") > 0): + database_name = "COH" + truth_path = sample_data_path.replace(".avi", "_dataFile.hdf5") # akutell für COHACE... + elif(str(sample_data_path).find("UBFC-PHYS") > 0): + database_name = "UB-Ph" + truth_path = sample_data_path.replace("vid_", "").replace(".avi","_dataFile.hdf5") + elif(str(sample_data_path).find("UBFC") > 0): + database_name = "UBFC" + truth_path = sample_data_path.replace("vid.avi", "dataFile.hdf5") + elif(str(sample_data_path).find("BP4D") > 0): + database_name = "BP4D" + truth_path = sample_data_path + "/BP_mmHg.txt" + else: + return print("Error in finding the ground truth signal...") + if database_name != "BP4D": + gound_truth_file = h5py.File(truth_path, "r") + pulse_truth = gound_truth_file["pulse"] ### range ground truth from 0 to 1 + pulse_truth = pulse_truth[0:dXsub_len] + else: + data = open(truth_path, 'r').read() + data = str(data).split("\n") + pulse_truth = np.array(list(map(float, data[0:-1]))) + mms = MinMaxScaler() + mean = pulse_truth.mean() + std = np.std(pulse_truth) + upper_limit = mean + std*3 + lower_limit = mean - std*3 + for x in range(0, len(pulse_truth)): + if pulse_truth[x] > upper_limit: + pulse_truth[x] = upper_limit + elif pulse_truth[x] < lower_limit: + pulse_truth[x] = lower_limit + + pulse_truth = np.array(signal.resample(pulse_truth, len(pulse_pred))) + pulse_truth = np.array(mms.fit_transform(pulse_truth.reshape(-1,1))).flatten() # normalization + pulse_truth = detrend(np.cumsum(pulse_truth), 100) + [b_pulse_tr, a_pulse_tr] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') + pulse_truth = scipy.signal.filtfilt(b_pulse_tr, a_pulse_tr, np.double(pulse_truth)) + pulse_truth = np.array(mms.fit_transform(pulse_truth.reshape(-1,1))).flatten() + ### same size ####### + if len(pulse_pred) > len(pulse_truth): + pulse_pred = pulse_pred[:len(pulse_truth)] + elif len(pulse_pred) < len(pulse_truth): + pulse_truth = pulse_truth[:len(pulse_pred)] + ########### Peaks ########### + working_data_pred, measures_pred = hp.process(pulse_pred, fs, calc_freq=True) + working_data_truth, measures_truth = hp.process(pulse_truth, fs, calc_freq=True) + peaks_pred = working_data_pred['peaklist'] + peaks_truth = working_data_truth['peaklist'] + + ######## name files ############# + if(str(sample_data_path).find("COHFACE") > 0): + nmr = str(sample_data_path).find("COHFACE") + nameStr = str(sample_data_path)[nmr + 7:].replace("\\", "-").replace("-data.avi", "") + elif(str(sample_data_path).find("UBFC-PHYS") > 0): + nmr = str(sample_data_path).find("UBFC-PHYS") + nameStr = str(sample_data_path)[nmr + 12:].replace("\\", "-").replace("vid_", "").replace(".avi", "") + elif(str(sample_data_path).find("UBFC") > 0): + nmr = str(sample_data_path).find("UBFC") + nameStr = str(sample_data_path)[nmr + 5:].replace("\\", "-").replace("vid.avi", "") + elif(str(sample_data_path).find("BP4D") > 0): + nmr = str(sample_data_path).find("BP4D") + nameStr = str(sample_data_path)[nmr + 5:].replace("\\", "-") + else: + raise ValueError + ########## Plot ################## + peaks_pred_new = [] + for peak in peaks_pred: + if (peak > 400 and peak <700): + peaks_pred_new.append(peak-400) + peaks_truth_new = [] + for peak in peaks_truth: + if (peak > 400 and peak <700): + peaks_truth_new.append(peak-400) + plt.figure() #subplot(211) + plt.plot(pulse_pred[400:700], "#E6001A", label='rPPG signal') + plt.plot(peaks_truth_new, pulse_truth[400:700][peaks_truth_new], "x", color="#005AA9") + plt.plot(peaks_pred_new, pulse_pred[400:700][peaks_pred_new], "x", color ='#E6001A') + plt.title('rPPG signal with ground truth') + plt.ylabel("normalized Signal [a.u.]") + plt.xlabel("time (samples)") + plt.plot(pulse_truth[400:700], '#005AA9', linewidth=0.9, label='ground truth') + plt.legend() + plt.savefig(database_name+ nameStr + "_both.svg", format="svg") + + plt.figure() + plt.subplot(211) + plt.plot(pulse_truth[400:700],"#004E8A", label='Ground truth') + plt.plot(peaks_truth_new, pulse_truth[400:700][peaks_truth_new], "x", color="#004E8A") + plt.ylabel("normalized Signal [a.u.]") + plt.title('Ground truth') + plt.subplot(212) + plt.plot(pulse_pred[400:700], "#004E8A",label='Prediction') + plt.plot(peaks_pred_new, pulse_pred[400:700][peaks_pred_new],"x", color="#004E8A") + plt.title("Predicted rPPG") + plt.ylabel("normalized Signal [a.u.]") + plt.xlabel("time (samples)") + plt.legend() + plt.savefig(database_name+ nameStr +".svg", format="svg") + + ########### IBI ############# + #ibi_truth = working_data_truth['RR_list_cor'] + #print(ibi_truth) + #ibi_pred = working_data_pred['RR_list_cor'] + #print(ibi_pred) + ######### HRV featurs ############## + #print("HRV Truth: ",measures_truth) + #print("HRV Pred: ", measures_pred) + ######### Metrics ############## + # MSE: + MAE = metrics.mean_absolute_error(pulse_truth, pulse_pred) + # RMSE: + RMSE = metrics.mean_squared_error(pulse_truth, pulse_pred, squared=False) + # Pearson correlation: + p = sc.pearsonr(pulse_truth, pulse_pred) + + ####### Logging ############# + if database_name != old_database: + counter_video += 1 + worksheet.write(counter_video,0, database_name) + worksheet.write(counter_video,1, nameStr) + worksheet.write(counter_video,2, measures_pred['bpm']) + worksheet.write(counter_video,3, measures_truth['bpm']) + worksheet.write(counter_video,4, p[0]) + worksheet.write(counter_video,5, measures_pred['ibi']) + worksheet.write(counter_video,6, measures_truth['ibi']) + worksheet.write(counter_video,7, measures_pred['sdnn']) + worksheet.write(counter_video,8, measures_truth['sdnn']) + worksheet.write(counter_video,9, measures_pred['rmssd']) + worksheet.write(counter_video,10, measures_truth['rmssd']) + worksheet.write(counter_video,11, measures_pred['pnn50']) + worksheet.write(counter_video,12, measures_truth['pnn50']) + worksheet.write(counter_video,13, measures_pred['lf_perc']) + worksheet.write(counter_video,14, measures_truth['lf_perc']) + worksheet.write(counter_video,15, measures_pred['hf_perc']) + worksheet.write(counter_video,16, measures_truth['hf_perc']) + worksheet.write(counter_video,17, measures_pred['p_total']) + worksheet.write(counter_video,18, measures_truth['p_total']) + worksheet.write(counter_video,19, measures_pred['lf/hf']) + worksheet.write(counter_video,20, measures_truth['lf/hf']) + worksheet.write(counter_video,21, measures_pred['sd1']) + worksheet.write(counter_video,22, measures_truth['sd1']) + worksheet.write(counter_video,23, measures_pred['sd2']) + worksheet.write(counter_video,24, measures_truth['sd2']) + worksheet.write(counter_video,25, MAE) + + counter_video += 1 + old_database = database_name + +if __name__ == "__main__": + #path_results = "D:/Databases/4)Results/Version5" #finalVersions" + path_results = "/home/quehl/Results/actualResults" + #data_dir = "C:/Users/sarah/Desktop"#\F001" + #data_dir = "D:/Databases/3)Testing/" + data_dir = '/mnt/share/StudiShare/sarah/Databases/Testing' + modelDir_names = glob(path_results +"/*") + testModel_names = [] + for dir in modelDir_names: + split = dir.split("\\") + testModel_names.append(split[-1]) + + video_path = glob(os.path.join(data_dir, "**/*", '*.avi'), recursive=True) + #video_path = glob(os.path.join(data_dir, "COHFACE/**/*", '*.avi'), recursive=True) + #video_path += glob(os.path.join(data_dir, "UBFC/**/*", '*.avi'), recursive=True) + ### BP4D #### + new_dirs = glob(os.path.join(data_dir, "BP4D/**/*")) + video_path = video_path + new_dirs + + #testModel_names=['PTS_CAN_Gauss2']#, 'PPTS_CAN_negPea_TE_sdnn_pnn50_lfhf'] + #save_dir = "D:/Databases/5)Evaluation/P_Evaluation_Mix2"#finalEvaluation" + save_dir = '/home/quehl/finalEvaluation/' + print("Models: ", testModel_names) + for test_name in testModel_names: + print("Current Modelname: ", test_name) + if str(test_name).find("3D_CAN") >=0: + model_name = "3D_CAN" + elif str(test_name).find("Hybrid_CAN") >= 0: + model_name = "Hybrid_CAN" + elif str(test_name).find("TS_CAN") >= 0: + model_name = "TS_CAN" + elif str(test_name).find("PPTS") >= 0: + model_name = "PPTS_CAN" + elif str(test_name).find("PTS") >= 0: + model_name = "PTS_CAN" + else: + if str(test_name).find("CAN") >= 0: + model_name = "CAN" + else: + raise Error("Model not found...") + + # neuer Ordner für Tests + os.chdir(save_dir) + try: + os.makedirs(str(test_name)) + except: + print("Directory exists...") + save_path = os.path.join(save_dir, str(test_name)) + os.chdir(save_path) + workbook = xlsxwriter.Workbook(test_name + ".xlsx") + worksheet = workbook.add_worksheet("Results") + write_header(worksheet) + predict_vitals(worksheet, test_name, model_name, video_path, path_results) + print("Ready with this model") + workbook.close() diff --git a/code/inference_preprocess.py b/code/inference_preprocess.py index 517049b..47c1dba 100644 --- a/code/inference_preprocess.py +++ b/code/inference_preprocess.py @@ -1,3 +1,5 @@ +from glob import glob +from importlib import import_module import numpy as np import cv2 from skimage.util import img_as_float @@ -6,6 +8,7 @@ import time import scipy.io from scipy.sparse import spdiags +from tensorflow.python.keras import backend as K def preprocess_raw_video(videoFilePath, dim=36): @@ -13,46 +16,112 @@ def preprocess_raw_video(videoFilePath, dim=36): # set up t = [] i = 0 - vidObj = cv2.VideoCapture(videoFilePath); + vidObj = cv2.VideoCapture(videoFilePath) + totalFrames = int(vidObj.get(cv2.CAP_PROP_FRAME_COUNT)) # get total frame size + fps = vidObj.get(cv2.CAP_PROP_FPS) + #print("fps: ", fps) + + #totalFrames = 2101 Xsub = np.zeros((totalFrames, dim, dim, 3), dtype = np.float32) height = vidObj.get(cv2.CAP_PROP_FRAME_HEIGHT) width = vidObj.get(cv2.CAP_PROP_FRAME_WIDTH) success, img = vidObj.read() dims = img.shape - print("Orignal Height", height) - print("Original width", width) + #print("Orignal Height", height) + #print("Original width", width) + ######################################################################### # Crop each frame size into dim x dim while success: t.append(vidObj.get(cv2.CAP_PROP_POS_MSEC))# current timestamp in milisecond - vidLxL = cv2.resize(img_as_float(img[:, int(width/2)-int(height/2 + 1):int(height/2)+int(width/2), :]), (dim, dim), interpolation = cv2.INTER_AREA) - vidLxL = cv2.rotate(vidLxL, cv2.ROTATE_90_CLOCKWISE) # rotate 90 degree + vidLxL = cv2.resize(img_as_float(img), (dim, dim), interpolation = cv2.INTER_AREA) #img[:, int(width/2)-int(height/2 + 1): int(height/2)+int(width/2), :]) + #vidLxL = cv2.rotate(vidLxL, cv2.ROTATE_90_CLOCKWISE) # rotate 90 degree vidLxL = cv2.cvtColor(vidLxL.astype('float32'), cv2.COLOR_BGR2RGB) vidLxL[vidLxL > 1] = 1 vidLxL[vidLxL < (1/255)] = 1/255 Xsub[i, :, :, :] = vidLxL success, img = vidObj.read() # read the next one i = i + 1 - plt.imshow(Xsub[0]) - plt.title('Sample Preprocessed Frame') - plt.show() + if i >= totalFrames: + break + # plt.imshow(Xsub[0]) + # plt.title('Sample Preprocessed Frame') + # plt.show() ######################################################################### # Normalized Frames in the motion branch normalized_len = len(t) - 1 + + #print("normalized Len") + #print(normalized_len) + dXsub = np.zeros((normalized_len, dim, dim, 3), dtype = np.float32) + for j in range(normalized_len - 1): + dXsub[j, :, :, :] = (Xsub[j+1, :, :, :] - Xsub[j, :, :, :]) / (Xsub[j+1, :, :, :] + Xsub[j, :, :, :]) + dXsub = dXsub / np.std(dXsub) + # plt.imshow(dXsub[0]) + # plt.title('Sample Preprocessed Frame') + # plt.show() + + ######################################################################### + # Normalize raw frames in the apperance branch + Xsub = Xsub - np.mean(Xsub) + Xsub = Xsub / np.std(Xsub) + Xsub = Xsub[:dXsub.shape[0], :, :, :] # -1 + ######################################################################### + # Plot an example of data after preprocess + dXsub = np.concatenate((dXsub, Xsub), axis = 3) + return dXsub, fps + +def preprocess_raw_frames(framePath, fps=25, dim=36): + ####### collect frames ######### + frames = glob(framePath + "/*.jpg") + ######################################################################### + # set up + i = 0 + frames = sorted(frames) + totalFrames = int(len(frames)) # get total frame size + print(totalFrames) + Xsub = np.zeros((totalFrames, dim, dim, 3), dtype = np.float32) + + ######################################################################### + # Crop each frame size into dim x dim + while i <= totalFrames: + img = cv2.imread(frames[i]) + print(img) + #t.append(vidObj.get(cv2.CAP_PROP_POS_MSEC))# current timestamp in milisecond + vidLxL = cv2.resize(img_as_float(img), (dim, dim), interpolation = cv2.INTER_AREA) + #vidLxL = cv2.rotate(vidLxL, cv2.ROTATE_90_CLOCKWISE) # rotate 90 degree + vidLxL = cv2.cvtColor(vidLxL.astype('float32'), cv2.COLOR_BGR2RGB) + vidLxL[vidLxL > 1] = 1 + vidLxL[vidLxL < (1/255)] = 1/255 + Xsub[i, :, :, :] = vidLxL + i = i + 1 + #plt.imshow(Xsub[0]) + #plt.title('Sample Preprocessed Frame') + #plt.show() + ######################################################################### + # Normalized Frames in the motion branch + normalized_len = totalFrames - 1 + + #print("normalized Len") + #print(normalized_len) dXsub = np.zeros((normalized_len, dim, dim, 3), dtype = np.float32) for j in range(normalized_len - 1): dXsub[j, :, :, :] = (Xsub[j+1, :, :, :] - Xsub[j, :, :, :]) / (Xsub[j+1, :, :, :] + Xsub[j, :, :, :]) dXsub = dXsub / np.std(dXsub) + # plt.imshow(dXsub[0]) + # plt.title('Sample Preprocessed Frame') + # plt.show() + ######################################################################### # Normalize raw frames in the apperance branch Xsub = Xsub - np.mean(Xsub) Xsub = Xsub / np.std(Xsub) - Xsub = Xsub[:totalFrames-1, :, :, :] + Xsub = Xsub[:dXsub.shape[0], :, :, :] # -1 ######################################################################### # Plot an example of data after preprocess - dXsub = np.concatenate((dXsub, Xsub), axis = 3); - return dXsub + dXsub = np.concatenate((dXsub, Xsub), axis = 3) + return dXsub, fps def detrend(signal, Lambda): """detrend(signal, Lambda) -> filtered_signal diff --git a/code/layer_output.py b/code/layer_output.py new file mode 100644 index 0000000..1dcc1f2 --- /dev/null +++ b/code/layer_output.py @@ -0,0 +1,99 @@ +import os +import scipy +from model import CAN_3D, PPTS_CAN, PTS_CAN, TS_CAN +from data_generator import DataGenerator +import matplotlib.pyplot as plt +import tensorflow as tf +import heartpy as hp +import numpy as np +from scipy.signal import butter +from inference_preprocess import preprocess_raw_video, detrend +from sklearn.preprocessing import MinMaxScaler + + + +def gaussian_loss(y_true, y_pred): + y_pred = tf.reshape(y_pred, (-1,)) + return -tf.reduce_sum(tf.abs(y_true*y_pred)) + +path_of_video_tr = ["D:/Databases/3)Testing/COHFACE/25/1/data_dataFile.hdf5"] + +model = PTS_CAN(10, 32, 64, (36,36,3), + dropout_rate1=0.25, dropout_rate2=0.5, nb_dense=128) #, parameter=['bpm', 'sdnn', 'pnn50', 'lf_hf'] +training_generator = DataGenerator(path_of_video_tr, 2100, (36, 36), + batch_size=1, frame_depth=10, + temporal="PTS_CAN", respiration=False, database_name="COHFACE", + time_error_loss=True) # truth_parameter=['bpm', 'sdnn', 'pnn50', 'lf_hf'] + +inp = model.input # input placeholder +#model.summary() +model_checkpoint = os.path.join("D:/Databases/4)Results/Version5/TS_CAN/cv_0_epoch24_model.hdf5")#PPTS_CAN_bpm_sdnn/cv_0_epoch24_model.hdf5") +model.load_weights(model_checkpoint) +#outputs = [layer.output for layer in model.layers] # all layer outputs +#functors = [K.function([inp, K.learning_phase()], [out]) for out in outputs] # evaluation functions + +# Testing +test = training_generator.data_generation(path_of_video_tr) +output = model(test) + +#### test new output ###### +# pred = tf.reshape(output[2], (-1)) +# truth = tf.reshape(test[1][2], (-1)) +# AE = tf.abs(pred - truth)/truth +# loss = tf.reduce_mean(AE) +# diff = pred - truth +#print(output[1]) +fs = 20 +mms = MinMaxScaler() +pulse_pred = np.array(tf.reshape(output[0], (-1,))) +pulse_pred2 = detrend(np.cumsum(pulse_pred), 100) +[b_pulse_pred, a_pulse_pred] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') +pulse_pred2 = scipy.signal.filtfilt(b_pulse_pred, a_pulse_pred, np.double(pulse_pred2)) +pulse_pred2 = np.array(mms.fit_transform(pulse_pred2.reshape(-1,1))).flatten() + +pulse_true = test[1][0] +working_data_pred, measures_pred= hp.process(pulse_pred, fs, calc_freq=True) +working_data_true, measures_true = hp.process(pulse_true, fs, calc_freq=True) +peaks_pred = working_data_pred['peaklist'] +peaks_true = working_data_true['peaklist'] +bin_arr = np.array(tf.reshape(output[1], (-1,))) +x_loss = np.where(bin_arr == 1)[0] + +mult = np.array(tf.reshape(output[1], (-1,))) * test[1][1] +y_loss =np.delete(mult, np.where(bin_arr ==0)) +plt.figure() +plt.subplot(211) +plt.title('TE loss function example') +plt.plot(pulse_pred, label='rPPG$_{out}$', linewidth=1, color="#B90F22") +plt.plot(peaks_pred, pulse_pred[peaks_pred], "x", color="#B90F22") +plt.plot(pulse_true, label='ground truth',linewidth=1, color="#004E8A") +plt.plot(peaks_true, pulse_true[peaks_true], "x", color="#004E8A") +plt.ylabel("rPPG [a.u.]") +plt.legend(loc="upper right") +plt.subplot(212) +plt.title('Corrensponding Loss') +plt.vlines(np.where(output[1] == 1), ymin=0, ymax=1 , label='binary$_{out}$', color="#B90F22")# +plt.hlines(0, 0, 505, color= "#B90F22") +plt.plot(test[1][1], label='ground truth',linewidth=1, color="#004E8A") +plt.plot(x_loss, y_loss , "x", color="k", label="Loss") +plt.ylabel("seconds") +plt.xlabel("time (samples)") +plt.legend(loc="upper right") +plt.show() + + +y_true = test[1][1] +y_pred = output[1] +loss = gaussian_loss(y_true, y_pred) +mult = y_true * tf.reshape(y_pred, (-1,)) +print(np.sum(mult)) +plt.plot(mult) +plt.title("Multiplication of y_true and y_pred") +plt.ylabel("[a.u.]") +plt.xlabel("time (samples)") +plt.show() + +#layer_outs = [func([test, 1.]) for func in functors] +#print(layer_outs) + + diff --git a/code/losses.py b/code/losses.py new file mode 100644 index 0000000..5a1b432 --- /dev/null +++ b/code/losses.py @@ -0,0 +1,116 @@ +###### LOSS FUNCTIONS ############### +# Defines different Loss Functions. +# Currently implemented: +# - Negative Pearson Coefficient + +from tkinter.tix import Y_REGION +import tensorflow as tf + +import tensorflow.keras.backend as K + +# Negative Pearson Coefficient +# x: truth rPPG y: predicted rPPG +def negPearsonLoss(x,y): + mean_x = tf.reduce_mean(x) + mean_y = tf.reduce_mean(y) + + x_1 = x - mean_x + y_1 = y - mean_y + + s_xy = tf.reduce_sum(tf.multiply(x_1,y_1)) + s_x = tf.reduce_sum(x_1**2) + s_y = tf.reduce_sum(y_1**2) + sx_sy = tf.sqrt(tf.multiply(s_x,s_y)) + + p = tf.divide(s_xy, sx_sy) + + negPearson_coeff = 1. - p + + return negPearson_coeff + + +def gaussian_loss(y_true, y_pred): + y_pred = tf.reshape(y_pred, (-1,)) + y_true = tf.reshape(y_true, (-1,)) + return -tf.reduce_sum(y_true*y_pred) + +def time_error_loss(y_true, y_pred): + y_pred = tf.reshape(y_pred, (-1,)) + y_true = tf.reshape(y_true, (-1,)) + return tf.reduce_sum(y_true*y_pred) + +def MAPE_parameter_loss(y_true, y_pred): + y_true = tf.reshape(y_true, (-1,)) + y_pred = tf.reshape(y_pred, (-1,)) + y_true = tf.maximum(y_true, 1e-6) + AE = tf.abs(y_true-y_pred)/y_true + return tf.reduce_mean(AE) + +#not working ---> PTS-CAN... +def negPearsonLoss_onlyPeaks(y_true, y_pred): + peaks_true = get_peaks(y_true) + peaks_pred = get_peaks(y_pred) + peaks_pred = peaks_pred[0:10] + peaks_true = peaks_true[0:10] + #peaks_true, peaks_pred = filt_peaks(peaks_true, peaks_pred) + peaks_true = tf.cast(peaks_true, tf.float32) + peaks_pred = tf.cast(peaks_pred, tf.float32) + + negPeaLoss = negPearsonLoss(peaks_true, peaks_pred) + return negPeaLoss + + +def get_peaks(y): + # y: (N,) + data_reshaped = tf.reshape(y, (1, -1, 1)) # (1, N, 1) + max_pooled_in_tensor = tf.nn.max_pool(data_reshaped, (20,), 1,'SAME') + + #maxima = tf.stop_gradient(tf.equal(data_reshaped,max_pooled_in_tensor)) # (1, N, 1) + #maxima = tf.cast(maxima, tf.float32) + #maxima = tf.squeeze(maxima) # (N,1) + #peaks = tf.where(maxima, name="Where") # now only the Peak Indices (A, 3) + #tf.no_gradient("Where") + #peaks = tf.reshape(peaks, (-1,)) # (A,1) + + return max_pooled_in_tensor + +# x: true y: prediction +# input: peaks of truth and prediction as tensor... +@tf.function +def filt_peaks(x,y): + def true_fn(): + return min + def false_fn(): + return tf.cast(-1, tf.int64) + max_offset = 10 + mask = tf.cast(tf.zeros(tf.size(x)),tf.bool) # tensor with size of x (truth data) + # check which peaks of truth are recognized in pred + min = 0 + min = tf.cast(min, tf.int64) + + def fn(item): + diff = tf.abs(x - item) # diff of truth data and item + min = tf.reduce_min(diff) # minimum of diff + min = tf.cond(tf.less(min, max_offset), true_fn, false_fn) + temp_mask = tf.equal(min, diff) + mask = tf.logical_or(mask, temp_mask) + return mask + mask = tf.map_fn(fn=lambda item: fn(item), elems=y) + # for item in y: # items of predicion + # diff = tf.abs(x - item) # diff of truth data and item + # min = tf.reduce_min(diff) # minimum of diff + # min = tf.cond(tf.less(min, max_offset), true_fn, false_fn) + # temp_mask = tf.equal(min, diff) + # mask = tf.logical_or(mask, temp_mask) + + x = tf.boolean_mask(x, mask) + # check if outliners are in pred + mask = tf.cast(tf.zeros(tf.size(y)), tf.bool) + for item in x: + diff = tf.abs(y - item) # diff of truth data and item + min = tf.reduce_min(diff) # minimum of diff + min = tf.cond(tf.less(min, max_offset), true_fn, false_fn) + temp_mask = tf.equal(min, diff) + mask = tf.logical_or(mask, temp_mask) + y = tf.boolean_mask(y,mask) + return x, y \ No newline at end of file diff --git a/code/model.py b/code/model.py index c41fe87..603b881 100755 --- a/code/model.py +++ b/code/model.py @@ -1,16 +1,16 @@ ''' Models for Multi-Task Temporal Shift Attention Networks for On-Device Contactless Vitals Measurement Author: Xin Liu +further developed: Sarah Quehl ''' - +from re import L, T +from numpy import float32 import tensorflow as tf -from tensorflow import keras from tensorflow.python.keras import backend as K from tensorflow.python.keras.layers import Conv2D, Conv3D, Input, AveragePooling2D, \ multiply, Dense, Dropout, Flatten, AveragePooling3D from tensorflow.python.keras.models import Model - class Attention_mask(tf.keras.layers.Layer): def call(self, x): xsum = K.sum(x, axis=1, keepdims=True) @@ -22,7 +22,6 @@ def get_config(self): config = super(Attention_mask, self).get_config() return config - class TSM(tf.keras.layers.Layer): def call(self, x, n_frame, fold_div=3): nt, h, w, c = x.shape @@ -60,10 +59,118 @@ def TSM_Cov2D(x, n_frame, nb_filters=128, kernel_size=(3, 3), activation='tanh', x = Conv2D(nb_filters, kernel_size, padding=padding, activation=activation)(x) return x +# own layer: +class ownLayer_binaryPeak(tf.keras.layers.Layer): + def call(self, x): + out = self.get_peaks(x) + + return out + + def get_peaks(self, y): + # y: (N,1) + data_reshaped = tf.reshape(y, (1, -1, 1)) # (1, N, 1) + max_pooled_in_tensor = tf.nn.max_pool(data_reshaped, (20,), 1,'SAME') + maxima = tf.equal(data_reshaped, max_pooled_in_tensor) # (1, N, 1) + maxima = tf.cast(maxima, tf.float32) + maxima = tf.reshape(maxima, (-1,1)) -# %% + return maxima + def get_config(self): + config = super(ownLayer_binaryPeak, self).get_config() + return config + +class ownLayer_parameter(tf.keras.layers.Layer): + def call(self, x, parameter): + rr = self.get_rr(x) + + f_bpm = lambda: self.get_HR(tf.cast(rr, dtype=float32)) + f_sdnn = lambda: self.get_sdnn(tf.cast(rr, dtype=float32)) + f_pnn50 = lambda: self.get_pNN50(tf.cast(rr, dtype=float32)) + f_lfhf = lambda: self.get_lf_hf(tf.cast(rr, dtype=float32)) + + result = [] + for item in parameter: + result_part = tf.case([(tf.equal(item,'bpm'), f_bpm), (tf.equal(item,'sdnn'), f_sdnn), (tf.equal(item,'pnn50'), f_pnn50), (tf.equal(item,'lf_hf'), f_lfhf)], default=f_bpm) + result.append(result_part) + + result = tf.convert_to_tensor(result) + result = tf.reshape(result, (-1,1)) + return result + + def get_rr(self, y): + # y: (N,1) + fs = 50 + fs = tf.cond(tf.less(tf.shape(tf.reshape(y, (-1,))),tf.convert_to_tensor(1300)), lambda: tf.cast(50, dtype=tf.int64), lambda: tf.cast(40, dtype=tf.int64)) + + indices = tf.where(tf.equal(tf.reshape(y, (-1,)),1)) + peak_locations = tf.squeeze(indices) + + def tf_diff_axis_0(a): + return a[1:]-a[:-1] + ibi_arr = tf_diff_axis_0(peak_locations)*fs + + mask = tf.logical_and(tf.greater_equal(ibi_arr,333),tf.less_equal(ibi_arr, 1500)) + mask.set_shape([None]) + + rr_arr = tf.boolean_mask(ibi_arr, mask) + + return rr_arr + + def get_HR(self, rr): + rr_mean = tf.reduce_mean(rr) + HR = 60000/rr_mean + return HR + + def get_sdnn(self, rr): + return tf.math.reduce_std(rr) + + def get_pNN50(self, rr): + def tf_diff_axis_0(a): + return a[1:]-a[:-1] + rr_diff = tf_diff_axis_0(rr) + + size = tf.cast(tf.reduce_sum(tf.ones(tf.size(rr))), dtype=tf.float32) + + mask = (tf.greater(tf.abs(rr_diff),50)) + + mask = tf.cast(mask, dtype=tf.int32) + nn50 = tf.cast(tf.math.reduce_sum(mask), dtype=tf.float32) + pNN50 = nn50/size + return pNN50 + + def get_lf_hf(self,rr): + data = tf.reshape(rr, (-1,)) + f_true = lambda: data + f_false = lambda: tf.convert_to_tensor([1], dtype=tf.float32) + + def custom_func(data): + frq = tf.cast(tf.abs(tf.signal.rfft(data)), tf.float32)/tf.cast(tf.size(data),tf.float32) + return tf.multiply(tf.pow(frq,2),tf.math.sqrt(tf.cast(2, tf.float32))) + data = tf.case([(tf.greater(tf.size(rr), 2),f_true), (tf.less(tf.size(rr), 2),f_false)]) + + frq = tf.keras.layers.Lambda(custom_func)(data) + + #frq = tf.cast(tf.abs(tf.signal.rfft(data)), tf.float32)/tf.cast(tf.size(data),tf.float32) + dt = tf.math.reduce_mean(data) / 1000 # in sec + t = tf.cast(tf.range(0, tf.size(frq)), tf.float32) + t = tf.cast(t, tf.float32)/(tf.cast(dt, tf.float32)*tf.cast(tf.size(frq)*2, tf.float32)) + + mask_lf = tf.cast(tf.logical_and(tf.greater_equal(t, 0.04), tf.less(t, 0.15)), tf.float32) + lf = tf.maximum(tf.reduce_sum(frq*mask_lf), 0.000001) + mask_hf = tf.cast(tf.logical_and(tf.greater_equal(t, 0.15), tf.less(t, 0.4)), tf.float32) + hf = tf.maximum(tf.reduce_sum(frq*mask_hf), 0.000001) + + lf_hf = lf/hf + + return lf_hf + + def get_config(self): + config = super(ownLayer_parameter, self).get_config() + return config +# %% +# DEEPPHYS???? def CAN(nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3), dropout_rate1=0.25, dropout_rate2=0.5, pool_size=(2, 2), nb_dense=128): diff_input = Input(shape=input_shape) @@ -105,15 +212,14 @@ def CAN(nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3), dropout_rate1 model = Model(inputs=[diff_input, rawf_input], outputs=out) return model - -# %% MT_CAN -def MT_CAN(nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3), dropout_rate1=0.25, dropout_rate2=0.5, +# %% TS_CAN --> Paper +def TS_CAN(n_frame, nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3), dropout_rate1=0.25, dropout_rate2=0.5, pool_size=(2, 2), nb_dense=128): diff_input = Input(shape=input_shape) rawf_input = Input(shape=input_shape) - d1 = Conv2D(nb_filters1, kernel_size, padding='same', activation='tanh')(diff_input) - d2 = Conv2D(nb_filters1, kernel_size, activation='tanh')(d1) + d1 = TSM_Cov2D(diff_input, n_frame, nb_filters1, kernel_size, padding='same', activation='tanh') + d2 = TSM_Cov2D(d1, n_frame, nb_filters1, kernel_size, padding='valid', activation='tanh') r1 = Conv2D(nb_filters1, kernel_size, padding='same', activation='tanh')(rawf_input) r2 = Conv2D(nb_filters1, kernel_size, activation='tanh')(r1) @@ -128,8 +234,8 @@ def MT_CAN(nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3), dropout_ra r3 = AveragePooling2D(pool_size)(r2) r4 = Dropout(dropout_rate1)(r3) - d5 = Conv2D(nb_filters2, kernel_size, padding='same', activation='tanh')(d4) - d6 = Conv2D(nb_filters2, kernel_size, activation='tanh')(d5) + d5 = TSM_Cov2D(d4, n_frame, nb_filters2, kernel_size, padding='same', activation='tanh') + d6 = TSM_Cov2D(d5, n_frame, nb_filters2, kernel_size, padding='valid', activation='tanh') r5 = Conv2D(nb_filters2, kernel_size, padding='same', activation='tanh')(r4) r6 = Conv2D(nb_filters2, kernel_size, activation='tanh')(r5) @@ -142,21 +248,14 @@ def MT_CAN(nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3), dropout_ra d8 = Dropout(dropout_rate1)(d7) d9 = Flatten()(d8) - d10_y = Dense(nb_dense, activation='tanh')(d9) - d11_y = Dropout(dropout_rate2)(d10_y) - out_y = Dense(1, name='output_1')(d11_y) - - d10_r = Dense(nb_dense, activation='tanh')(d9) - d11_r = Dropout(dropout_rate2)(d10_r) - out_r = Dense(1, name='output_2')(d11_r) - - model = Model(inputs=[diff_input, rawf_input], outputs=[out_y, out_r]) + d10 = Dense(nb_dense, activation='tanh')(d9) + d11 = Dropout(dropout_rate2)(d10) + out = Dense(1)(d11) + model = Model(inputs=[diff_input, rawf_input], outputs=out) return model - -# %% TS_CAN - -def TS_CAN(n_frame, nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3), dropout_rate1=0.25, dropout_rate2=0.5, +#%% PTS_CAN --> Advanced TS_CAN with binary output signal +def PTS_CAN(n_frame, nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3), dropout_rate1=0.25, dropout_rate2=0.5, pool_size=(2, 2), nb_dense=128): diff_input = Input(shape=input_shape) rawf_input = Input(shape=input_shape) @@ -193,15 +292,15 @@ def TS_CAN(n_frame, nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3), d d9 = Flatten()(d8) d10 = Dense(nb_dense, activation='tanh')(d9) d11 = Dropout(dropout_rate2)(d10) - out = Dense(1)(d11) - model = Model(inputs=[diff_input, rawf_input], outputs=out) - return model + out1 = Dense(1, name='output_1')(d11) + out_peaks = ownLayer_binaryPeak(name='output_2')(out1) + model = Model(inputs=[diff_input, rawf_input], outputs=[out1, out_peaks]) + return model -# %% MTTS-CAN - -def MTTS_CAN(n_frame, nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3), dropout_rate1=0.25, - dropout_rate2=0.5, pool_size=(2, 2), nb_dense=128): +# Advanced PTS-CAN: with additional parameter calculation +def PPTS_CAN(n_frame, nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3), dropout_rate1=0.25, dropout_rate2=0.5, + pool_size=(2, 2), nb_dense=128, parameter=None): diff_input = Input(shape=input_shape) rawf_input = Input(shape=input_shape) @@ -235,20 +334,17 @@ def MTTS_CAN(n_frame, nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3), d8 = Dropout(dropout_rate1)(d7) d9 = Flatten()(d8) + d10 = Dense(nb_dense, activation='tanh')(d9) + d11 = Dropout(dropout_rate2)(d10) + out1 = Dense(1, name='output_1')(d11) + out_peaks = ownLayer_binaryPeak(name='output_2')(out1) + out_params = ownLayer_parameter(trainable=False, name='output_3')(out_peaks, parameter) + #out_params = ownLayer_parameter(name='output_3')(out_peaks, parameter) - d10_y = Dense(nb_dense, activation='tanh')(d9) - d11_y = Dropout(dropout_rate2)(d10_y) - out_y = Dense(1, name='output_1')(d11_y) - - d10_r = Dense(nb_dense, activation='tanh')(d9) - d11_r = Dropout(dropout_rate2)(d10_r) - out_r = Dense(1, name='output_2')(d11_r) - - model = Model(inputs=[diff_input, rawf_input], outputs=[out_y, out_r]) + model = Model(inputs=[diff_input, rawf_input], outputs=[out1, out_peaks, out_params]) return model - -# %% +# %% --> PhysNet mit AttentionModule def CAN_3D(n_frame, nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3, 3), dropout_rate1=0.25, dropout_rate2=0.5, pool_size=(2, 2, 2), nb_dense=128): diff_input = Input(shape=input_shape) @@ -286,59 +382,7 @@ def CAN_3D(n_frame, nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3, 3) model = Model(inputs=[diff_input, rawf_input], outputs=out) return model - -# input_shape = (36, 36, 10, 3) -# model = DeepPhy_3DCNN(10, 32, 64, input_shape) -# print('==========================') - - # %% -def MT_CAN_3D(n_frame, nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3, 3), dropout_rate1=0.25, - dropout_rate2=0.5, pool_size=(2, 2, 2), nb_dense=128): - diff_input = Input(shape=input_shape) - rawf_input = Input(shape=input_shape) - - d1 = Conv3D(nb_filters1, kernel_size, padding='same', activation='tanh')(diff_input) - d2 = Conv3D(nb_filters1, kernel_size, activation='tanh')(d1) - - # Appearance Branch - r1 = Conv3D(nb_filters1, kernel_size, padding='same', activation='tanh')(rawf_input) - r2 = Conv3D(nb_filters1, kernel_size, activation='tanh')(r1) - g1 = Conv3D(1, (1, 1, 1), padding='same', activation='sigmoid')(r2) - g1 = Attention_mask()(g1) - gated1 = multiply([d2, g1]) - - d3 = AveragePooling3D(pool_size)(gated1) - d4 = Dropout(dropout_rate1)(d3) - d5 = Conv3D(nb_filters2, kernel_size, padding='same', activation='tanh')(d4) - d6 = Conv3D(nb_filters2, kernel_size, activation='tanh')(d5) - - r3 = AveragePooling3D(pool_size)(r2) - r4 = Dropout(dropout_rate1)(r3) - r5 = Conv3D(nb_filters2, kernel_size, padding='same', activation='tanh')(r4) - r6 = Conv3D(nb_filters2, kernel_size, activation='tanh')(r5) - g2 = Conv3D(1, (1, 1, 1), padding='same', activation='sigmoid')(r6) - g2 = Attention_mask()(g2) - gated2 = multiply([d6, g2]) - d7 = AveragePooling3D(pool_size)(gated2) - d8 = Dropout(dropout_rate1)(d7) - - d9 = Flatten()(d8) - d10_y = Dense(nb_dense, activation='tanh')(d9) - d11_y = Dropout(dropout_rate2)(d10_y) - out_y = Dense(n_frame, name='output_1')(d11_y) - - d10_r = Dense(nb_dense, activation='tanh')(d9) - d11_r = Dropout(dropout_rate2)(d10_r) - out_r = Dense(n_frame, name='output_2')(d11_r) - - model = Model(inputs=[diff_input, rawf_input], outputs=[out_y, out_r]) - - return model - - -# %% - def Hybrid_CAN(n_frame, nb_filters1, nb_filters2, input_shape_1, input_shape_2, kernel_size_1=(3, 3, 3), kernel_size_2=(3, 3), dropout_rate1=0.25, dropout_rate2=0.5, pool_size_1=(2, 2, 2), pool_size_2=(2, 2), nb_dense=128): @@ -391,6 +435,145 @@ def Hybrid_CAN(n_frame, nb_filters1, nb_filters2, input_shape_1, input_shape_2, model = Model(inputs=[diff_input, rawf_input], outputs=out) return model +####### Multi Task Models (Blood Volume Pulse and Respiration Rate) +# %% MT_CAN +def MT_CAN(nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3), dropout_rate1=0.25, dropout_rate2=0.5, + pool_size=(2, 2), nb_dense=128): + diff_input = Input(shape=input_shape) + rawf_input = Input(shape=input_shape) + + d1 = Conv2D(nb_filters1, kernel_size, padding='same', activation='tanh')(diff_input) + d2 = Conv2D(nb_filters1, kernel_size, activation='tanh')(d1) + + r1 = Conv2D(nb_filters1, kernel_size, padding='same', activation='tanh')(rawf_input) + r2 = Conv2D(nb_filters1, kernel_size, activation='tanh')(r1) + + g1 = Conv2D(1, (1, 1), padding='same', activation='sigmoid')(r2) + g1 = Attention_mask()(g1) + gated1 = multiply([d2, g1]) + + d3 = AveragePooling2D(pool_size)(gated1) + d4 = Dropout(dropout_rate1)(d3) + + r3 = AveragePooling2D(pool_size)(r2) + r4 = Dropout(dropout_rate1)(r3) + + d5 = Conv2D(nb_filters2, kernel_size, padding='same', activation='tanh')(d4) + d6 = Conv2D(nb_filters2, kernel_size, activation='tanh')(d5) + + r5 = Conv2D(nb_filters2, kernel_size, padding='same', activation='tanh')(r4) + r6 = Conv2D(nb_filters2, kernel_size, activation='tanh')(r5) + + g2 = Conv2D(1, (1, 1), padding='same', activation='sigmoid')(r6) + g2 = Attention_mask()(g2) + gated2 = multiply([d6, g2]) + + d7 = AveragePooling2D(pool_size)(gated2) + d8 = Dropout(dropout_rate1)(d7) + + d9 = Flatten()(d8) + d10_y = Dense(nb_dense, activation='tanh')(d9) + d11_y = Dropout(dropout_rate2)(d10_y) + out_y = Dense(1, name='output_1')(d11_y) + + d10_r = Dense(nb_dense, activation='tanh')(d9) + d11_r = Dropout(dropout_rate2)(d10_r) + out_r = Dense(1, name='output_2')(d11_r) + + model = Model(inputs=[diff_input, rawf_input], outputs=[out_y, out_r]) + return model + +# %% MTTS-CAN +def MTTS_CAN(n_frame, nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3), dropout_rate1=0.25, + dropout_rate2=0.5, pool_size=(2, 2), nb_dense=128): + diff_input = Input(shape=input_shape) + rawf_input = Input(shape=input_shape) + + d1 = TSM_Cov2D(diff_input, n_frame, nb_filters1, kernel_size, padding='same', activation='tanh') + d2 = TSM_Cov2D(d1, n_frame, nb_filters1, kernel_size, padding='valid', activation='tanh') + + r1 = Conv2D(nb_filters1, kernel_size, padding='same', activation='tanh')(rawf_input) + r2 = Conv2D(nb_filters1, kernel_size, activation='tanh')(r1) + + g1 = Conv2D(1, (1, 1), padding='same', activation='sigmoid')(r2) + g1 = Attention_mask()(g1) + gated1 = multiply([d2, g1]) + + d3 = AveragePooling2D(pool_size)(gated1) + d4 = Dropout(dropout_rate1)(d3) + + r3 = AveragePooling2D(pool_size)(r2) + r4 = Dropout(dropout_rate1)(r3) + + d5 = TSM_Cov2D(d4, n_frame, nb_filters2, kernel_size, padding='same', activation='tanh') + d6 = TSM_Cov2D(d5, n_frame, nb_filters2, kernel_size, padding='valid', activation='tanh') + + r5 = Conv2D(nb_filters2, kernel_size, padding='same', activation='tanh')(r4) + r6 = Conv2D(nb_filters2, kernel_size, activation='tanh')(r5) + + g2 = Conv2D(1, (1, 1), padding='same', activation='sigmoid')(r6) + g2 = Attention_mask()(g2) + gated2 = multiply([d6, g2]) + + d7 = AveragePooling2D(pool_size)(gated2) + d8 = Dropout(dropout_rate1)(d7) + + d9 = Flatten()(d8) + + d10_y = Dense(nb_dense, activation='tanh')(d9) + d11_y = Dropout(dropout_rate2)(d10_y) + out_y = Dense(1, name='output_1')(d11_y) + + d10_r = Dense(nb_dense, activation='tanh')(d9) + d11_r = Dropout(dropout_rate2)(d10_r) + out_r = Dense(1, name='output_2')(d11_r) + + model = Model(inputs=[diff_input, rawf_input], outputs=[out_y, out_r]) + return model + +# %% +def MT_CAN_3D(n_frame, nb_filters1, nb_filters2, input_shape, kernel_size=(3, 3, 3), dropout_rate1=0.25, + dropout_rate2=0.5, pool_size=(2, 2, 2), nb_dense=128): + diff_input = Input(shape=input_shape) + rawf_input = Input(shape=input_shape) + + d1 = Conv3D(nb_filters1, kernel_size, padding='same', activation='tanh')(diff_input) + d2 = Conv3D(nb_filters1, kernel_size, activation='tanh')(d1) + + # Appearance Branch + r1 = Conv3D(nb_filters1, kernel_size, padding='same', activation='tanh')(rawf_input) + r2 = Conv3D(nb_filters1, kernel_size, activation='tanh')(r1) + g1 = Conv3D(1, (1, 1, 1), padding='same', activation='sigmoid')(r2) + g1 = Attention_mask()(g1) + gated1 = multiply([d2, g1]) + + d3 = AveragePooling3D(pool_size)(gated1) + d4 = Dropout(dropout_rate1)(d3) + d5 = Conv3D(nb_filters2, kernel_size, padding='same', activation='tanh')(d4) + d6 = Conv3D(nb_filters2, kernel_size, activation='tanh')(d5) + + r3 = AveragePooling3D(pool_size)(r2) + r4 = Dropout(dropout_rate1)(r3) + r5 = Conv3D(nb_filters2, kernel_size, padding='same', activation='tanh')(r4) + r6 = Conv3D(nb_filters2, kernel_size, activation='tanh')(r5) + g2 = Conv3D(1, (1, 1, 1), padding='same', activation='sigmoid')(r6) + g2 = Attention_mask()(g2) + gated2 = multiply([d6, g2]) + d7 = AveragePooling3D(pool_size)(gated2) + d8 = Dropout(dropout_rate1)(d7) + + d9 = Flatten()(d8) + d10_y = Dense(nb_dense, activation='tanh')(d9) + d11_y = Dropout(dropout_rate2)(d10_y) + out_y = Dense(n_frame, name='output_1')(d11_y) + + d10_r = Dense(nb_dense, activation='tanh')(d9) + d11_r = Dropout(dropout_rate2)(d10_r) + out_r = Dense(n_frame, name='output_2')(d11_r) + + model = Model(inputs=[diff_input, rawf_input], outputs=[out_y, out_r]) + + return model # %% def MT_Hybrid_CAN(n_frame, nb_filters1, nb_filters2, input_shape_1, input_shape_2, kernel_size_1=(3, 3, 3), @@ -450,9 +633,8 @@ def MT_Hybrid_CAN(n_frame, nb_filters1, nb_filters2, input_shape_1, input_shape_ model = Model(inputs=[diff_input, rawf_input], outputs=[out_y, out_r]) return model - # %% -class HeartBeat(keras.callbacks.Callback): +class HeartBeat(tf.keras.callbacks.Callback): def __init__(self, train_gen, test_gen, args, cv_split, save_dir): super(HeartBeat, self).__init__() self.train_gen = train_gen diff --git a/code/model_evaluation.py b/code/model_evaluation.py new file mode 100644 index 0000000..396ee1a --- /dev/null +++ b/code/model_evaluation.py @@ -0,0 +1,342 @@ +from aifc import Error +from operator import mod +import numpy as np +import scipy.io +import xlsxwriter +from model import CAN, CAN_3D, PPTS_CAN, PTS_CAN, TS_CAN, Hybrid_CAN +import h5py +import os +import matplotlib.pyplot as plt +from scipy.signal import butter +from inference_preprocess import preprocess_raw_video, detrend +from sklearn.preprocessing import MinMaxScaler +import pandas as pd +from sklearn import metrics +import scipy.stats as sc +from glob import glob +import tensorflow as tf +from tensorflow.python.framework import ops + +import heartpy as hp + +def prepare_3D_CAN(dXsub): + frame_depth = 10 + num_window = int(dXsub.shape[0]) - frame_depth + 1 + tempX = np.array([dXsub[f:f + frame_depth, :, :, :] # (491, 10, 36, 36 ,6) (169, 10, 36, 36, 6) + for f in range(num_window)]) + tempX = np.swapaxes(tempX, 1, 3) # (169, 36, 36, 10, 6) + tempX = np.swapaxes(tempX, 1, 2) # (169, 36, 36, 10, 6) + return tempX + +def prepare_Hybrid_CAN(dXsub): + frame_depth = 10 + num_window = int(dXsub.shape[0]) - frame_depth + 1 + tempX = np.array([dXsub[f:f + frame_depth, :, :, :] # (169, 10, 36, 36, 6) + for f in range(num_window)]) + tempX = np.swapaxes(tempX, 1, 3) # (169, 36, 36, 10, 6) + tempX = np.swapaxes(tempX, 1, 2) # (169, 36, 36, 10, 6) + motion_data = tempX[:, :, :, :, :3] + apperance_data = np.average(tempX[:, :, :, :, -3:], axis=-2) + return motion_data, apperance_data + +def predict_vitals(workBook, test_name, model_name, video_path, path_results): + mms = MinMaxScaler() + img_rows = 36 + img_cols = 36 + frame_depth = 10 + batch_size = 100 + model_checkpoint = None + try: + model_checkpoint = os.path.join(path_results, test_name, "cv_0_epoch24_model.hdf5") + except: + model_checkpoint = os.path.join(path_results, test_name, "cv_0_epoch23_model.hdf5") + batch_size = batch_size + sample_data_path = video_path + print("path: ",sample_data_path) + dXsub, fs = preprocess_raw_video(sample_data_path, dim=36) + print('dXsub shape', dXsub.shape, "fs: ", fs) + + if model_name == "PPTS_CAN": + dXsub_len = (dXsub.shape[0] // (frame_depth*10)) * (frame_depth*10) + dXsub = dXsub[:dXsub_len, :, :, :] + + else: + dXsub_len = (dXsub.shape[0] // frame_depth) * frame_depth + dXsub = dXsub[:dXsub_len, :, :, :] + + if model_name == "TS_CAN": + model = TS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3)) + elif model_name == "3D_CAN": + model = CAN_3D(frame_depth, 32, 64, (img_rows, img_cols, frame_depth, 3)) + dXsub = prepare_3D_CAN(dXsub) + dXsub_len = (dXsub.shape[0] // (frame_depth)) * (frame_depth) + dXsub = dXsub[:dXsub_len, :, :, :,:] + elif model_name == "CAN": + model = CAN(32, 64, (img_rows, img_cols, 3)) + elif model_name == "Hybrid_CAN": + model = Hybrid_CAN(frame_depth, 32, 64, (img_rows, img_cols, frame_depth, 3), + (img_rows, img_cols, 3)) + dXsub1, dXsub2 = prepare_Hybrid_CAN(dXsub) + dXsub_len1 = (dXsub1.shape[0] // (frame_depth)) * (frame_depth) + dXsub1 = dXsub1[:dXsub_len1, :, :, :,:] + dXsub_len2 = (dXsub2.shape[0] // (frame_depth)) * (frame_depth) + dXsub2 = dXsub2[:dXsub_len2, :, :, :] + elif model_name == "PTS_CAN": + model = PTS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3)) + elif model_name == "PPTS_CAN": + model = PPTS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3), parameter=['bpm', 'sdnn']) + else: + raise NotImplementedError + + + model.load_weights(model_checkpoint) + if model_name == "3D_CAN": + yptest = model.predict((dXsub[:, :, :,: , :3], dXsub[:, :, :, : , -3:])) + #yptest = model((dXsub[:, :, :,: , :3], dXsub[:, :, :, : , -3:]), training=False) + elif model_name == "Hybrid_CAN": + #yptest = model.predict((dXsub1, dXsub2), batch_size=batch_size, verbose=1) + yptest = model.predict((dXsub1, dXsub2)) + else: + yptest = model((dXsub[:, :, :, :3], dXsub[:, :, :, -3:]), training=False) #, verbose=1) + # yptest = model.predict((dXsub[:, :, :, :3], dXsub[:, :, :, -3:])) + if model_name == "3D_CAN" or model_name == "Hybrid_CAN": + pulse_pred = yptest[:,0] + elif model_name != "PTS_CAN" and model_name != "PPTS_CAN": + pulse_pred = yptest + + else: + pulse_pred = yptest[0] + + pulse_pred = detrend(np.cumsum(pulse_pred), 100) + [b_pulse_pred, a_pulse_pred] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') + pulse_pred = scipy.signal.filtfilt(b_pulse_pred, a_pulse_pred, np.double(pulse_pred)) + pulse_pred = np.array(mms.fit_transform(pulse_pred.reshape(-1,1))).flatten() + + ##### ground truth data resampled ####### + if(str(sample_data_path).find("COHFACE") > 0): + truth_path = sample_data_path.replace(".avi", "_dataFile.hdf5") + elif(str(sample_data_path).find("UBFC-PHYS") > 0): + truth_path = sample_data_path.replace("vid_", "").replace(".avi","_dataFile.hdf5") + elif(str(sample_data_path).find("UBFC") > 0): + truth_path = sample_data_path.replace("vid.avi", "dataFile.hdf5") + else: + return print("Error in finding the ground truth signal...") + + gound_truth_file = h5py.File(truth_path, "r") + pulse_truth = gound_truth_file["pulse"] ### range ground truth from 0 to 1 + pulse_truth = pulse_truth[0:dXsub_len] + pulse_truth = detrend(np.cumsum(pulse_truth), 100) + [b_pulse_tr, a_pulse_tr] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') + pulse_truth = scipy.signal.filtfilt(b_pulse_tr, a_pulse_tr, np.double(pulse_truth)) + pulse_truth = np.array(mms.fit_transform(pulse_truth.reshape(-1,1))).flatten() + ### same size ####### + if len(pulse_pred) > len(pulse_truth): + pulse_pred = pulse_pred[:len(pulse_truth)] + elif len(pulse_pred) < len(pulse_truth): + pulse_truth = pulse_truth[:len(pulse_pred)] + ########### Peaks ########### + working_data_pred, measures_pred = hp.process(pulse_pred, fs, calc_freq=True) + working_data_truth, measures_truth = hp.process(pulse_truth, fs, calc_freq=True) + peaks_pred = working_data_pred['peaklist'] + peaks_truth = working_data_truth['peaklist'] + + ######## name files ############# + nameStrAll = str(sample_data_path).split("/") + nameStr = "" + if(str(sample_data_path).find("COHFACE") > 0): + for item in nameStrAll[4:6]: + nameStr += item + "-" + elif(str(sample_data_path).find("UBFC-PHYS") > 0): + nameStr = str(nameStrAll[5]).replace("vid", "").replace(".avi", "") + elif(str(sample_data_path).find("UBFC") > 0): + nmr = str(sample_data_path).find("UBFC") + nameStr = str(sample_data_path)[nmr + 5:].replace("\\", "-").replace("vid.avi", "").replace("/", "") + else: + raise ValueError + + ########## Plot ################## + peaks_pred_new = [] + for peak in peaks_pred: + if (peak > 400 and peak <700): + peaks_pred_new.append(peak-400) + peaks_truth_new = [] + for peak in peaks_truth: + if (peak > 400 and peak <700): + peaks_truth_new.append(peak-400) + plt.figure() #subplot(211) + plt.plot(pulse_pred[400:700], "#E6001A", label='rPPG signal') + plt.plot(peaks_truth_new, pulse_truth[400:700][peaks_truth_new], "x", color="#005AA9") + plt.plot(peaks_pred_new, pulse_pred[400:700][peaks_pred_new], "x", color ='#E6001A') + plt.title('rPPG signal with ground truth') + plt.ylabel("normalized Signal [a.u.]") + plt.xlabel("time (samples)") + plt.plot(pulse_truth[400:700], '#005AA9', linewidth=0.9, label='ground truth') + plt.legend() + plt.savefig(nameStr + "_both.svg", format="svg") + + plt.figure() + plt.subplot(211) + plt.plot(pulse_truth[400:700],"#004E8A", label='Ground truth') + plt.plot(peaks_truth_new, pulse_truth[400:700][peaks_truth_new], "x", color="#004E8A") + plt.ylabel("normalized Signal [a.u.]") + plt.title('Ground truth') + plt.subplot(212) + plt.plot(pulse_pred[400:700], "#004E8A",label='Prediction') + plt.plot(peaks_pred_new, pulse_pred[400:700][peaks_pred_new],"x", color="#004E8A") + plt.title("Predicted rPPG") + plt.ylabel("normalized Signal [a.u.]") + plt.xlabel("time (samples)") + plt.legend() + plt.savefig(nameStr) + + ########### IBI ############# + #ibi_truth = working_data_truth['RR_list_cor'] + #print(ibi_truth) + #ibi_pred = working_data_pred['RR_list_cor'] + #print(ibi_pred) + ######### HRV featurs ############## + #print("HRV Truth: ",measures_truth) + #print("HRV Pred: ", measures_pred) + ######### Metrics ############## + # MSE: + MAE = metrics.mean_absolute_error(pulse_truth, pulse_pred) + MSE = metrics.mean_squared_error(pulse_truth, pulse_pred) + # RMSE: + RMSE = metrics.mean_squared_error(pulse_truth, pulse_pred, squared=False) + # Pearson correlation: + p = sc.pearsonr(pulse_truth, pulse_pred) + + ####### Logging ############# + worksheet = workBook.add_worksheet(nameStr) + worksheet.write(0,0, video_path) + worksheet.write(1,0, "MAE") + worksheet.write(1,1, MAE) + worksheet.write(2,0, "RMSE") + worksheet.write(2,1, RMSE) + worksheet.write(3,0, "p") + worksheet.write(3,1, p[0]) + worksheet.write(5,1, "Truth") + worksheet.write(5,2, "Prediction") + worksheet.write(6,0, "bpm") + worksheet.write(6,1, measures_truth["bpm"]) + worksheet.write(6,2, measures_pred["bpm"]) + worksheet.write(7,0, "sdnn") + worksheet.write(7,1, measures_truth["sdnn"]) + worksheet.write(7,2, measures_pred["sdnn"]) + worksheet.write(8,0, "rmssd") + worksheet.write(8,1, measures_truth["rmssd"]) + worksheet.write(8,2, measures_pred["rmssd"]) + worksheet.write(9,0, "pnn50") + worksheet.write(9,1, measures_truth["pnn50"]) + worksheet.write(9,2, measures_pred["pnn50"]) + worksheet.write(10,0, "lf/hf") + worksheet.write(10,1, measures_truth["lf/hf"]) + worksheet.write(10,2, measures_pred["lf/hf"]) + worksheet.write(11,0, "ibi Average") + worksheet.write(11,1, measures_truth["ibi"]) + worksheet.write(11,2, measures_pred["ibi"]) + + worksheet.write(13,0, "pulse_truth") + col = 0 + for val in pulse_truth: + worksheet.write(14, col, val) + col += 1 + worksheet.write(15,0, "pulse_predict") + col = 0 + for val in pulse_pred: + worksheet.write(15, col, val) + col += 1 + +if __name__ == "__main__": + path_results = "D:/Databases/4)Results/Version5" + dir_names = glob(path_results + "/P*") + test_names = [] + for dir in dir_names: + split = dir.split("\\") + test_names.append(split[len(split)-1]) + # video_path = ["D:/Databases/1)Training/COHFACE/5/1/data.avi", + # "D:/Databases/1)Training/COHFACE/10/2/data.avi", "D:/Databases/1)Training/UBFC-PHYS/s5/vid_s5_T1.avi", + # "D:/Databases/1)Training/COHFACE/6/0/data.avi", + # "D:/Databases/1)Training/UBFC-PHYS/s13/vid_s13_T3.avi", + + # "D:/Databases/2)Validation/UBFC-PHYS/s40/vid_s40_T2.avi", "D:/Databases/2)Validation/UBFC-PHYS/s44/vid_s44_T1.avi", + # "D:/Databases/2)Validation/COHFACE/38/0/data.avi", "D:/Databases/2)Validation/UBFC-PHYS/s38/vid_s38_T1.avi", + # "D:/Databases/2)Validation/COHFACE/34/2/data.avi"] + # video_path = ["D:/Databases/1)Training/COHFACE/5/1/data.avi", + # "D:/Databases/1)Training/COHFACE/10/2/data.avi", "C:/Users/sarah/OneDrive/Desktop/UBFC/DATASET_2/subject3/vid.avi", + # "D:/Databases/1)Training/COHFACE/6/0/data.avi", + # "C:/Users/sarah/OneDrive/Desktop/UBFC/DATASET_2/subject15/vid.avi", + + # "C:/Users/sarah/OneDrive/Desktop/UBFC/DATASET_2/subject34/vid.avi", "C:/Users/sarah/OneDrive/Desktop/UBFC/DATASET_2/subject38/vid.avi", + # "D:/Databases/2)Validation/COHFACE/38/0/data.avi", "C:/Users/sarah/OneDrive/Desktop/UBFC/DATASET_2/subject41/vid.avi", + # "D:/Databases/2)Validation/COHFACE/34/2/data.avi"] + + video_path = [\ + "D:/Databases/1)Training/COHFACE/5/1/data.avi", + "D:/Databases/1)Training/COHFACE/10/2/data.avi", + "D:/Databases/2)Validation/COHFACE/38/0/data.avi", + + "D:/Databases/1)Training/UBFC-PHYS/s5/vid_s5_T1.avi", + "D:/Databases/1)Training/UBFC-PHYS/s13/vid_s13_T3.avi", + "D:/Databases/2)Validation/UBFC-PHYS/s38/vid_s38_T1.avi", + + "C:/Users/sarah/OneDrive/Desktop/UBFC/DATASET_2/subject3/vid.avi", + "C:/Users/sarah/OneDrive/Desktop/UBFC/DATASET_2/subject9/vid.avi", + "C:/Users/sarah/OneDrive/Desktop/UBFC/DATASET_2/subject40/vid.avi"] + video_path = [\ + "D:/Databases/3)Testing/COHFACE/21/1/data.avi", + "D:/Databases/3)Testing/COHFACE/25/2/data.avi", + "D:/Databases/3)Testing/COHFACE/28/0/data.avi", + + "D:/Databases/3)Testing/UBFC/subject42/vid.avi", + "D:/Databases/3)Testing/UBFC/subject44/vid.avi", + "D:/Databases/3)Testing/UBFC/subject47/vid.avi"] + + test_names = ['PPTS_CAN_all','PPTS_CAN_sdnn_pnn50_lfhf2', 'PTS_CAN_TE2'] + + save_dir = "D:/Databases/5)Evaluation/Test" + print("Models: ", test_names) + + + for test_name in test_names: + tf.keras.backend.clear_session() + tf.autograph.set_verbosity(10) + ops.reset_default_graph() + print("Current Modelname: ", test_name) + if str(test_name).find("3D_CAN") >=0: + model_name = "3D_CAN" + continue + elif str(test_name).find("Hybrid_CAN") >= 0: + model_name = "Hybrid_CAN" + continue + elif str(test_name).find("PPTS") >= 0: + model_name = "PPTS_CAN" + elif str(test_name).find("PTS") >= 0: + model_name = "PTS_CAN" + elif str(test_name).find("TS_CAN") >= 0: + model_name = "TS_CAN" + else: + if str(test_name).find("CAN") >= 0: + model_name = "CAN" + continue + else: + raise Error("Model not found...") + + # neuer Ordner für Tests + os.chdir(save_dir) + try: + os.makedirs(str(test_name)) + except: + print("Directory exists...") + save_path = os.path.join(save_dir, str(test_name)) + os.chdir(save_path) + workbook = xlsxwriter.Workbook(test_name + ".xlsx") + for vid in video_path: + predict_vitals(workbook, test_name, model_name, vid, path_results) + print("Ready with this model") + workbook.close() + +#python code/predict_vitals_new.py --video_path "D:\Databases\1)Training\COHFACE\1\1\data.avi" --trained_model ./cv_0_epoch24_model.hdf5 +#./rPPG-checkpoints/testCohFace1/cv_0_epoch24_model.hdf5 +#./rPPG-checkpoints/test1/cv_0_epoch04_model.hdf5' +#python code/predict_vitals_new.py --video_path "D:\Databases\1)Training\UBFC-PHYS\s1\vid_s1_T1.avi" --trained_model ./cv_0_epoch24_model.hdf5 \ No newline at end of file diff --git a/code/pre_process.py b/code/pre_process.py index 0516865..6f0efb3 100644 --- a/code/pre_process.py +++ b/code/pre_process.py @@ -1,9 +1,11 @@ import glob +import itertools import os import h5py import numpy as np import scipy.io +import pandas as pd def get_nframe_video(path): @@ -12,6 +14,12 @@ def get_nframe_video(path): nframe_per_video = temp_dysub.shape[0] return nframe_per_video +def get_nframe_video_(path): + temp_f1 = h5py.File(path, 'r') + temp_data = np.array(temp_f1["data"]) + nframe_per_video = temp_data.shape[0] + return nframe_per_video + def get_nframe_video_val(path): temp_f1 = scipy.io.loadmat(path) @@ -20,10 +28,12 @@ def get_nframe_video_val(path): return nframe_per_video -def split_subj(data_dir, cv_split, subNum): - f3 = h5py.File(data_dir + '/M.mat', 'r') - M = np.transpose(np.array(f3["M"])).astype(np.bool) - subTrain = subNum[~M[:, cv_split]].tolist() +def split_subj(data_dir, cv_split, subNum): # trennen der Daten innerhalb 1 Subjekts... + print(subNum) + f3 = h5py.File( data_dir +'/s1/bvp_s1_T1.csv', 'r')# "/testSub.mat" + # f4 = pd.read_csv(data_dir + '/s1/bvp_s1_T1.csv') + M = np.transpose(np.array(f3["M"])).astype(np.bool) #? wieso als bool? + subTrain = subNum[~M[:, cv_split]].tolist() # wieso nur cv_split? subTest = subNum[M[:, cv_split]].tolist() return subTrain, subTest @@ -45,3 +55,184 @@ def sort_video_list(data_dir, taskList, subTrain): x = sorted(x, key=take_last_ele) final.append(x) return final + + +def sort_video_list_(data_dir, taskList, subTrain, database_name, train): + final = [] + if database_name == "UBFC_PHYS": + if train: + for p in subTrain: + + x = glob.glob(os.path.join(data_dir, '1)Training/UBFC-PHYS/s' + str(p), 'vid_s*')) + x = sorted(x) + #x = sorted(x, key=take_last_ele) + final.append(x) + else: + for p in subTrain: + x = glob.glob(os.path.join(data_dir, '2)Validation/UBFC-PHYS/s' + str(p), 'vid_s*')) + x = sorted(x) + #x = sorted(x, key=take_last_ele) + final.append(x) + + elif database_name == "COHFACE": + if train: + for p in subTrain: + for t in taskList: + x = glob.glob(os.path.join(data_dir, '1)Training/COHFACE/', str(p), str(t), 'data.avi')) + x = sorted(x) + #x = sorted(x, key=take_last_ele) + final.append(x) + else: + for p in subTrain: + for t in taskList: + x = glob.glob(os.path.join(data_dir, '2)Validation/COHFACE/', str(p), str(t), 'data.avi')) + x = sorted(x) + #x = sorted(x, key=take_last_ele) + final.append(x) + elif database_name == "UBFC": + if train: + x = glob.glob(os.path.join(data_dir, "**/", 'vid.avi'), recursive=True) + x = sorted(x) + #x = sorted(x, key=take_last_ele) + final.append(x) + else: + print("not implemented yet.") + return final + +def sort_dataFile_list_(data_dir, subTrain, database_name, trainMode): + if database_name == "UBFC_PHYS": + final = dataFiles_UBFC_PHYS(data_dir, subTrain, trainMode, mode=1) + final = list(itertools.chain(*final)) + elif database_name == "COHFACE": + taskList = [0, 1, 2, 3] + final = dataFile_COHFACE(data_dir, taskList, subTrain, trainMode) + final = list(itertools.chain(*final)) + elif database_name == "UBFC": + final = dataFile_UBFC(data_dir, trainMode) + final = list(itertools.chain(*final)) + elif database_name == "MIX1": + for database in subTrain.keys(): + final = [] + if(str(database).find("UBFC") >= 0): + finalPart1 = dataFiles_UBFC_PHYS(data_dir, subTrain[database], trainMode, mode=0) + finalPart1 = list(itertools.chain(*finalPart1)) + elif str(database).find("COHFACE") >= 0: + taskList = [0, 1, 2, 3] + finalPart2 = dataFile_COHFACE(data_dir, taskList, subTrain[database], trainMode) + finalPart2 = list(itertools.chain(*finalPart2)) + else: + raise NotImplementedError + final = finalPart1 + finalPart2 + elif database_name == "MIX2": + for database in subTrain.keys(): + final = [] + if(str(database).find("UBFC") >= 0): + finalPart1 = dataFile_UBFC(data_dir, trainMode) + finalPart1 = list(itertools.chain(*finalPart1)) + elif str(database).find("COHFACE") >= 0: + taskList = [0, 1, 2, 3] + finalPart2 = dataFile_COHFACE(data_dir, taskList, subTrain[database], trainMode) + finalPart2 = list(itertools.chain(*finalPart2)) + else: + raise NotImplementedError + final = finalPart1 + finalPart2 + else: + print("not implemented yet.") + return final + +def dataFile_COHFACE(data_dir, taskList, subTrain, train): + final = [] + if train: + for p in subTrain: + for t in taskList: + x = glob.glob(os.path.join(data_dir, '1)Training/COHFACE', str(p), str(t), '*dataFile.hdf5')) + x = sorted(x) + #x = sorted(x, key=take_last_ele) + final.append(x) + else: + for p in subTrain: + for t in taskList: + x = glob.glob(os.path.join(data_dir, '2)Validation/COHFACE/', str(p), str(t), '*dataFile.hdf5')) + x = sorted(x) + #x = sorted(x, key=take_last_ele) + final.append(x) + return final + +def dataFiles_UBFC_PHYS(data_dir, subTrain, train, mode): + final = [] + if mode == 1: + if train: + for p in subTrain: + x = glob.glob(os.path.join(data_dir, '1)Training/UBFC-PHYS/s' + str(p), "s" + str(p) + "*")) + x = sorted(x) + #x = sorted(x, key=take_last_ele) + final.append(x) + else: + for p in subTrain: + x = glob.glob(os.path.join(data_dir, '2)Validation/UBFC-PHYS/s' + str(p), "s" + str(p) + "*")) + x = sorted(x) + #x = sorted(x, key=take_last_ele) + final.append(x) + else: + if train: + for p in subTrain: + x = glob.glob(os.path.join(data_dir, '1)Training/UBFC-PHYS/' + str(p), str(p) + "*")) + x = sorted(x) + #x = sorted(x, key=take_last_ele) + final.append(x) + else: + for p in subTrain: + x = glob.glob(os.path.join(data_dir, '2)Validation/UBFC-PHYS/' + str(p), str(p) + "*")) + x = sorted(x) + #x = sorted(x, key=take_last_ele) + final.append(x) + return final + +def dataFile_UBFC(data_dir, train): + final = [] + if train: + x = glob.glob(os.path.join(data_dir,'1)Training/UBFC', "**/", 'dataFile.hdf5'), recursive=True) + x = sorted(x) + #x = sorted(x, key=take_last_ele) + final.append(x) + else: + x = glob.glob(os.path.join(data_dir,'2)Validation/UBFC', "**/", 'dataFile.hdf5'), recursive=True) + x = sorted(x) + #x = sorted(x, key=take_last_ele) + final.append(x) + return final + + +def split_subj_(data_dir, database): # trennen der Daten innerhalb 1 Subjekts... + if database == "UBFC_PHYS": + subTrain = np.array(range(1, 37)).tolist() #,37)).tolist() + subTest = np.array(range(37,57)).tolist() + elif database == "COHFACE": + subTrain = np.array(range(1, 33)).tolist()# 33)).tolist() + subTest = np.array(range(32,41)).tolist() # 41)).tolist() + elif database == "UBFC": + subTrain = np.array(range(1,34)) + subTest = np.array([range(34,42)]) + else: + print("This Database isn't implemented yet.") + return subTrain, subTest + +def collect_subj(data_dir, database_name): # collecting all subject out of data_dir.. + + path_tr = os.path.join(data_dir, "1)Training") + path_val = os.path.join(data_dir, "2)Validation") + mix1 = ['COHFACE', 'UBFC-PHYS'] + mix2 = ['COHFACE', 'UBFC'] + if database_name == "MIX1": + mix = mix1 + elif database_name == "MIX2": + mix = mix2 + subTrain = {} + subTest = {} + for database in mix: + subj_tr = os.listdir(os.path.join(path_tr, database)) + subTrain[database] = subj_tr + subj_val = os.listdir(os.path.join(path_val, database)) + subTest[database] = subj_val + + return subTrain, subTest diff --git a/code/predict_vitals.py b/code/predict_vitals.py index 1a2b74e..b7bdab1 100644 --- a/code/predict_vitals.py +++ b/code/predict_vitals.py @@ -1,63 +1,136 @@ +from scipy import signal import tensorflow as tf import numpy as np import scipy.io -import os import sys import argparse sys.path.append('../') -from model import Attention_mask, MTTS_CAN +from model import TS_CAN import h5py import matplotlib.pyplot as plt from scipy.signal import butter from inference_preprocess import preprocess_raw_video, detrend +from hrvanalysis import get_time_domain_features, get_frequency_domain_features def predict_vitals(args): img_rows = 36 img_cols = 36 frame_depth = 10 - model_checkpoint = './mtts_can.hdf5' + #model_checkpoint = './mtts_can.hdf5' + model_checkpoint = args.trained_model batch_size = args.batch_size - fs = args.sampling_rate sample_data_path = args.video_path + print("path: ",sample_data_path) - dXsub = preprocess_raw_video(sample_data_path, dim=36) + dXsub, fs = preprocess_raw_video(sample_data_path, dim=36) print('dXsub shape', dXsub.shape) dXsub_len = (dXsub.shape[0] // frame_depth) * frame_depth dXsub = dXsub[:dXsub_len, :, :, :] - model = MTTS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3)) + model = TS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3)) model.load_weights(model_checkpoint) yptest = model.predict((dXsub[:, :, :, :3], dXsub[:, :, :, -3:]), batch_size=batch_size, verbose=1) - pulse_pred = yptest[0] + pulse_pred = yptest#[0] + pulse_pred = (pulse_pred - pulse_pred.min())/(pulse_pred.max() - pulse_pred.min()) * 2 -1 pulse_pred = detrend(np.cumsum(pulse_pred), 100) [b_pulse, a_pulse] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') pulse_pred = scipy.signal.filtfilt(b_pulse, a_pulse, np.double(pulse_pred)) + + ##### ground truth data resampled ####### + if(str(sample_data_path).find("COHFACE") >= 0): + truth_path = args.video_path.replace(".avi", "_dataFile.hdf5") # akutell für COHACE... + elif(str(sample_data_path).find("UBFC-PHYS")>= 0): + truth_path = args.video_path.replace("vid_", "").replace(".avi","_dataFile.hdf5") + elif(str(sample_data_path).find("UBFC")>= 0): + truth_path = args.video_path.replace("vid.avi","dataFile.hdf5") + else: + return print("Error in finding the ground truth signal...") + gound_truth_file = h5py.File(truth_path, "r") + pulse_truth = gound_truth_file["pulse"] + pulse_truth = detrend(np.cumsum(pulse_truth), 100) + [b_pulse_tr, a_pulse_tr] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') + pulse_truth = scipy.signal.filtfilt(b_pulse_tr, a_pulse_tr, np.double(pulse_truth)) + ### range ground truth from -1 to 1 + pulse_truth = (pulse_truth - pulse_truth.min())/(pulse_truth.max() - pulse_truth.min()) * 2 -1 + pulse_truth = pulse_truth[0: dXsub_len] + + #pulse_pred = pulse_pred[5:] + ########### Peaks ########### + peaks_truth, peaks_ = np.array(signal.find_peaks(pulse_truth, prominence=0.5)) + peaks_pred, b = np.array(signal.find_peaks(pulse_pred, prominence=0.2)) - resp_pred = yptest[1] - resp_pred = detrend(np.cumsum(resp_pred), 100) - [b_resp, a_resp] = butter(1, [0.08 / fs * 2, 0.5 / fs * 2], btype='bandpass') - resp_pred = scipy.signal.filtfilt(b_resp, a_resp, np.double(resp_pred)) + ######## x-axis: time ######### + duration_vid = dXsub_len/fs + x_axis = np.linspace(0, duration_vid, dXsub_len) ########## Plot ################## - plt.subplot(211) - plt.plot(pulse_pred) + plt.figure() #subplot(211) + plt.plot(x_axis, pulse_pred, label='Prediction', color ='#E6001A') + plt.plot(x_axis[peaks_truth], pulse_truth[peaks_truth], "x", color="#721085") + plt.plot(x_axis[peaks_pred], pulse_pred[peaks_pred], "x", color ='#E6001A') plt.title('Pulse Prediction') + plt.xlabel("time (s)") + plt.ylabel("normalized Signal [a.u.]") + plt.plot(x_axis, pulse_truth, label='ground truth', color="#721085") + plt.legend() + + plt.figure() + plt.subplot(211) + plt.plot(x_axis, pulse_truth, label='Ground truth') + plt.plot(x_axis[peaks_truth], pulse_truth[peaks_truth], "x") + plt.ylabel("normalized Signal [a.u.]") + plt.xlabel("time (s)") + plt.title('Ground truth') plt.subplot(212) - plt.plot(resp_pred) - plt.title('Resp Prediction') + plt.plot(x_axis, pulse_pred, label='Prediction') + plt.plot(x_axis[peaks_pred], pulse_pred[peaks_pred], "x") + plt.title("Prediction") + plt.ylabel("normalized Signal [a.u.]") + plt.xlabel("time (s)") + plt.legend() plt.show() + ########### IBI ############# + ibi_truth = np.diff(peaks_truth)*(1000/fs) + print(ibi_truth) + ibi_pred = np.diff(peaks_pred)*(1000/fs) + print(ibi_pred) + ######### HRV featurs ############## + time_domain_features = get_time_domain_features(ibi_truth) + time_domain_features_pred = get_time_domain_features(ibi_pred) + print(time_domain_features) + print(time_domain_features_pred) + freq_domain_features = get_frequency_domain_features(ibi_truth) + freq_domain_features_pred = get_frequency_domain_features(ibi_pred) + print(freq_domain_features) + print(freq_domain_features_pred) + ####### Logging ############# + # neuer Ordner für Tests + file = open(str(sample_data_path).replace(".avi", "_result.txt"),"w") + file.write("LogFile\n\n") + file.write("IBI: "), file.write(str(ibi_pred)) + file.write("\nTime-domain features: "), file.write(str(time_domain_features_pred)) + file.write("\nFrequency-domain features: "), file.write(str(freq_domain_features_pred)) + + file.write("\n\n\nGround truth infos!") + file.write("\nTime-domain features: "), file.write(str(time_domain_features)) + file.write("\nFrequency-domain features: "), file.write(str(freq_domain_features)) if __name__ == "__main__": parser = argparse.ArgumentParser() parser.add_argument('--video_path', type=str, help='processed video path') - parser.add_argument('--sampling_rate', type=int, default = 30, help='sampling rate of your video') parser.add_argument('--batch_size', type=int, default = 100, help='batch size (multiplier of 10)') + parser.add_argument('--trained_model', type=str, default = "D:/Databases/4)Results/Version4/TS_Databases/TS_CAN_COHFACE_2GPU/cv_0_epoch24_model.hdf5", help='path to trained model') + args = parser.parse_args() predict_vitals(args) + +#python code/predict_vitals.py --video_path "C:\Users\sarah\OneDrive\Desktop\UBFC\DATASET_2\subject4\vid.avi" --trained_model "D:\Databases\4)Results\Version4\TS_Databases\TS_CAN_COHFACE_2GPU\cv_0_epoch23_model.hdf5" +#./rPPG-checkpoints/testCohFace1/cv_0_epoch24_model.hdf5 +#./rPPG-checkpoints/test1/cv_0_epoch04_model.hdf5' \ No newline at end of file diff --git a/code/predict_vitals_comparison.py b/code/predict_vitals_comparison.py new file mode 100644 index 0000000..9448b64 --- /dev/null +++ b/code/predict_vitals_comparison.py @@ -0,0 +1,271 @@ +from msilib.schema import Error +import numpy as np +import scipy.io +import sys +import argparse +from losses import filt_peaks +sys.path.append('../') +from model import CAN_3D, PPTS_CAN, PTS_CAN, Attention_mask, MTTS_CAN, TS_CAN +import h5py +import matplotlib.pyplot as plt +from scipy.signal import butter +from inference_preprocess import preprocess_raw_video, detrend +from sklearn.preprocessing import MinMaxScaler + +import heartpy as hp +import os + +def prepare_3D_CAN(dXsub): + frame_depth = 10 + num_window = int(dXsub.shape[0]) - frame_depth + 1 + tempX = np.array([dXsub[f:f + frame_depth, :, :, :] # (491, 10, 36, 36 ,6) (169, 10, 36, 36, 6) + for f in range(num_window)]) + tempX = np.swapaxes(tempX, 1, 3) # (169, 36, 36, 10, 6) + tempX = np.swapaxes(tempX, 1, 2) # (169, 36, 36, 10, 6) + return tempX + +def filt_peaks(x,y): + max_offset = 10 + mask = [] + + for index in range(len(y)-1): + item = y[index] + diff = np.abs(x - item) # diff of truth data and item + min = np.min(diff) # minimum of diff + if min >= max_offset: + mask.append(index) + y = np.delete(y, mask) + mask = [] + for index in range(len(x)-1): + item = x[index] + diff = np.abs(y - item) # diff of truth data and item + min = np.min(diff) # minimum of diff + if min >= max_offset: + mask.append(index) + x = np.delete(x, mask) + return x, y + +def temp_loss(x,y): + x = np.array(x) + y = np.array(y) + diff = np.abs(x-y) + lossframes = np.sum(diff) + loss = lossframes*0.05 + return loss + +def predict_vitals(args): + mms = MinMaxScaler() + img_rows = 36 + img_cols = 36 + frame_depth = 10 + batch_size = args.batch_size + sample_data_path = args.video_path + print("path: ",sample_data_path) + + ts_can_COH = os.path.join(args.trained_model, '3D_CAN_MIX2\\cv_0_epoch24_model.hdf5' ) + ts_can_UB_Ph = os.path.join(args.trained_model,"PPTS_CAN_negPea_TE_sdnn_pnn50_lfhf/cv_0_epoch24_model.hdf5") + ts_can_MIX1 = os.path.join(args.trained_model, "TS_CAN/cv_0_epoch24_model.hdf5") + ts_can_UBFC = os.path.join(args.trained_model, 'TS_CAN\\cv_0_epoch24_model.hdf5' ) + ts_can_MIX2 = os.path.join(args.trained_model,"PPTS_CAN_negPea_gauss_sdnn_pnn50_lfhf/cv_0_epoch24_model.hdf5") + + dXsub, fs = preprocess_raw_video(sample_data_path, dim=36) + print("PROCESSES", fs) + dXsub_len = (dXsub.shape[0] // frame_depth) * frame_depth + dXsub = dXsub[:dXsub_len, :, :, :] + + model_COH = CAN_3D(frame_depth, 32, 64, (img_rows, img_cols, frame_depth, 3)) + model_COH.load_weights(ts_can_COH) + model_MIX1 = TS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3)) + model_MIX1.load_weights(ts_can_MIX1) + model_UB_Ph = PPTS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3), parameter=['bpm', 'sdnn']) + model_UB_Ph.load_weights(ts_can_UB_Ph) + model_MIX2 = PPTS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3), parameter=['bpm', 'sdnn']) + model_MIX2.load_weights(ts_can_MIX2) + model_UBFC = TS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3)) + model_UBFC.load_weights(ts_can_UBFC) + dXsub_3D = prepare_3D_CAN(dXsub) + dXsub_len_3D = (dXsub.shape[0] // (frame_depth)) * (frame_depth) + dXsub_3D = dXsub_3D[:dXsub_len_3D, :, :, :,:] + yptest_COH = model_COH.predict((dXsub_3D[:, :, :,: , :3], dXsub_3D[:, :, :, : , -3:]), verbose=1) + + #yptest_COH = model_COH((dXsub[:, :, :, :3], dXsub[:, :, :, -3:]), training=False) + yptest_COH = yptest_COH[:,0] + yptest_MIX1 = model_MIX1((dXsub[:, :, :, :3], dXsub[:, :, :, -3:]), training=False) + yptest_UB_Ph = model_UB_Ph((dXsub[:, :, :, :3], dXsub[:, :, :, -3:]), training=False) + yptest_UB_Ph = yptest_UB_Ph[0] + yptest_MIX2 = model_MIX2((dXsub[:, :, :, :3], dXsub[:, :, :, -3:]), training=False) + yptest_MIX2 = yptest_MIX2[0] + yptest_UBFC = model_UBFC((dXsub[:, :, :, :3], dXsub[:, :, :, -3:]), training=False) + + pulse_pred_COH = detrend(np.cumsum(yptest_COH), 100) + [b_pulse_pred, a_pulse_pred] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') + pulse_pred_COH = scipy.signal.filtfilt(b_pulse_pred, a_pulse_pred, np.double(pulse_pred_COH)) + pulse_pred_COH = np.array(mms.fit_transform(pulse_pred_COH.reshape(-1,1))).flatten() + + pulse_pred_MIX1 = detrend(np.cumsum(yptest_MIX1), 100) + [b_pulse_pred, a_pulse_pred] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') + pulse_pred_MIX1 = scipy.signal.filtfilt(b_pulse_pred, a_pulse_pred, np.double(pulse_pred_MIX1)) + pulse_pred_MIX1 = np.array(mms.fit_transform(pulse_pred_MIX1.reshape(-1,1))).flatten() + + pulse_pred_UB_Ph = detrend(np.cumsum(yptest_UB_Ph), 100) + [b_pulse_pred, a_pulse_pred] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') + pulse_pred_UB_Ph = scipy.signal.filtfilt(b_pulse_pred, a_pulse_pred, np.double(pulse_pred_UB_Ph)) + pulse_pred_UB_Ph = np.array(mms.fit_transform(pulse_pred_UB_Ph.reshape(-1,1))).flatten() + + pulse_pred_MIX2 = detrend(np.cumsum(yptest_MIX2), 100) + [b_pulse_pred, a_pulse_pred] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') + pulse_pred_MIX2 = scipy.signal.filtfilt(b_pulse_pred, a_pulse_pred, np.double(pulse_pred_MIX2)) + pulse_pred_MIX2 = np.array(mms.fit_transform(pulse_pred_MIX2.reshape(-1,1))).flatten() + + pulse_pred_UBFC = detrend(np.cumsum(yptest_UBFC), 100) + [b_pulse_pred, a_pulse_pred] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') + pulse_pred_UBFC = scipy.signal.filtfilt(b_pulse_pred, a_pulse_pred, np.double(pulse_pred_UBFC)) + pulse_pred_UBFC = np.array(mms.fit_transform(pulse_pred_UBFC.reshape(-1,1))).flatten() + + verk_path = sample_data_path.replace(".avi", "_CHROM.txt") + pulse_pred_verk = open(verk_path, 'r').read() + pulse_pred_verk = str(pulse_pred_verk).split("\n") + pulse_pred_verk = pulse_pred_verk[:-1] + pulse_pred_verk = np.array(list(map(float, pulse_pred_verk))) + + mms = MinMaxScaler() + mean = pulse_pred_verk.mean() + std = np.std(pulse_pred_verk) + upper_limit = mean + std*3 + lower_limit = mean - std*3 + for x in range(0, len(pulse_pred_verk)): + if pulse_pred_verk[x] > upper_limit: + pulse_pred_verk[x] = upper_limit + elif pulse_pred_verk[x] < lower_limit: + pulse_pred_verk[x] = lower_limit + pulse_pred_verk = np.array(mms.fit_transform(pulse_pred_verk.reshape(-1,1))).flatten() # normalization + + ##### ground truth data resampled ####### + if(str(sample_data_path).find("COHFACE") >=0): + truth_path = args.video_path.replace(".avi", "_dataFile.hdf5") + elif(str(sample_data_path).find("UBFC-PHYS") >= 0): + truth_path = args.video_path.replace("vid_", "").replace(".avi","_dataFile.hdf5") + elif(str(sample_data_path).find("UBFC") > 0): + truth_path = sample_data_path.replace("vid.avi", "dataFile.hdf5") + else: + raise ValueError("Error in finding the ground truth signal...") + gound_truth_file = h5py.File(truth_path, "r") + pulse_truth = gound_truth_file["pulse"] ### range ground truth from 0 to 1 + pulse_truth = detrend(np.cumsum(pulse_truth), 100) + [b_pulse_tr, a_pulse_tr] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') + pulse_truth = scipy.signal.filtfilt(b_pulse_tr, a_pulse_tr, np.double(pulse_truth)) + pulse_truth = np.array(mms.fit_transform(pulse_truth.reshape(-1,1))).flatten() + + ########### Peaks ########### + working_data_pred_COH, measures_pred_COH = hp.process(pulse_pred_COH, fs, calc_freq=True) + working_data_pred_MIX1, measures_pred_MIX1 = hp.process(pulse_pred_MIX1, fs, calc_freq=True) + working_data_pred_UB_Ph, measures_pred_UB_Ph = hp.process(pulse_pred_UB_Ph, fs, calc_freq=True) + working_data_pred_MIX2, measures_pred_MIX2 = hp.process(pulse_pred_MIX2, fs, calc_freq=True) + working_data_pred_UBFC, measures_pred_UBFC = hp.process(pulse_pred_UBFC, fs, calc_freq=True) + working_data_Verk, measures_verk = hp.process(pulse_pred_verk, fs, calc_freq=True) + working_data_truth, measures_truth = hp.process(pulse_truth, fs, calc_freq=True) + + peaks_pred_COH = working_data_pred_COH['peaklist'] + peaks_pred_MIX1 = working_data_pred_MIX1['peaklist'] + peaks_pred_UB_Ph = working_data_pred_UB_Ph['peaklist'] + peaks_pred_MIX2 = working_data_pred_MIX2['peaklist'] + peaks_pred_UBFC = working_data_pred_UBFC['peaklist'] + peaks_truth = working_data_truth['peaklist'] + peaks_verk = working_data_Verk['peaklist'] + + ############## loss ############### + #peakCOH, peak_true = filt_peaks(peaks_pred_COH, peaks_truth) + #loss_coh = temp_loss(peak_true, peakCOH) + + ########## Plot ################## + print("FIGURE") + plt.figure() #subplot(211) + plt.plot(pulse_pred_COH, linewidth=1.1, color="#E6001A", label="$\mathrm{3D-CAN}$") + plt.plot(peaks_pred_COH, pulse_pred_COH[peaks_pred_COH], "x", color="#E6001A") + plt.plot(pulse_pred_UB_Ph,linewidth=1.1, color="#721085", label='$\mathrm{PPTS-CAN}_{\mathrm{negPea/TE/sdnn/pNN50/lfhf}}$') + plt.plot(peaks_pred_UB_Ph, pulse_pred_UB_Ph[peaks_pred_UB_Ph], "x", color="#721085") + #.plot(pulse_pred_UBFC, "-.", color="dimgrey", linewidth=1.1, label='$\mathrm{TS-CAN}_{\mathrm{UBFC}}$') + #plt.plot(peaks_pred_UBFC, pulse_pred_UBFC[peaks_pred_UBFC], "x", color="dimgrey") + plt.plot(pulse_pred_MIX2,color="#F5A300", linewidth=1.1, label='$\mathrm{PPTS-CAN}_{\mathrm{negPea/ownGauss/sdnn/pNN50/lfhf}}$') + plt.plot(peaks_pred_MIX2, pulse_pred_MIX2[peaks_pred_MIX2], "x", color="#F5A300") + plt.plot(pulse_pred_verk,color="#99C000", linewidth=1.1, label='$\mathrm{CHROM}$') + plt.plot(peaks_verk, pulse_pred_verk[peaks_verk], "x", color="#99C000") + + plt.ylabel("normalized Signal [a.u.]") + plt.xlabel("time (samples)") + + plt.plot(peaks_truth, pulse_truth[peaks_truth], "x", color="#005AA9") + plt.plot(pulse_truth, "#005AA9", label='ground truth') + plt.legend(loc="lower right") + plt.show() + + print("\n3D: ", measures_pred_COH) + print("\nPTS: ", measures_pred_UB_Ph) + print("\nPPTS: ", measures_pred_MIX2) + print("\nVerk: ", measures_verk) + print("\ntruth: ", measures_truth) + + # plt.figure() + # plt.subplot(211) + # plt.plot(pulse_truth, label='Ground truth') + # plt.plot(peaks_truth, pulse_truth[peaks_truth], "x") + # plt.ylabel("normalized Signal") + # plt.title('Ground truth') + # plt.subplot(212) + # plt.plot(pulse_pred_COHFACE, label='Prediction') + # plt.plot(peaks_pred_COH, pulse_pred_COHFACE[peaks_pred_COH], "x") + # plt.title("Prediction") + # plt.ylabel("normalized Signal") + # plt.legend() + # plt.show() + + ########### IBI ############# + # ibi_truth = working_data_truth['RR_list_cor'] + # print(ibi_truth) + # ibi_pred_COH = working_data_pred_COH['RR_list_cor'] + # print(ibi_pred_COH) + # ibi_pred_MIX = working_data_pred_MIX['RR_list_cor'] + # print(ibi_pred_MIX) + # ibi_pred_UBFC = working_data_pred_UBFC['RR_list_cor'] + # print(ibi_pred_UBFC) + # ######### HRV featurs ############## + # print("HRV Truth: ",measures_truth) + # print("HRV Pred COHFACE: ", measures_pred_COH) + # print("HRV Pred MIX: ", measures_pred_MIX) + # print("HRV Pred UBFC: ", measures_pred_UBFC) + ####### Logging ############# + # neuer Ordner für Tests + # file = open(str(sample_data_path).replace(".avi", "comparisonALL_result.txt"),"w") + # file.write("LogFile\n\n") + # file.write("\nCOHFACE:") + # file.write("\nIBI: "), file.write(str(ibi_pred_COH)) + # file.write("\nHR and HRVfeatures: "), file.write(str(measures_pred_COH)) + + # file.write("\nMIX:") + # file.write("\nIBI: "), file.write(str(ibi_pred_MIX)) + # file.write("\nHR and HRVfeatures: "), file.write(str(measures_pred_MIX)) + + # file.write("\nUBFC-PHYS:") + # file.write("\nIBI: "), file.write(str(ibi_pred_UBFC)) + # file.write("\nHR and HRVfeatures: "), file.write(str(measures_pred_UBFC)) + + # file.write("\n\n\nGround truth infos!") + # file.write("\nHR and HRV features: "), file.write(str(measures_truth)) + +if __name__ == "__main__": + + parser = argparse.ArgumentParser() + parser.add_argument('--video_path', type=str, help='processed video path') + parser.add_argument('--batch_size', type=int, default = 100, help='batch size (multiplier of 10)') + parser.add_argument('--trained_model', type=str, default = './rPPG-checkpoints/testCohFace1/cv_0_epoch24_model.hdf5', help='path to trained model') + + args = parser.parse_args() + + predict_vitals(args) + + +#python code/predict_vitals_comparison.py --video_path "D:/Databases/1)Training/COHFACE/5/1/data.avi" --trained_model "D:\Databases\4)Results\Version4\TS_Databases" +#python code/predict_vitals_comparison.py --video_path "D:/Databases/1)Training/UBFC-PHYS/s5/vid_s5_T1.avi" --trained_model "D:\Databases\4)Results\Version4\TS_Databases" +#python code/predict_vitals_comparison.py --video_path "C:/Users/sarah/OneDrive/Desktop/UBFC/DATASET_2/subject3/vid.avi" --trained_model "D:\Databases\4)Results\Version4\TS_Databases" + +#python code/predict_vitals_comparison.py --video_path "D:/Databases/3)Testing/UBFC/subject42/vid.avi" --trained_model "D:\Databases\4)Results\Version5\" \ No newline at end of file diff --git a/code/predict_vitals_new.py b/code/predict_vitals_new.py new file mode 100644 index 0000000..669499f --- /dev/null +++ b/code/predict_vitals_new.py @@ -0,0 +1,135 @@ +from scipy import signal +import tensorflow as tf +import numpy as np +import scipy.io +import sys +import argparse +sys.path.append('../') +from model import Attention_mask, MTTS_CAN, TS_CAN +import h5py +import matplotlib.pyplot as plt +from scipy.signal import butter +from inference_preprocess import preprocess_raw_video, detrend +from sklearn.preprocessing import MinMaxScaler + +import heartpy as hp + +def predict_vitals(args): + mms = MinMaxScaler() + img_rows = 36 + img_cols = 36 + frame_depth = 10 + #model_checkpoint = './mtts_can.hdf5' + model_checkpoint = args.trained_model + batch_size = args.batch_size + sample_data_path = args.video_path + print("path: ",sample_data_path) + + dXsub, fs = preprocess_raw_video(sample_data_path, dim=36) + print('dXsub shape', dXsub.shape) + + dXsub_len = (dXsub.shape[0] // frame_depth) * frame_depth + dXsub = dXsub[:dXsub_len, :, :, :] + + model = TS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3)) + model.load_weights(model_checkpoint) + + yptest = model.predict((dXsub[:, :, :, :3], dXsub[:, :, :, -3:]), batch_size=batch_size, verbose=1) + + pulse_pred = yptest#[0] + pulse_pred = detrend(np.cumsum(pulse_pred), 100) + [b_pulse_pred, a_pulse_pred] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') + pulse_pred = scipy.signal.filtfilt(b_pulse_pred, a_pulse_pred, np.double(pulse_pred)) + pulse_pred = np.array(mms.fit_transform(pulse_pred.reshape(-1,1))).flatten() + + ##### ground truth data resampled ####### + print("\n", sample_data_path) + print(str(sample_data_path).find("COHFACE")) + print(str(sample_data_path).find("UBFC-PHYS")) + if(str(sample_data_path).find("COHFACE") > 0): + truth_path = args.video_path.replace(".avi", "_dataFile.hdf5") # akutell für COHACE... + elif(str(sample_data_path).find("UBFC-PHYS") > 0): + print("OK") + truth_path = args.video_path.replace("vid_", "") + "_dataFile.hdf5" + else: + return("Error in finding the ground truth signal...") + print(truth_path) + gound_truth_file = h5py.File(truth_path, "r") + pulse_truth = gound_truth_file["pulse"] ### range ground truth from 0 to 1 + pulse_truth = detrend(np.cumsum(pulse_truth), 100) + [b_pulse_tr, a_pulse_tr] = butter(1, [0.75 / fs * 2, 2.5 / fs * 2], btype='bandpass') + pulse_truth = scipy.signal.filtfilt(b_pulse_tr, a_pulse_tr, np.double(pulse_truth)) + pulse_truth = np.array(mms.fit_transform(pulse_truth.reshape(-1,1))).flatten() + pulse_truth = pulse_truth[0: dXsub_len] + + ########### Peaks ########### + working_data_pred, measures_pred = hp.process(pulse_pred, fs, calc_freq=True) + working_data_truth, measures_truth = hp.process(pulse_truth, fs, calc_freq=True) + peaks_pred = working_data_pred['peaklist'] + peaks_truth = working_data_truth['peaklist'] + + ######## x-axis: time ######### + duration_vid = dXsub_len/fs + x_axis = np.linspace(0, duration_vid, dXsub_len) + + ########## Plot ################## + plt.figure() #subplot(211) + plt.plot(x_axis, pulse_pred, label='Prediction', color ='#E6001A') + plt.plot(x_axis[peaks_truth], pulse_truth[peaks_truth], "x", color="#721085") + plt.plot(x_axis[peaks_pred], pulse_pred[peaks_pred], "x", color ='#E6001A') + plt.title('Pulse Prediction') + plt.xlabel("time (s)") + plt.ylabel("normalized Signal [a.u.]") + plt.plot(x_axis, pulse_truth, label='ground truth', color="#721085") + plt.legend() + + plt.figure() + plt.subplot(211) + plt.plot(x_axis, pulse_truth, label='Ground truth') + plt.plot(x_axis[peaks_truth], pulse_truth[peaks_truth], "x") + plt.ylabel("normalized Signal [a.u.]") + plt.xlabel("time (s)") + plt.title('Ground truth') + plt.subplot(212) + plt.plot(x_axis, pulse_pred, label='Prediction') + plt.plot(x_axis[peaks_pred], pulse_pred[peaks_pred], "x") + plt.title("Prediction") + plt.ylabel("normalized Signal [a.u.]") + plt.xlabel("time (s)") + plt.legend() + plt.show() + + ########### IBI ############# + ibi_truth = working_data_truth['RR_list_cor'] + print(ibi_truth) + ibi_pred = working_data_pred['RR_list_cor'] + print(ibi_pred) + ######### HRV featurs ############## + print("HRV Truth: ",measures_truth) + print("HRV Pred: ", measures_pred) + ####### Logging ############# + # neuer Ordner für Tests + file = open(str(sample_data_path).replace(".avi", "_result.txt"),"w") + file.write("LogFile\n\n") + file.write("IBI: "), file.write(str(ibi_pred)) + file.write("\nHR and HRVfeatures: "), file.write(str(measures_pred)) + + file.write("\n\n\nGround truth infos!") + file.write("\nHR and HRV features: "), file.write(str(measures_truth)) + +if __name__ == "__main__": + + parser = argparse.ArgumentParser() + parser.add_argument('--video_path', type=str, help='processed video path') + parser.add_argument('--batch_size', type=int, default = 100, help='batch size (multiplier of 10)') + parser.add_argument('--trained_model', type=str, default = './rPPG-checkpoints/testCohFace1/cv_0_epoch24_model.hdf5', help='path to trained model') + + args = parser.parse_args() + + predict_vitals(args) + + +#python code/predict_vitals_new.py --video_path "D:\Databases\1)Training\COHFACE\1\1\data.avi" --trained_model ./cv_0_epoch24_model.hdf5 +#./rPPG-checkpoints/testCohFace1/cv_0_epoch24_model.hdf5 +#./rPPG-checkpoints/test1/cv_0_epoch04_model.hdf5' +#python code/predict_vitals_new.py --video_path "D:\Databases\1)Training\UBFC-PHYS\s1\vid_s1_T1.avi" --trained_model ./cv_0_epoch24_model.hdf5 \ No newline at end of file diff --git a/code/predict_vitals_oneVideo.py b/code/predict_vitals_oneVideo.py new file mode 100644 index 0000000..e5328cc --- /dev/null +++ b/code/predict_vitals_oneVideo.py @@ -0,0 +1,180 @@ +''' +Script to predict a video over a defined model. +Subsequent saving of the model outputs and +analyses of the heartpy module. +input example: +python code/predict_vitals_oneVideo.py --video_path "path-to-video" + --trained_model "D:\Databases\4)Results\Version5\TS_CAN\cv_0_epoch24_model.hdf5" + --model_name "TS_CAN" +PPTS_CAN example: + --parameter "bpm,sdnn" +''' + +from aifc import Error +import argparse +import os +import numpy as np +import scipy.io +from model import CAN, CAN_3D, PPTS_CAN, PTS_CAN, TS_CAN, Hybrid_CAN +import matplotlib.pyplot as plt +from scipy.signal import butter +from inference_preprocess import preprocess_raw_frames, preprocess_raw_video, detrend +from sklearn.preprocessing import MinMaxScaler +import heartpy as hp + +def prepare_3D_CAN(dXsub): + frame_depth = 10 + num_window = int(dXsub.shape[0]) - frame_depth + 1 + tempX = np.array([dXsub[f:f + frame_depth, :, :, :] # (491, 10, 36, 36 ,6) (169, 10, 36, 36, 6) + for f in range(num_window)]) + tempX = np.swapaxes(tempX, 1, 3) # (169, 36, 36, 10, 6) + tempX = np.swapaxes(tempX, 1, 2) # (169, 36, 36, 10, 6) + return tempX + +def prepare_Hybrid_CAN(dXsub): + frame_depth = 10 + num_window = int(dXsub.shape[0]) - frame_depth + 1 + tempX = np.array([dXsub[f:f + frame_depth, :, :, :] # (169, 10, 36, 36, 6) + for f in range(num_window)]) + tempX = np.swapaxes(tempX, 1, 3) # (169, 36, 36, 10, 6) + tempX = np.swapaxes(tempX, 1, 2) # (169, 36, 36, 10, 6) + motion_data = tempX[:, :, :, :, :3] + apperance_data = np.average(tempX[:, :, :, :, -3:], axis=-2) + return motion_data, apperance_data + +def predict_vitals(args): + model_checkpoint = args.trained_model + sample_data_path = args.video_path + model_name = args.model_name + save_dir = args.save_dir + parameter_PPTS = str(args.parameter).split(",") + mms = MinMaxScaler() + img_rows = 36 + img_cols = 36 + frame_depth = 10 + # neuer Ordner für Tests + os.chdir(save_dir) + try: + os.makedirs(str(model_name)) + except: + print("Directory exists...") + save_dir = os.path.join(save_dir, str(model_name)) + + print("path: ",sample_data_path) + if sample_data_path[-4:] == ".avi": + dXsub, fps = preprocess_raw_video(sample_data_path, dim=36) + elif sample_data_path[-4:] == ".mp4": + dXsub, fps = preprocess_raw_video(sample_data_path, dim=36) + else: + dXsub, fps = preprocess_raw_frames(sample_data_path, dim=36) + print('dXsub shape', dXsub.shape, "fps: ", fps) + + if model_name == "PPTS_CAN": + dXsub_len = (dXsub.shape[0] // (frame_depth*10)) * (frame_depth*10) + dXsub = dXsub[:dXsub_len, :, :, :] + else: + dXsub_len = (dXsub.shape[0] // frame_depth) * frame_depth + dXsub = dXsub[:dXsub_len, :, :, :] + + if model_name == "TS_CAN": + model = TS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3)) + elif model_name == "3D_CAN": + model = CAN_3D(frame_depth, 32, 64, (img_rows, img_cols, frame_depth, 3)) + dXsub = prepare_3D_CAN(dXsub) + dXsub_len = (dXsub.shape[0] // (frame_depth)) * (frame_depth) + dXsub = dXsub[:dXsub_len, :, :, :,:] + elif model_name == "CAN": + model = CAN(32, 64, (img_rows, img_cols, 3)) + elif model_name == "Hybrid_CAN": + model = Hybrid_CAN(frame_depth, 32, 64, (img_rows, img_cols, frame_depth, 3), + (img_rows, img_cols, 3)) + dXsub1, dXsub2 = prepare_Hybrid_CAN(dXsub) + dXsub_len1 = (dXsub1.shape[0] // (frame_depth)) * (frame_depth) + dXsub1 = dXsub1[:dXsub_len1, :, :, :,:] + dXsub_len2 = (dXsub2.shape[0] // (frame_depth)) * (frame_depth) + dXsub2 = dXsub2[:dXsub_len2, :, :, :] + elif model_name == "PTS_CAN": + model = PTS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3)) + elif model_name == "PPTS_CAN": + model = PPTS_CAN(frame_depth, 32, 64, (img_rows, img_cols, 3), parameter=parameter_PPTS) + else: + raise NotImplementedError + + + model.load_weights(model_checkpoint) + if model_name == "3D_CAN": + yptest = model.predict((dXsub[:, :, :,: , :3], dXsub[:, :, :, : , -3:])) + elif model_name == "Hybrid_CAN": + yptest = model.predict((dXsub1, dXsub2)) + else: + yptest = model((dXsub[:, :, :, :3], dXsub[:, :, :, -3:]), training=False) + if model_name == "3D_CAN" or model_name == "Hybrid_CAN": + pulse_pred = yptest[:,0] + elif model_name != "PTS_CAN" and model_name != "PPTS_CAN": + pulse_pred = yptest + + else: + pulse_pred = yptest[0] + + if model_name == "PPTS_CAN": + parameter_out = yptest[2] + + pulse_pred = detrend(np.cumsum(pulse_pred), 100) + [b_pulse_pred, a_pulse_pred] = butter(1, [0.75 / fps * 2, 2.5 / fps * 2], btype='bandpass') + pulse_pred = scipy.signal.filtfilt(b_pulse_pred, a_pulse_pred, np.double(pulse_pred)) + pulse_pred = np.array(mms.fit_transform(pulse_pred.reshape(-1,1))).flatten() + + ########### Peaks ########### + working_data_pred, measures_pred = hp.process(pulse_pred, fps, calc_freq=True) + peaks_pred = working_data_pred['peaklist'] + + ######## x-axis: time ######### + duration_vid = dXsub_len/fps + x_axis = np.linspace(0, duration_vid, dXsub_len) + + ########## Plot ################## + peaks_pred_new = [] + for peak in peaks_pred: + if (peak > 400 and peak <700): + peaks_pred_new.append(peak-400) + + path_plot = save_dir + "/plot.png" + print(path_plot) + + plt.figure() #subplot(211) + plt.plot(x_axis, pulse_pred, "#E6001A", label='rPPG signal') + plt.plot(x_axis[peaks_pred], pulse_pred[peaks_pred], "x", color ='#E6001A') + plt.title('rPPG signal') + plt.ylabel("normalized Signal [a.u.]") + plt.xlabel("time (s)") + plt.legend() + plt.savefig(path_plot) + + file_rPPG = open(str(save_dir) + "/rPPG_out.txt","w") + for value in pulse_pred: + file_rPPG.write(str(value))# + file_rPPG.write("\n") + file_rPPG.close() + if model_name =="PPTS_CAN": + file_parameter = open(str(save_dir) + "/parameter_out.txt","w") + file_parameter.write(str(parameter_out)) + file_parameter.close() + file_hrAnalysis = open(str(save_dir) + "/HRVAnalysis.txt","w") + file_hrAnalysis.write(str(measures_pred)) + file_hrAnalysis.close() + + + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument('--video_path', type=str, help='processed video path') + parser.add_argument('--save_dir', type=str, help='save dir path') + parser.add_argument('--trained_model', type=str, default = "D:/Databases/4)Results/Version4/TS_Databases/TS_CAN_COHFACE_2GPU/cv_0_epoch24_model.hdf5", help='path to trained model') + parser.add_argument('--model_name', type=str, help='name of model (TS_CAN, PTS_CAN,...)') + parser.add_argument('--parameter', type=str, help='parameter for PPTS_CAN ("bpm,sdnn")') + + args = parser.parse_args() + + predict_vitals(args) + +#python code/predict_vitals_oneVideo.py --video_path "C:\Users\sarah\OneDrive\Desktop\UBFC\DATASET_2\subject4\vid.avi" --trained_model "D:\Databases\4)Results\Version5\TS_CAN\cv_0_epoch24_model.hdf5" --model_name "TS_CAN" --save_dir "D:\Databases\Test" \ No newline at end of file diff --git a/code/prepare_databases.py b/code/prepare_databases.py new file mode 100644 index 0000000..ae7bd6e --- /dev/null +++ b/code/prepare_databases.py @@ -0,0 +1,248 @@ +from pre_process import sort_video_list_, split_subj_ +from inference_preprocess import preprocess_raw_video + +import h5py +import itertools +import numpy as np +from scipy import signal +import os +from sklearn.preprocessing import MinMaxScaler +import heartpy as hp +import pandas as pd + +def hr_analysis(path, hr_discrete, frames_vid, fps_vid, nr=None): + mms = MinMaxScaler() + mean = hr_discrete.mean() + std = np.std(hr_discrete) + + upper_limit = mean + std*3 + lower_limit = mean - std*3 + + for x in range(0, len(hr_discrete)): + if hr_discrete[x] > upper_limit: + hr_discrete[x] = upper_limit + elif hr_discrete[x] < lower_limit: + hr_discrete[x] = lower_limit + + hrdata_res = np.array(signal.resample(hr_discrete, frames_vid)) + hrdata_res = np.array(mms.fit_transform(hrdata_res.reshape(-1,1))).flatten() # normalization + + working_data, measures = hp.process(hrdata_res, fps_vid, calc_freq=True) + + plot = hp.plotter(working_data, measures, show=False, title = 'Heart Rate Signal and Peak Detection') + if nr ==None: + path = path.replace('vid', 'plot_truthData').replace('.avi', '.jpg').replace('data', 'plot_truthData') + else: + path = path.replace('vid', 'plot_truthData_' + str(nr)).replace('.avi', '.jpg').replace('data', 'plot_truthData') + plot.savefig(path) + + return working_data, measures + +def dataSet_preprocess(vid, name): + if name== "COHFACE": + if os.path.exists(str(vid).replace(".avi", "_vid.hdf5")): + os.remove(str(vid).replace(".avi", "_vid.hdf5")) + print("deleted") + if os.path.exists(str(vid).replace(".avi", "_dataFile.hdf5")): + os.remove(str(vid).replace(".avi", "_dataFile.hdf5")) + print("deleted") + + dXsub, fps = preprocess_raw_video(vid, 36) + print(dXsub.shape, " fps: ", fps) + nframesPerVideo = dXsub.shape[0] + + # ground truth data: + truth_data_path = str(vid).replace(".avi", ".hdf5") + hf = h5py.File(truth_data_path, 'r') + pulse = np.array(hf['pulse']) + + hf.close() # close the hdf5 file + + return nframesPerVideo, fps, dXsub, pulse + + elif name=="UBFC_PHYS": + if os.path.exists(str(vid).replace(".avi", "_data.hdf5")): + os.remove(str(vid).replace(".avi", "_data.hdf5")) + print("deleted") + if os.path.exists(str(vid).replace(".avi", "_dataFile.hdf5").replace('vid_', '')): + os.rename(str(vid).replace(".avi", "_dataFile.hdf5").replace('vid_', ''), str(vid).replace(".avi", "_dataFile.hdf5").replace('vid_', 'ALL_')) + print("renamed") + if os.path.exists(str(vid).replace(".avi", "_dataFileAll.hdf5").replace('vid_', '')): + os.remove(str(vid).replace(".avi", "_dataFileAll.hdf5").replace('vid_', '')) + print("deleted") + + dXsub, fps = preprocess_raw_video(vid, 36) + print(dXsub.shape, " fps: ", fps) + nframesPerVideo = dXsub.shape[0] + + # ground truth data: + truth_data_path = str(vid).replace(".avi", ".csv").replace('vid', 'bvp') + hf = pd.read_csv(truth_data_path) + pulse = np.array(hf) + + return nframesPerVideo, fps, dXsub, pulse + if name == "UBFC": + if os.path.exists(str(vid).replace(".avi", "_vid.hdf5")): + os.remove(str(vid).replace(".avi", "_vid.hdf5")) + print("deleted") + if os.path.exists(str(vid).replace(".avi", "_dataFile.hdf5")): + os.remove(str(vid).replace(".avi", "_dataFile.hdf5")) + print("deleted") + + dXsub, fps = preprocess_raw_video(vid, 36) + print(dXsub.shape, " fps: ", fps) + nframesPerVideo = dXsub.shape[0] + + # ground truth data: + truth_data_path = str(vid).replace(".avi", ".txt").replace('vid', 'ground_truth') + data = open(truth_data_path, 'r').read() + data = str(data).split(" ") + data = data[1:] + pulse = np.array(list(map(float, data[0: nframesPerVideo]))) + + return nframesPerVideo, fps, dXsub, pulse + +def process_save(nframesPerVideo, fps, dXsub, pulse, vid): + # HR and HRV analysis + working_data, measures = hr_analysis(vid, pulse, nframesPerVideo, fps) + + # Data for H5PY + pulse_res = working_data['hr'] # resampled and normalized HR + nn_list = working_data['RR_list_cor'] # nn-intervals + parameter = str(measures) # HR and HRV Parameter + + peak_list = working_data['peaklist'] + if not isinstance(peak_list, list): + peak_list = peak_list.tolist() + removed = working_data['removed_beats'] + for item in removed: + peak_list.remove(item) # list with position of the peaks + + ##### save data ###### + newPath_name = str(vid).replace(".avi", "_dataFile.hdf5") + if (str(vid).find("UBFC") >=0): + newPath_name = newPath_name.replace("vid_", "") + data_file = h5py.File(newPath_name, 'a') + data_file.create_dataset('data', data=dXsub) # write the data to hdf5 file + data_file.create_dataset('pulse', data=pulse_res) + data_file.create_dataset('peaklist', data=peak_list) + data_file.create_dataset('nn', data=nn_list) + data_file.create_dataset('parameter', data=parameter) + data_file.close() + +def process_save_UBFC(nframesPerVideo, fps, dXsub, pulse, vid): + ##### split in 3 one minute parts ########### + frames_per_dataPacket = int(fps*60) + frames_per_pulsePacket = int(64*60) + dXsub_1 = dXsub[0:frames_per_dataPacket,:,:,:] + dXsub_2 = dXsub[frames_per_dataPacket: frames_per_dataPacket*2,:,:,:] + dXsub_3 = dXsub[frames_per_dataPacket*2:,:,:,:] + pulse_1 = pulse[0:frames_per_pulsePacket,:] + pulse_2 = pulse[frames_per_pulsePacket:frames_per_pulsePacket*2,:] + pulse_3 = pulse[frames_per_pulsePacket*2:,:] + # HR and HRV analysis + working_data1, measures1 = hr_analysis(vid, pulse_1, dXsub_1.shape[0], fps, nr=1) + working_data2, measures2 = hr_analysis(vid, pulse_2, dXsub_2.shape[0], fps, nr=2) + working_data3, measures3 = hr_analysis(vid, pulse_3, dXsub_3.shape[0], fps, nr=3) + + # Data for H5PY + pulse_res1 = working_data1['hr'] # resampled and normalized HR + pulse_res2 = working_data2['hr'] # resampled and normalized HR + pulse_res3 = working_data3['hr'] # resampled and normalized HR + nn_list1 = working_data1['RR_list_cor'] # nn-intervals + nn_list2 = working_data2['RR_list_cor'] # nn-intervals + nn_list3 = working_data3['RR_list_cor'] # nn-intervals + parameter1 = str(measures1) # HR and HRV Parameter + parameter2 = str(measures2) # HR and HRV Parameter + parameter3 = str(measures3) # HR and HRV Parameter + + peak_list1 = working_data1['peaklist'] + if not isinstance(peak_list1, list): + peak_list1 = peak_list1.tolist() + removed1 = working_data1['removed_beats'] + for item in removed1: + peak_list1.remove(item) # list with position of the peaks + + peak_list2 = working_data2['peaklist'] + if not isinstance(peak_list2, list): + peak_list2 = peak_list2.tolist() + removed2 = working_data2['removed_beats'] + for item in removed2: + peak_list2.remove(item) # list with position of the peaks + + peak_list3 = working_data3['peaklist'] + if not isinstance(peak_list3, list): + peak_list3 = peak_list3.tolist() + removed3 = working_data3['removed_beats'] + for item in removed3: + peak_list3.remove(item) # list with position of the peaks + + + newPath_name1 = str(vid).replace(".avi", "_0_dataFile.hdf5").replace('vid_', '') + newPath_name2 = str(vid).replace(".avi", "_1_dataFile.hdf5").replace('vid_', '') + newPath_name3 = str(vid).replace(".avi", "_2_dataFile.hdf5").replace('vid_', '') + data_file1 = h5py.File(newPath_name1, 'a') + data_file1.create_dataset('data', data=dXsub_1) # write the data to hdf5 file + data_file1.create_dataset('pulse', data=pulse_res1) + data_file1.create_dataset('peaklist', data=peak_list1) + data_file1.create_dataset('nn', data=nn_list1) + data_file1.create_dataset('parameter', data=parameter1) + data_file1.close() + data_file2 = h5py.File(newPath_name2, 'a') + data_file2.create_dataset('data', data=dXsub_2) # write the data to hdf5 file + data_file2.create_dataset('pulse', data=pulse_res2) + data_file2.create_dataset('peaklist', data=peak_list2) + data_file2.create_dataset('nn', data=nn_list2) + data_file2.create_dataset('parameter', data=parameter2) + data_file2.close() + data_file3 = h5py.File(newPath_name3, 'a') + data_file3.create_dataset('data', data=dXsub_3) # write the data to hdf5 file + data_file3.create_dataset('pulse', data=pulse_res3) + data_file3.create_dataset('peaklist', data=peak_list3) + data_file3.create_dataset('nn', data=nn_list3) + data_file3.create_dataset('parameter', data=parameter3) + data_file3.close() + +def build_h5py(vid, name): + print("Dataset: ", name) + print("Current: ", vid) + nframesPerVideo, fps, dXsub, pulse = dataSet_preprocess(vid, name) + if name != "UBFC_PHYS": + process_save(nframesPerVideo, fps, dXsub, pulse, vid) + else: + process_save_UBFC(nframesPerVideo, fps, dXsub, pulse, vid) + + print("next") + + +def prepare_database(name, tasks, data_dir): + if name == "UBFC_PHYS": + taskList = list(range(1, tasks+1)) + elif name == "COHFACE": + taskList = list(range(0, tasks)) + elif name == "UBFC": + taskList = [0] + else: + print("Not implemented yet") + subTrain, subTest = split_subj_(data_dir, name) + print("subTrain: ", subTrain) + print("subTest: ", subTest) + video_path_list_tr = sort_video_list_(data_dir, taskList, subTrain, name, True) + video_path_list_test = sort_video_list_(data_dir, taskList, subTest, name, False) + video_path_list_tr = list(itertools.chain(*video_path_list_tr)) + video_path_list_test = list(itertools.chain(*video_path_list_test)) + + for vid in video_path_list_tr: + build_h5py(vid, name) + for vid in video_path_list_test: + build_h5py(vid, name) + + + +#data_dir = "D:/Databases" +#prepare_database("COHFACE", 4, data_dir) +data_dir = "/mnt/share/StudiShare/sarah/Databases/" + +data_dir = "C:/Users/sarah/OneDrive/Desktop/UBFC/DATASET_2" +prepare_database("UBFC", 1, data_dir) + diff --git a/code/train.py b/code/train.py index e1c6b33..979ea5c 100755 --- a/code/train.py +++ b/code/train.py @@ -1,27 +1,37 @@ ''' -Training Script for Multi-Task Temporal Shift Attention Networks for On-Device Contactless Vitals Measurement +Training Script for Multi-Task Temporal Shift Attention Networks for On-Device +Contactless Vitals Measurement Author: Xin Liu, Daniel McDuff + +Further development: Sarah Quehl ''' # %% from __future__ import print_function import argparse -import itertools import json import os +os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2' +from xmlrpc.client import boolean +from losses import negPearsonLoss, gaussian_loss, MAPE_parameter_loss, time_error_loss import numpy as np import scipy.io import tensorflow as tf - +from tensorflow.python.keras.optimizers import adadelta_v2 from data_generator import DataGenerator from model import HeartBeat, CAN, CAN_3D, Hybrid_CAN, TS_CAN, MTTS_CAN, \ - MT_Hybrid_CAN, MT_CAN_3D, MT_CAN -from pre_process import get_nframe_video, split_subj, sort_video_list + MT_Hybrid_CAN, MT_CAN_3D, MT_CAN, PTS_CAN, PPTS_CAN +from pre_process import split_subj_, sort_dataFile_list_, collect_subj np.random.seed(100) # for reproducibility -tf.test.is_gpu_available() +print("START!") +list_gpu = tf.config.list_physical_devices('GPU') +#tf.config.experimental.set_memory_growth(list_gpu[0], enable=True) +#tf.config.experimental.set_memory_growth(list_gpu[1], enable=True) +print(list_gpu) tf.keras.backend.clear_session() +tf.autograph.set_verbosity(10) print(tf.__version__) # %% @@ -30,7 +40,7 @@ parser.add_argument('-exp', '--exp_name', type=str, help='experiment name') parser.add_argument('-i', '--data_dir', type=str, help='Location for the dataset') -parser.add_argument('-o', '--save_dir', type=str, default='./rPPG-checkpoints', +parser.add_argument('-o', '--save_dir', type=str, default='/home/quehl/Results/', help='Location for parameter checkpoints and samples') parser.add_argument('-a', '--nb_filters1', type=int, default=32, help='number of convolutional filters to use') @@ -41,7 +51,7 @@ parser.add_argument('-d', '--dropout_rate2', type=float, default=0.5, help='dropout rates') parser.add_argument('-l', '--lr', type=float, default=1.0, - help='learning rate') + help='learning rate') parser.add_argument('-e', '--nb_dense', type=int, default=128, help='number of dense units') parser.add_argument('-f', '--cv_split', type=int, default=0, @@ -52,23 +62,24 @@ help='nb_task') parser.add_argument('-fd', '--frame_depth', type=int, default=10, help='frame_depth for CAN_3D, TS_CAN, Hybrid_CAN') -parser.add_argument('-temp', '--temporal', type=str, default='MTTS_CAN', +parser.add_argument('-temp', '--temporal', type=str, default='PTS_CAN', help='CAN, MT_CAN, CAN_3D, MT_CAN_3D, Hybrid_CAN, \ - MT_Hybrid_CAN, TS_CAN, MTTS_CAN ') + MT_Hybrid_CAN, TS_CAN, MTTS_CAN. PTS_CAN ') parser.add_argument('-save', '--save_all', type=int, default=1, help='save all or not') parser.add_argument('-resp', '--respiration', type=int, default=0, help='train with resp or not') +parser.add_argument('-database', '--database_name', type=str, + default="MIX2", help='Which database') +parser.add_argument('-lf1', '--loss_function1', type=str, default="MSE", help="MSE,NegPea") +parser.add_argument('-lf2', '--loss_function2', type=str, default="MSE", help="MSE,NegPea, Gauss_Peak, time_Error") +parser.add_argument('-min', '--decrease_database', type=boolean, default=False) +parser.add_argument('-ml', '--maxFrames_video', type=int, default=2050, help="frames") +parser.add_argument('-p', '--parameter', default=None, help="bpm, sdnn, pnn50, lfhf") args = parser.parse_args() print('input args:\n', json.dumps(vars(args), indent=4, separators=(',', ':'))) # pretty print args -# %% Spliting Data - -print('Spliting Data...') -subNum = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 20, 21, 22, 23, 25, 26, 27]) -taskList = list(range(1, args.nb_task+1)) - # %% Training @@ -79,39 +90,49 @@ def train(args, subTrain, subTest, cv_split, img_rows=36, img_cols=36): print('subTest', subTest) input_shape = (img_rows, img_cols, 3) + maxLen_video = args.maxFrames_video - path_of_video_tr = sort_video_list(args.data_dir, taskList, subTrain) - path_of_video_test = sort_video_list(args.data_dir, taskList, subTest) - path_of_video_tr = list(itertools.chain(*path_of_video_tr)) # Fllaten the list - path_of_video_test = list(itertools.chain(*path_of_video_test)) + path_of_video_tr = sort_dataFile_list_(args.data_dir, subTrain, args.database_name, trainMode=True) + path_of_video_test = sort_dataFile_list_(args.data_dir, subTest, args.database_name, trainMode=False) - print('sample path: ', path_of_video_tr[0]) - nframe_per_video = get_nframe_video(path_of_video_tr[0]) - print('Trian Length: ', len(path_of_video_tr)) + #nframe_per_video = get_nframe_video_(path_of_video_tr[0]) + print('Train Length: ', len(path_of_video_tr)) print('Test Length: ', len(path_of_video_test)) - print('nframe_per_video', nframe_per_video) - - strategy = tf.distribute.MirroredStrategy() + if len(list_gpu) > 1: + print("Using MultiWorkerMirroredStrategy") + strategy = tf.distribute.MultiWorkerMirroredStrategy() + else: + print("Using MirroredStrategy") + strategy = tf.distribute.MirroredStrategy() print('Number of devices: {}'.format(strategy.num_replicas_in_sync)) with strategy.scope(): - if strategy.num_replicas_in_sync == 4: - print("Using 4 GPUs for training") - if args.temporal == 'CAN' or args.temporal == 'MT_CAN': - args.batch_size = 32 - elif args.temporal == 'CAN_3D' or args.temporal == 'MT_CAN_3D': - args.batch_size = 12 - elif args.temporal == 'TS_CAN' or args.temporal == 'MTTS_CAN': - args.batch_size = 32 - elif args.temporal == 'Hybrid_CAN' or args.temporal == 'MT_Hybrid_CAN': - args.batch_size = 16 - else: - raise ValueError('Unsupported Model Type!') - elif strategy.num_replicas_in_sync == 8: + + if args.temporal == 'CAN' or args.temporal == 'MT_CAN': + args.batch_size = 16 + elif args.temporal == 'CAN_3D' or args.temporal == 'MT_CAN_3D': + args.batch_size = 2 + elif args.temporal == 'TS_CAN' or args.temporal == 'MTTS_CAN'\ + or args.temporal == 'PTS_CAN': + args.batch_size = 12#32 + elif args.temporal == 'PPTS_CAN': + args.batch_size = 2 + elif args.temporal == 'Hybrid_CAN' or args.temporal == 'MT_Hybrid_CAN': + args.batch_size = 2# 16 + else: + raise ValueError('Unsupported Model Type!') + + if strategy.num_replicas_in_sync == 8: print('Using 8 GPUs for training!') args.batch_size = args.batch_size * 2 elif strategy.num_replicas_in_sync == 2: + print('Using 2 GPUs for training!') args.batch_size = args.batch_size // 2 + elif strategy.num_replicas_in_sync == 1: + print('Using 1 GPU for training!') + args.batch_size = 1#4 + elif strategy.num_replicas_in_sync == 4: + print("Using 4 GPUs for training") else: raise Exception('Only supporting 4 GPUs or 8 GPUs now. Please adjust learning rate in the training script!') @@ -139,6 +160,18 @@ def train(args, subTrain, subTest, cv_split, img_rows=36, img_cols=36): input_shape = (img_rows, img_cols, 3) model = TS_CAN(args.frame_depth, args.nb_filters1, args.nb_filters2, input_shape, dropout_rate1=args.dropout_rate1, dropout_rate2=args.dropout_rate2, nb_dense=args.nb_dense) + elif args.temporal == 'PTS_CAN': + print('Using PTS_CAN: with PeakLocation!') + input_shape = (img_rows, img_cols, 3) + model = PTS_CAN(args.frame_depth, args.nb_filters1, args.nb_filters2, input_shape, + dropout_rate1=args.dropout_rate1, dropout_rate2=args.dropout_rate2, nb_dense=args.nb_dense) + elif args.temporal == 'PPTS_CAN': + print('Using PPTS_CAN: with PeakLocation!') + input_shape = (img_rows, img_cols, 3) + args.parameter = str(args.parameter).split(",") + print(args.parameter) + model = PPTS_CAN(args.frame_depth, args.nb_filters1, args.nb_filters2, input_shape, + dropout_rate1=args.dropout_rate1, dropout_rate2=args.dropout_rate2, nb_dense=args.nb_dense, parameter=args.parameter) elif args.temporal == 'MTTS_CAN': print('Using MTTS_CAN!') input_shape = (img_rows, img_cols, 3) @@ -162,24 +195,125 @@ def train(args, subTrain, subTest, cv_split, img_rows=36, img_cols=36): nb_dense=args.nb_dense) else: raise ValueError('Unsupported Model Type!') + + model.summary() - optimizer = tf.keras.optimizers.Adadelta(learning_rate=args.lr) + optimizer = adadelta_v2.Adadelta(learning_rate=args.lr) if args.temporal == 'MTTS_CAN' or args.temporal == 'MT_Hybrid_CAN' or args.temporal == 'MT_CAN_3D' or \ args.temporal == 'MT_CAN': losses = {"output_1": "mean_squared_error", "output_2": "mean_squared_error"} loss_weights = {"output_1": 1.0, "output_2": 1.0} model.compile(loss=losses, loss_weights=loss_weights, optimizer=optimizer) + + elif args.temporal == 'PTS_CAN': + # output 1: rPPG Signal + if args.loss_function1 == "MSE": + loss1 = 'mean_squared_error' + loss_weights1 = 1 + elif args.loss_function1 == "NegPea": + loss1 = negPearsonLoss + loss_weights1 = 1 + elif args.loss_function1 == "MSE_negPea": + loss1a = 'mean_squared_error' + loss1b = negPearsonLoss + loss1 = [loss1a, loss1b] + loss_weights1 = [1,1] + elif args.loss_function1 == "Gauss_Peak": + raise NotImplementedError + # output 2: Gaussdistribution around peak locations or TimeError + if args.loss_function2 == "MSE": + loss2 = 'mean_squared_error' + loss_weights2 = 1 + elif args.loss_function2 == "NegPea": + loss2 = negPearsonLoss + loss_weights2 = 1 + raise NotImplementedError + elif args.loss_function2 == "Gauss_Peak": + loss2 = gaussian_loss + loss_weights2 = 1 + elif args.loss_function2 == "time_Error": + loss2 = time_error_loss + loss_weights2 = 1 + + losses = {"output_1": loss1, "output_2": loss2} + loss_weights = {"output_1": loss_weights1, "output_2": loss_weights2} + model.compile(loss=losses, loss_weights=loss_weights, optimizer=optimizer) + + elif args.temporal == 'PPTS_CAN': + # output 1: rPPG Signal + if args.loss_function1 == "MSE": + loss1 = 'mean_squared_error' + loss_weights1 = 1 + elif args.loss_function1 == "NegPea": + loss1 = negPearsonLoss + loss_weights1 = 1 + elif args.loss_function1 == "MSE_negPea": + loss1a = 'mean_squared_error' + loss1b = negPearsonLoss + loss1 = [loss1a, loss1b] + loss_weights1 = [1,1] + elif args.loss_function1 == "Gauss_Peak": + raise NotImplementedError + # output 2: Gaussdistribution around peak locations or TimeError + if args.loss_function2 == "MSE": + loss2 = 'mean_squared_error' + loss_weights2 = 1 + elif args.loss_function2 == "NegPea": + loss2 = negPearsonLoss + loss_weights2 = 1 + raise NotImplementedError + elif args.loss_function2 == "Gauss_Peak": + loss2 = gaussian_loss + loss_weights2 = 1 + elif args.loss_function2 == "time_Error": + loss2 = time_error_loss + loss_weights2 = 1 + else: + raise NotImplementedError + # output 3: different Parameter + loss3 = MAPE_parameter_loss + loss_weights3 = 1 + + losses = {"output_1": loss1, "output_3": loss3}#"output_2": loss2 + loss_weights = {"output_1": loss_weights1, "output_3": loss_weights3}# "output_2": loss_weights2, + model.compile(loss=losses, loss_weights=loss_weights, optimizer=optimizer) + else: - model.compile(loss='mean_squared_error', optimizer=optimizer) + if args.loss_function1 == "MSE": + model.compile(loss='mean_squared_error', optimizer=optimizer) + elif args.loss_function1 == "negPea": + print("negative Pearson Loss ") + loss = negPearsonLoss + model.compile(loss=loss, optimizer=optimizer) + elif args.loss_function1 == "MSE_negPea": + loss1 = 'mean_squared_error' + loss2 = negPearsonLoss + losses = [loss1, loss2] + loss_weights = [1,1] + model.compile(loss= losses, loss_weights=loss_weights, optimizer=optimizer) + else: + return ValueError('Unsupported Loss Function') + print('learning rate: ', args.lr) + print('batch size: ', args.batch_size) + + if args.loss_function2 == "time_Error": + timeError = True + else: + timeError = False # %% Create data genener - training_generator = DataGenerator(path_of_video_tr, nframe_per_video, (img_rows, img_cols), + training_generator = DataGenerator(path_of_video_tr, maxLen_video, (img_rows, img_cols), batch_size=args.batch_size, frame_depth=args.frame_depth, - temporal=args.temporal, respiration=args.respiration) - validation_generator = DataGenerator(path_of_video_test, nframe_per_video, (img_rows, img_cols), + temporal=args.temporal, respiration=args.respiration, + database_name=args.database_name, time_error_loss=timeError, + truth_parameter=args.parameter) + validation_generator = DataGenerator(path_of_video_test, maxLen_video, (img_rows, img_cols), batch_size=args.batch_size, frame_depth=args.frame_depth, - temporal=args.temporal, respiration=args.respiration) + temporal=args.temporal, respiration=args.respiration, + database_name=args.database_name, time_error_loss=timeError, + truth_parameter=args.parameter) + # %% Checkpoint Folders checkpoint_folder = str(os.path.join(args.save_dir, args.exp_name)) if not os.path.exists(checkpoint_folder): @@ -196,10 +330,10 @@ def train(args, subTrain, subTest, cv_split, img_rows=36, img_cols=36): save_best_only=False, verbose=1) csv_logger = tf.keras.callbacks.CSVLogger(filename=cv_split_path + '_train_loss_log.csv') hb_callback = HeartBeat(training_generator, validation_generator, args, str(cv_split), checkpoint_folder) - + # %% Model Training and Saving Results - history = model.fit(x=training_generator, validation_data=validation_generator, epochs=args.nb_epoch, verbose=1, - shuffle=False, callbacks=[csv_logger, save_best_callback, hb_callback], validation_freq=4) + history = model.fit(x=training_generator, validation_data=validation_generator, epochs=args.nb_epoch, + verbose=1, shuffle=True, callbacks=[csv_logger, save_best_callback, hb_callback], validation_freq=4) val_loss_history = history.history['val_loss'] val_loss = np.array(val_loss_history) @@ -218,13 +352,17 @@ def train(args, subTrain, subTest, cv_split, img_rows=36, img_cols=36): print('****************************************') print('Start saving predicitions from the last epoch') - training_generator = DataGenerator(path_of_video_tr, nframe_per_video, (img_rows, img_cols), + training_generator = DataGenerator(path_of_video_tr, maxLen_video, (img_rows, img_cols), batch_size=args.batch_size, frame_depth=args.frame_depth, - temporal=args.temporal, respiration=args.respiration, shuffle=False) + temporal=args.temporal, respiration=args.respiration, shuffle=False, + database_name=args.database_name, time_error_loss=timeError, + truth_parameter=args.parameter) - validation_generator = DataGenerator(path_of_video_test, nframe_per_video, (img_rows, img_cols), + validation_generator = DataGenerator(path_of_video_test, maxLen_video, (img_rows, img_cols), batch_size=args.batch_size, frame_depth=args.frame_depth, - temporal=args.temporal, respiration=args.respiration, shuffle=False) + temporal=args.temporal, respiration=args.respiration, shuffle=False, + database_name=args.database_name, time_error_loss=timeError, + truth_parameter=args.parameter) yptrain = model.predict(training_generator, verbose=1) scipy.io.savemat(checkpoint_folder + '/yptrain_best_' + '_cv' + str(cv_split) + '.mat', @@ -233,11 +371,47 @@ def train(args, subTrain, subTest, cv_split, img_rows=36, img_cols=36): scipy.io.savemat(checkpoint_folder + '/yptest_best_' + '_cv' + str(cv_split) + '.mat', mdict={'yptest': yptest}) + file = open(checkpoint_folder + "/log.txt","w") + file.write("LogFile\n\n") + file.write("Name: "), file.write(args.exp_name) + file.write("\nModel: "), file.write(args.temporal) + file.write("\nBatch Size: "), file.write(str(args.batch_size)) + file.write("\nLoss Function (output1): "), file.write(args.loss_function1) + file.write("\nLoss Function (output2): "), file.write(args.loss_function2) + file.write("\nLoss Function (output3): "), file.write("MAPE") + file.write("\nMax Frames Video: "), file.write(str(args.maxFrames_video)) + file.write("\nLearningrate: "), file.write(str(args.lr)) + file.write("\nTrain Subjects: "), file.write(str(subTrain)) + file.write("\nValidation Subjects: "), file.write(str(subTest)) + file.close() + print('Finish saving the results from the last epoch') # %% Training print('Using Split ', str(args.cv_split)) -subTrain, subTest = split_subj(args.data_dir, args.cv_split, subNum) +print("DatabaseName: ", args.database_name) +# Mix1: COHFACE and UBFC-Phys +# Mix2: COHFACE and UBFC-rPPG +if args.database_name != "MIX1" and args.database_name != "MIX2": + subTrain, subTest = split_subj_(args.data_dir, args.database_name) +else: + subTrain, subTest = collect_subj(args.data_dir, args.database_name) + +if args.decrease_database == True: + if args.database_name == "COHFACE": + subTrain = subTrain[0:10] + subTest = subTest[0:3] + elif args.database_name == "UBFC_PHYS": + subTrain = subTrain[0:25] + subTest = subTest[0:10] + elif args.database_name == "MIX1": + for key in subTrain.keys(): + subTrain[key] = subTrain[key][0:6] + for key in subTest.keys(): + subTest[key] = subTest[key][0:3] + train(args, subTrain, subTest, args.cv_split) + + diff --git a/cv_0_epoch24_model.hdf5 b/cv_0_epoch24_model.hdf5 new file mode 100644 index 0000000..c5be08b Binary files /dev/null and b/cv_0_epoch24_model.hdf5 differ