11858895593

Committed 15 Nov 2024 03:23PM UTC coverage: 0.783% (-0.006%) from 0.789%

Build # 11858895593

Build Type

push

github

Committed by

Pablo1990

Commit Message

Wound indentation and create video

Run Details

0 of 1025 branches covered (0.0%)

Branch coverage included in aggregate %.

0 of 39 new or added lines in 3 files covered. (0.0%)

58 existing lines in 2 files now uncovered.

51 of 5490 relevant lines covered (0.93%)

0.01 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/src/pyVertexModel/analysis/analyse_simulation.py

import os
import pickle
import cv2

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit

from src.pyVertexModel.algorithm.vertexModel import VertexModel, logger
from src.pyVertexModel.util.utils import load_state, load_variables, save_variables


def analyse_simulation(folder):
    """
    Analyse the simulation results
    :param folder:
    :return:
    """

    # Check if the pkl file exists
    if not os.path.exists(os.path.join(folder, 'features_per_time.pkl')):
        # return None, None, None, None
        vModel = VertexModel(create_output_folder=False)

        features_per_time = []
        features_per_time_all_cells = []

        # Go through all the files in the folder
        all_files = os.listdir(folder)

        # Sort files by date
        all_files = sorted(all_files, key=lambda x: os.path.getmtime(os.path.join(folder, x)))

        # if all_files has less than 20 files, return None
        if len(all_files) < 20:
            return None, None, None, None

        for file_id, file in enumerate(all_files):
            if file.endswith('.pkl') and not file.__contains__('data_step_before_remodelling') and not file.__contains__('recoil'):
                # Load the state of the model
                load_state(vModel, os.path.join(folder, file))

                # Analyse the simulation
                all_cells, avg_cells = vModel.analyse_vertex_model()
                features_per_time_all_cells.append(all_cells)
                features_per_time.append(avg_cells)

                # # Create a temporary directory to store the images
                # temp_dir = os.path.join(folder, 'images')
                # if not os.path.exists(temp_dir):
                #     os.mkdir(temp_dir)
                # vModel.screenshot(temp_dir)

                # temp_dir = os.path.join(folder, 'images_wound_edge')
                # if not os.path.exists(temp_dir):
                #     os.mkdir(temp_dir)
                # _, debris_cells = vModel.geo.compute_wound_centre()
                # list_of_cell_distances_top = vModel.geo.compute_cell_distance_to_wound(debris_cells, location_filter=0)
                # alive_cells = [cell.ID for cell in vModel.geo.Cells if cell.AliveStatus == 1]
                # wound_edge_cells = []
                # for cell_num, cell_id in enumerate(alive_cells):
                #     if list_of_cell_distances_top[cell_num] == 1:
                #         wound_edge_cells.append(cell_id)
                # vModel.screenshot(temp_dir, wound_edge_cells)

        if not features_per_time:
            return None, None, None, None

        # Export to xlsx
        features_per_time_all_cells_df = pd.DataFrame(np.concatenate(features_per_time_all_cells),
                                                      columns=features_per_time_all_cells[0].columns)
        features_per_time_all_cells_df.sort_values(by='time', inplace=True)
        features_per_time_all_cells_df.to_excel(os.path.join(folder, 'features_per_time_all_cells.xlsx'))

        features_per_time_df = pd.DataFrame(features_per_time)
        features_per_time_df.sort_values(by='time', inplace=True)
        features_per_time_df.to_excel(os.path.join(folder, 'features_per_time.xlsx'))

        # Obtain pre-wound features
        try:
            pre_wound_features = features_per_time_df['time'][features_per_time_df['time'] < vModel.set.TInitAblation]
            pre_wound_features = features_per_time_df[features_per_time_df['time'] ==
                                                      pre_wound_features.iloc[-1]]
        except Exception as e:
            pre_wound_features = features_per_time_df.iloc[0]

        # Obtain post-wound features
        post_wound_features = features_per_time_df[features_per_time_df['time'] >= vModel.set.TInitAblation]

        if not post_wound_features.empty:
            # Reset time to ablation time.
            post_wound_features.loc[:, 'time'] = post_wound_features['time'] - vModel.set.TInitAblation

            # Compare post-wound features with pre-wound features in percentage
            for feature in post_wound_features.columns:
                if np.any(np.isnan(pre_wound_features[feature])) or np.any(np.isnan(post_wound_features[feature])):
                    continue

                if feature == 'time':
                    continue
                post_wound_features.loc[:, feature] = (post_wound_features[feature] / np.array(
                    pre_wound_features[feature])) * 100

            # Export to xlsx
            post_wound_features.to_excel(os.path.join(folder, 'post_wound_features.xlsx'))

            important_features = calculate_important_features(post_wound_features)
        else:
            important_features = {
                'max_recoiling_top': np.nan,
                'max_recoiling_time_top': np.nan,
                'min_height_change': np.nan,
                'min_height_change_time': np.nan,
            }

        # Export to xlsx
        df = pd.DataFrame([important_features])
        df.to_excel(os.path.join(folder, 'important_features.xlsx'))

        # Save dataframes to a single pkl
        with open(os.path.join(folder, 'features_per_time.pkl'), 'wb') as f:
            pickle.dump(features_per_time_df, f)
            pickle.dump(important_features, f)
            pickle.dump(post_wound_features, f)
            pickle.dump(features_per_time_all_cells_df, f)

    else:
        # Load dataframes from pkl
        with open(os.path.join(folder, 'features_per_time.pkl'), 'rb') as f:
            features_per_time_df = pickle.load(f)
            important_features = pickle.load(f)
            post_wound_features = pickle.load(f)
            features_per_time_all_cells_df = pickle.load(f)

        important_features = calculate_important_features(post_wound_features)

    # Plot wound area top evolution over time and save it to a file
    plot_feature(folder, post_wound_features, name='wound_area_top')
    plot_feature(folder, post_wound_features, name='num_cells_wound_edge_top')
    plot_feature(folder, post_wound_features, name='wound_indentation_top')
    plot_feature(folder, post_wound_features, name='wound_indentation_bottom')

    return features_per_time_df, post_wound_features, important_features, features_per_time_all_cells_df


def plot_feature(folder, post_wound_features, name='wound_area_top'):
    """
    Plot a feature and save it to a file
    :param folder:
    :param post_wound_features:
    :param name:
    :return:
    """
    plt.figure()
    plt.plot(post_wound_features['time'], post_wound_features[name])
    plt.xlabel('Time (h)')
    plt.ylabel(name)
    # Change axis limits
    if np.max(post_wound_features['time']) > 60:
        plt.xlim([0, np.max(post_wound_features['time'])])
    else:
        plt.xlim([0, 60])
    plt.ylim([0, 250])
    plt.savefig(os.path.join(folder, name + '.png'))
    plt.close()


def calculate_important_features(post_wound_features):
    """
    Calculate important features from the post-wound features
    :param post_wound_features:
    :return:
    """
    # Obtain important features for post-wound
    if not post_wound_features['wound_area_top'].empty and post_wound_features['time'].iloc[-1] > 4:
        important_features = {
            'max_recoiling_top': np.max(post_wound_features['wound_area_top']),
            'max_recoiling_time_top': np.array(post_wound_features['time'])[
                np.argmax(post_wound_features['wound_area_top'])],
            'min_recoiling_top': np.min(post_wound_features['wound_area_top']),
            'min_recoiling_time_top': np.array(post_wound_features['time'])[
                np.argmin(post_wound_features['wound_area_top'])],
            'min_height_change': np.min(post_wound_features['wound_height']),
            'min_height_change_time': np.array(post_wound_features['time'])[
                np.argmin(post_wound_features['wound_height'])],
            'last_area_top': post_wound_features['wound_area_top'].iloc[-1],
            'last_area_time_top': post_wound_features['time'].iloc[-1],
        }

        # Extrapolate features to a given time
        times_to_extrapolate = {3.0, 6.0, 9.0, 12.0, 15.0, 21.0, 30.0, 36.0, 45.0, 51.0, 60.0}
        columns_to_extrapolate = {'wound_area_top', 'wound_height'}  # post_wound_features.columns
        for feature in columns_to_extrapolate:
            for time in times_to_extrapolate:
                # Extrapolate results to a given time
                important_features[feature + '_extrapolated_' + str(time)] = np.interp(time,
                                                                                       post_wound_features['time'],
                                                                                       post_wound_features[feature])

        # # Get ratio from area the first time to the other times
        # for time in times_to_extrapolate:
        #     if time != 6.0:
        #         important_features['ratio_area_top_' + str(time)] = (
        #                 important_features['wound_area_top_extrapolated_' + str(time)] /
        #                 important_features['wound_area_top_extrapolated_6.0'])

    else:
        important_features = {
            'max_recoiling_top': np.nan,
            'max_recoiling_time_top': np.nan,
            'min_height_change': np.nan,
            'min_height_change_time': np.nan,
        }

    return important_features


def analyse_edge_recoil(file_name_v_model, type_of_ablation='recoil_edge_info_apical', n_ablations=2, location_filter=0, t_end=0.5):
    """
    Analyse how much an edge recoil if we ablate an edge of a cell
    :param type_of_ablation:
    :param t_end: Time to iterate after the ablation
    :param file_name_v_model: file nae of the Vertex model
    :param n_ablations: Number of ablations to perform
    :param location_filter: Location filter
    :return:
    """

    v_model = VertexModel(create_output_folder=False)
    load_state(v_model, file_name_v_model)

    # Cells to ablate
    # cell_to_ablate = np.random.choice(possible_cells_to_ablate, 1)
    cell_to_ablate = [v_model.geo.Cells[0]]

    #Pick the neighbouring cell to ablate
    neighbours = cell_to_ablate[0].compute_neighbours(location_filter)

    # Random order of neighbours
    np.random.seed(0)
    np.random.shuffle(neighbours)

    list_of_dicts_to_save = []
    for num_ablation in range(n_ablations):
        load_state(v_model, file_name_v_model)
        try:
            vars = load_variables(file_name_v_model.replace('before_ablation.pkl', type_of_ablation + '.pkl'))
            list_of_dicts_to_save_loaded = vars['recoiling_info_df_apical']

            cell_to_ablate_ID = list_of_dicts_to_save_loaded['cell_to_ablate'][num_ablation]
            neighbour_to_ablate_ID = list_of_dicts_to_save_loaded['neighbour_to_ablate'][num_ablation]
            edge_length_init = list_of_dicts_to_save_loaded['edge_length_init'][num_ablation]
            edge_length_final = list_of_dicts_to_save_loaded['edge_length_final'][num_ablation]
            if 'edge_length_final_normalized' in list_of_dicts_to_save_loaded:
                edge_length_final_normalized = list_of_dicts_to_save_loaded['edge_length_final_normalized'][
                    num_ablation]
            else:
                edge_length_final_normalized = (edge_length_final - edge_length_init) / edge_length_init

            initial_recoil = list_of_dicts_to_save_loaded['initial_recoil_in_s'][num_ablation]
            K = list_of_dicts_to_save_loaded['K'][num_ablation]
            scutoid_face = list_of_dicts_to_save_loaded['scutoid_face'][num_ablation]
            distance_to_centre = list_of_dicts_to_save_loaded['distance_to_centre'][num_ablation]
            if 'time_steps' in list_of_dicts_to_save_loaded:
                time_steps = list_of_dicts_to_save_loaded['time_steps'][num_ablation]
            else:
                time_steps = np.arange(0, len(edge_length_final)) * 6

            if edge_length_final[0] == 0:
                # Remove the first element
                edge_length_final = edge_length_final[1:]
                time_steps = time_steps[1:]
        except Exception as e:
            logger.info('Performing the analysis...' + str(e))
            # Change name of folder and create it
            if type_of_ablation == 'recoil_info_apical':
                v_model.set.OutputFolder = v_model.set.OutputFolder + '_ablation_' + str(num_ablation)
            else:
                v_model.set.OutputFolder = v_model.set.OutputFolder + '_ablation_edge_' + str(num_ablation)

            if not os.path.exists(v_model.set.OutputFolder):
                os.mkdir(v_model.set.OutputFolder)

            neighbour_to_ablate = [neighbours[num_ablation]]

            # Calculate if the cell is neighbour on both sides
            scutoid_face = None
            neighbours_other_side = []
            if location_filter == 0:
                neighbours_other_side = cell_to_ablate[0].compute_neighbours(location_filter=2)
                scutoid_face = np.nan
            elif location_filter == 2:
                neighbours_other_side = cell_to_ablate[0].compute_neighbours(location_filter=0)
                scutoid_face = np.nan

            if scutoid_face is not None:
                if neighbour_to_ablate[0] in neighbours_other_side:
                    scutoid_face = True
                else:
                    scutoid_face = False

            # Get the centre of the tissue
            centre_of_tissue = v_model.geo.compute_centre_of_tissue()
            neighbour_to_ablate_cell = [cell for cell in v_model.geo.Cells if cell.ID == neighbour_to_ablate[0]][0]
            distance_to_centre = np.mean([cell_to_ablate[0].compute_distance_to_centre(centre_of_tissue),
                                          neighbour_to_ablate_cell.compute_distance_to_centre(centre_of_tissue)])

            # Pick the neighbour and put it in the list
            cells_to_ablate = [cell_to_ablate[0].ID, neighbour_to_ablate[0]]

            # Get the edge that share both cells
            edge_length_init = v_model.geo.get_edge_length(cells_to_ablate, location_filter)

            # Ablate the edge
            v_model.set.ablation = True
            v_model.geo.cellsToAblate = cells_to_ablate
            v_model.set.TInitAblation = v_model.t
            if type_of_ablation == 'recoil_info_apical':
                v_model.geo.ablate_cells(v_model.set, v_model.t, combine_cells=False)
                v_model.geo.y_ablated = []
            elif type_of_ablation == 'recoil_edge_info_apical':
                v_model.geo.y_ablated = v_model.geo.ablate_edge(v_model.set, v_model.t, domain=location_filter,
                                                                adjacent_surface=False)

            # Relax the system
            initial_time = v_model.t
            v_model.set.tend = v_model.t + t_end
            if type_of_ablation == 'recoil_info_apical':
                v_model.set.dt = 0.005
            elif type_of_ablation == 'recoil_edge_info_apical':
                v_model.set.dt = 0.005

            v_model.set.Remodelling = False

            v_model.set.dt0 = v_model.set.dt
            if type_of_ablation == 'recoil_edge_info_apical':
                v_model.set.RemodelingFrequency = 0.05
            else:
                v_model.set.RemodelingFrequency = 100
            v_model.set.ablation = False
            v_model.set.export_images = True
            if v_model.set.export_images and not os.path.exists(v_model.set.OutputFolder + '/images'):
                os.mkdir(v_model.set.OutputFolder + '/images')
            edge_length_final_normalized = []
            edge_length_final = []
            recoil_speed = []
            time_steps = []

            # if os.path.exists(v_model.set.OutputFolder):
            #     list_of_files = os.listdir(v_model.set.OutputFolder)
            #     # Get file modification times and sort files by date
            #     files_with_dates = [(file, os.path.getmtime(os.path.join(v_model.set.OutputFolder, file))) for file in
            #                         list_of_files]
            #     files_with_dates.sort(key=lambda x: x[1])
            #     for file in files_with_dates:
            #         load_state(v_model, os.path.join(v_model.set.OutputFolder, file[0]))
            #         compute_edge_length_v_model(cells_to_ablate, edge_length_final, edge_length_final_normalized,
            #                                     edge_length_init, initial_time, location_filter, recoil_speed,
            #                                     time_steps,
            #                                     v_model)

            while v_model.t <= v_model.set.tend and not v_model.didNotConverge:
                gr = v_model.single_iteration()

                compute_edge_length_v_model(cells_to_ablate, edge_length_final, edge_length_final_normalized,
                                            edge_length_init, initial_time, location_filter, recoil_speed, time_steps,
                                            v_model)

                if np.isnan(gr):
                    break

            cell_to_ablate_ID = cell_to_ablate[0].ID
            neighbour_to_ablate_ID = neighbour_to_ablate[0]

        K, initial_recoil, error_bars = fit_ablation_equation(edge_length_final, time_steps)

        # Generate a plot with the edge length final and the fit for each ablation
        plt.figure()
        plt.plot(time_steps, edge_length_final_normalized, 'o')
        # Plot fit line of the Kelvin-Voigt model
        plt.plot(time_steps, recoil_model(np.array(time_steps), initial_recoil, K), 'r')
        plt.xlabel('Time (s)')
        plt.ylabel('Edge length final')
        plt.title('Ablation fit - ' + str(cell_to_ablate_ID) + ' ' + str(neighbour_to_ablate_ID))

        # Save plot
        if type_of_ablation == 'recoil_info_apical':
            plt.savefig(
                os.path.join(file_name_v_model.replace('before_ablation.pkl', 'ablation_fit_' + str(num_ablation) + '.png'))
            )
        elif type_of_ablation == 'recoil_edge_info_apical':
            plt.savefig(
                os.path.join(file_name_v_model.replace('before_ablation.pkl', 'ablation_edge_fit_' + str(num_ablation) + '.png'))
            )
        plt.close()

        # Save the results
        dict_to_save = {
            'cell_to_ablate': cell_to_ablate_ID,
            'neighbour_to_ablate': neighbour_to_ablate_ID,
            'edge_length_init': edge_length_init,
            'edge_length_final': edge_length_final,
            'edge_length_final_normalized': edge_length_final_normalized,
            'initial_recoil_in_s': initial_recoil,
            'K': K,
            'scutoid_face': scutoid_face,
            'location_filter': location_filter,
            'distance_to_centre': distance_to_centre,
            'time_steps': time_steps,
        }
        list_of_dicts_to_save.append(dict_to_save)

    recoiling_info_df_apical = pd.DataFrame(list_of_dicts_to_save)
    recoiling_info_df_apical.to_excel(file_name_v_model.replace('before_ablation.pkl', type_of_ablation+'.xlsx'))
    save_variables({'recoiling_info_df_apical': recoiling_info_df_apical},
                   file_name_v_model.replace('before_ablation.pkl', type_of_ablation+'.pkl'))

    return list_of_dicts_to_save


def recoil_model(x, initial_recoil, K):
    """
    Model of the recoil based on a Kelvin-Voigt model
    :param x:
    :param initial_recoil:
    :param K:
    :return:   Recoil
    """
    return (initial_recoil / K) * (1 - np.exp(-K * x))


def fit_ablation_equation(edge_length_final_normalized, time_steps):
    """
    Fit the ablation equation. Thanks to Veronika Lachina.
    :param edge_length_final_normalized:
    :param time_steps:
    :return:    K, initial_recoil
    """

    # Normalize the edge length
    edge_length_init = edge_length_final_normalized[0]
    edge_length_final_normalized = (edge_length_final_normalized - edge_length_init) / edge_length_init

    # Fit the model to the data
    [params, covariance] = curve_fit(recoil_model, time_steps, edge_length_final_normalized,
                                     p0=[0.00001, 3], bounds=(0, np.inf))

    # Get the error
    error_bars = np.sqrt(np.diag(covariance))

    initial_recoil, K = params
    return K, initial_recoil, error_bars


def compute_edge_length_v_model(cells_to_ablate, edge_length_final, edge_length_final_normalized, edge_length_init,
                                initial_time, location_filter, recoil_speed, time_steps, v_model):
    """
    Compute the edge length of the edge that share the cells_to_ablate
    :param cells_to_ablate:
    :param edge_length_final:
    :param edge_length_final_normalized:
    :param edge_length_init:
    :param initial_time:
    :param location_filter:
    :param recoil_speed:
    :param time_steps:
    :param v_model:
    :return:
    """
    if v_model.t == initial_time:
        return
    # Get the edge length
    edge_length_final.append(v_model.geo.get_edge_length(cells_to_ablate, location_filter))
    edge_length_final_normalized.append((edge_length_final[-1] - edge_length_init) / edge_length_init)
    print('Edge length final: ', edge_length_final[-1])
    # In seconds. 1 t = 1 minute = 60 seconds
    time_steps.append((v_model.t - initial_time) * 60)
    # Calculate the recoil
    recoil_speed.append(edge_length_final_normalized[-1] / time_steps[-1])

def create_video(folder, name_containing_images='top_'):
    """
    Create a video of the images in the folder
    :param folder:
    :return:
    """

    # Get the images
    images = [img for img in os.listdir(folder) if img.endswith(".png") and name_containing_images in img]
    images.sort(key=lambda x: int(x.split('_')[-1].split('.')[0]))

    # Determine the width and height from the first image
    image_path = os.path.join(folder, images[0])
    frame = cv2.imread(image_path)
    height, width, layers = frame.shape

    # Define the codec and create VideoWriter object
    fourcc = cv2.VideoWriter_fourcc(*'mp4v')  # Use 'mp4v' for MP4
    video = cv2.VideoWriter(os.path.join(folder, name_containing_images + 'video.mp4'), fourcc, 7, (width, height))

    for image in images:
        video.write(cv2.imread(os.path.join(folder, image)))

    cv2.destroyAllWindows()
    video.release()

1	import os	×
2	import pickle	×
3	import cv2	×
4
5	import numpy as np	×
6	import pandas as pd	×
7	from matplotlib import pyplot as plt	×
8	from scipy.optimize import curve_fit	×
9
10	from src.pyVertexModel.algorithm.vertexModel import VertexModel, logger	×
11	from src.pyVertexModel.util.utils import load_state, load_variables, save_variables	×
12
13
14	def analyse_simulation(folder):	×
15	"""
16	Analyse the simulation results
17	:param folder:
18	:return:
19	"""
20
21	# Check if the pkl file exists
22	if not os.path.exists(os.path.join(folder, 'features_per_time.pkl')):	×
23	# return None, None, None, None
24	vModel = VertexModel(create_output_folder=False)	×
25
26	features_per_time = []	×
27	features_per_time_all_cells = []	×
28
29	# Go through all the files in the folder
30	all_files = os.listdir(folder)	×
31
32	# Sort files by date
33	all_files = sorted(all_files, key=lambda x: os.path.getmtime(os.path.join(folder, x)))	×
34
35	# if all_files has less than 20 files, return None
36	if len(all_files) < 20:	×
37	return None, None, None, None	×
38
39	for file_id, file in enumerate(all_files):	×
40	if file.endswith('.pkl') and not file.__contains__('data_step_before_remodelling') and not file.__contains__('recoil'):	×
41	# Load the state of the model
42	load_state(vModel, os.path.join(folder, file))	×
43
44	# Analyse the simulation
45	all_cells, avg_cells = vModel.analyse_vertex_model()	×
46	features_per_time_all_cells.append(all_cells)	×
47	features_per_time.append(avg_cells)	×
48
49	# # Create a temporary directory to store the images
50	# temp_dir = os.path.join(folder, 'images')
51	# if not os.path.exists(temp_dir):
52	# os.mkdir(temp_dir)
53	# vModel.screenshot(temp_dir)
54
55	# temp_dir = os.path.join(folder, 'images_wound_edge')
56	# if not os.path.exists(temp_dir):
57	# os.mkdir(temp_dir)
58	# _, debris_cells = vModel.geo.compute_wound_centre()
59	# list_of_cell_distances_top = vModel.geo.compute_cell_distance_to_wound(debris_cells, location_filter=0)
60	# alive_cells = [cell.ID for cell in vModel.geo.Cells if cell.AliveStatus == 1]
61	# wound_edge_cells = []
62	# for cell_num, cell_id in enumerate(alive_cells):
63	# if list_of_cell_distances_top[cell_num] == 1:
64	# wound_edge_cells.append(cell_id)
65	# vModel.screenshot(temp_dir, wound_edge_cells)
66
67	if not features_per_time:	×
68	return None, None, None, None	×
69
70	# Export to xlsx
71	features_per_time_all_cells_df = pd.DataFrame(np.concatenate(features_per_time_all_cells),	×
72	columns=features_per_time_all_cells[0].columns)
73	features_per_time_all_cells_df.sort_values(by='time', inplace=True)	×
74	features_per_time_all_cells_df.to_excel(os.path.join(folder, 'features_per_time_all_cells.xlsx'))	×
75
76	features_per_time_df = pd.DataFrame(features_per_time)	×
77	features_per_time_df.sort_values(by='time', inplace=True)	×
78	features_per_time_df.to_excel(os.path.join(folder, 'features_per_time.xlsx'))	×
79
80	# Obtain pre-wound features
81	try:	×
82	pre_wound_features = features_per_time_df['time'][features_per_time_df['time'] < vModel.set.TInitAblation]	×
83	pre_wound_features = features_per_time_df[features_per_time_df['time'] ==	×
84	pre_wound_features.iloc[-1]]
85	except Exception as e:	×
86	pre_wound_features = features_per_time_df.iloc[0]	×
87
88	# Obtain post-wound features
89	post_wound_features = features_per_time_df[features_per_time_df['time'] >= vModel.set.TInitAblation]	×
90
91	if not post_wound_features.empty:	×
92	# Reset time to ablation time.
93	post_wound_features.loc[:, 'time'] = post_wound_features['time'] - vModel.set.TInitAblation	×
94
95	# Compare post-wound features with pre-wound features in percentage
96	for feature in post_wound_features.columns:	×
97	if np.any(np.isnan(pre_wound_features[feature])) or np.any(np.isnan(post_wound_features[feature])):	×
98	continue	×
99
100	if feature == 'time':	×
101	continue	×
102	post_wound_features.loc[:, feature] = (post_wound_features[feature] / np.array(	×
103	pre_wound_features[feature])) * 100
104
105	# Export to xlsx
106	post_wound_features.to_excel(os.path.join(folder, 'post_wound_features.xlsx'))	×
107
108	important_features = calculate_important_features(post_wound_features)	×
109	else:
110	important_features = {	×
111	'max_recoiling_top': np.nan,
112	'max_recoiling_time_top': np.nan,
113	'min_height_change': np.nan,
114	'min_height_change_time': np.nan,
115	}
116
117	# Export to xlsx
118	df = pd.DataFrame([important_features])	×
119	df.to_excel(os.path.join(folder, 'important_features.xlsx'))	×
120
121	# Save dataframes to a single pkl
122	with open(os.path.join(folder, 'features_per_time.pkl'), 'wb') as f:	×
123	pickle.dump(features_per_time_df, f)	×
124	pickle.dump(important_features, f)	×
125	pickle.dump(post_wound_features, f)	×
126	pickle.dump(features_per_time_all_cells_df, f)	×
127
128	else:
129	# Load dataframes from pkl
130	with open(os.path.join(folder, 'features_per_time.pkl'), 'rb') as f:	×
131	features_per_time_df = pickle.load(f)	×
132	important_features = pickle.load(f)	×
133	post_wound_features = pickle.load(f)	×
134	features_per_time_all_cells_df = pickle.load(f)	×
135
136	important_features = calculate_important_features(post_wound_features)	×
137
138	# Plot wound area top evolution over time and save it to a file
139	plot_feature(folder, post_wound_features, name='wound_area_top')	×
140	plot_feature(folder, post_wound_features, name='num_cells_wound_edge_top')	×
NEW 141	plot_feature(folder, post_wound_features, name='wound_indentation_top')	×
NEW 142	plot_feature(folder, post_wound_features, name='wound_indentation_bottom')	×
143
144	return features_per_time_df, post_wound_features, important_features, features_per_time_all_cells_df	×
145
146
147	def plot_feature(folder, post_wound_features, name='wound_area_top'):	×
148	"""
149	Plot a feature and save it to a file
150	:param folder:
151	:param post_wound_features:
152	:param name:
153	:return:
154	"""
155	plt.figure()	×
156	plt.plot(post_wound_features['time'], post_wound_features[name])	×
157	plt.xlabel('Time (h)')	×
158	plt.ylabel(name)	×
159	# Change axis limits
160	if np.max(post_wound_features['time']) > 60:	×
161	plt.xlim([0, np.max(post_wound_features['time'])])	×
162	else:
163	plt.xlim([0, 60])	×
164	plt.ylim([0, 250])	×
165	plt.savefig(os.path.join(folder, name + '.png'))	×
166	plt.close()	×
167
168
169	def calculate_important_features(post_wound_features):	×
170	"""
171	Calculate important features from the post-wound features
172	:param post_wound_features:
173	:return:
174	"""
175	# Obtain important features for post-wound
176	if not post_wound_features['wound_area_top'].empty and post_wound_features['time'].iloc[-1] > 4:	×
177	important_features = {	×
178	'max_recoiling_top': np.max(post_wound_features['wound_area_top']),
179	'max_recoiling_time_top': np.array(post_wound_features['time'])[
180	np.argmax(post_wound_features['wound_area_top'])],
181	'min_recoiling_top': np.min(post_wound_features['wound_area_top']),
182	'min_recoiling_time_top': np.array(post_wound_features['time'])[
183	np.argmin(post_wound_features['wound_area_top'])],
184	'min_height_change': np.min(post_wound_features['wound_height']),
185	'min_height_change_time': np.array(post_wound_features['time'])[
186	np.argmin(post_wound_features['wound_height'])],
187	'last_area_top': post_wound_features['wound_area_top'].iloc[-1],
188	'last_area_time_top': post_wound_features['time'].iloc[-1],
189	}
190
191	# Extrapolate features to a given time
192	times_to_extrapolate = {3.0, 6.0, 9.0, 12.0, 15.0, 21.0, 30.0, 36.0, 45.0, 51.0, 60.0}	×
193	columns_to_extrapolate = {'wound_area_top', 'wound_height'} # post_wound_features.columns	×
194	for feature in columns_to_extrapolate:	×
195	for time in times_to_extrapolate:	×
196	# Extrapolate results to a given time
197	important_features[feature + '_extrapolated_' + str(time)] = np.interp(time,	×
198	post_wound_features['time'],
199	post_wound_features[feature])
200
201	# # Get ratio from area the first time to the other times
202	# for time in times_to_extrapolate:
203	# if time != 6.0:
204	# important_features['ratio_area_top_' + str(time)] = (
205	# important_features['wound_area_top_extrapolated_' + str(time)] /
206	# important_features['wound_area_top_extrapolated_6.0'])
207
208	else:
209	important_features = {	×
210	'max_recoiling_top': np.nan,
211	'max_recoiling_time_top': np.nan,
212	'min_height_change': np.nan,
213	'min_height_change_time': np.nan,
214	}
215
216	return important_features	×
217
218
219	def analyse_edge_recoil(file_name_v_model, type_of_ablation='recoil_edge_info_apical', n_ablations=2, location_filter=0, t_end=0.5):	×
220	"""
221	Analyse how much an edge recoil if we ablate an edge of a cell
222	:param type_of_ablation:
223	:param t_end: Time to iterate after the ablation
224	:param file_name_v_model: file nae of the Vertex model
225	:param n_ablations: Number of ablations to perform
226	:param location_filter: Location filter
227	:return:
228	"""
229
230	v_model = VertexModel(create_output_folder=False)	×
231	load_state(v_model, file_name_v_model)	×
232
233	# Cells to ablate
234	# cell_to_ablate = np.random.choice(possible_cells_to_ablate, 1)
235	cell_to_ablate = [v_model.geo.Cells[0]]	×
236
237	#Pick the neighbouring cell to ablate
238	neighbours = cell_to_ablate[0].compute_neighbours(location_filter)	×
239
240	# Random order of neighbours
241	np.random.seed(0)	×
242	np.random.shuffle(neighbours)	×
243
244	list_of_dicts_to_save = []	×
245	for num_ablation in range(n_ablations):	×
246	load_state(v_model, file_name_v_model)	×
247	try:	×
248	vars = load_variables(file_name_v_model.replace('before_ablation.pkl', type_of_ablation + '.pkl'))	×
249	list_of_dicts_to_save_loaded = vars['recoiling_info_df_apical']	×
250
251	cell_to_ablate_ID = list_of_dicts_to_save_loaded['cell_to_ablate'][num_ablation]	×
252	neighbour_to_ablate_ID = list_of_dicts_to_save_loaded['neighbour_to_ablate'][num_ablation]	×
253	edge_length_init = list_of_dicts_to_save_loaded['edge_length_init'][num_ablation]	×
254	edge_length_final = list_of_dicts_to_save_loaded['edge_length_final'][num_ablation]	×
255	if 'edge_length_final_normalized' in list_of_dicts_to_save_loaded:	×
256	edge_length_final_normalized = list_of_dicts_to_save_loaded['edge_length_final_normalized'][	×
257	num_ablation]
258	else:
259	edge_length_final_normalized = (edge_length_final - edge_length_init) / edge_length_init	×
260
261	initial_recoil = list_of_dicts_to_save_loaded['initial_recoil_in_s'][num_ablation]	×
262	K = list_of_dicts_to_save_loaded['K'][num_ablation]	×
263	scutoid_face = list_of_dicts_to_save_loaded['scutoid_face'][num_ablation]	×
264	distance_to_centre = list_of_dicts_to_save_loaded['distance_to_centre'][num_ablation]	×
265	if 'time_steps' in list_of_dicts_to_save_loaded:	×
266	time_steps = list_of_dicts_to_save_loaded['time_steps'][num_ablation]	×
267	else:
268	time_steps = np.arange(0, len(edge_length_final)) * 6	×
269
270	if edge_length_final[0] == 0:	×
271	# Remove the first element
272	edge_length_final = edge_length_final[1:]	×
273	time_steps = time_steps[1:]	×
274	except Exception as e:	×
275	logger.info('Performing the analysis...' + str(e))	×
276	# Change name of folder and create it
277	if type_of_ablation == 'recoil_info_apical':	×
278	v_model.set.OutputFolder = v_model.set.OutputFolder + '_ablation_' + str(num_ablation)	×
279	else:
280	v_model.set.OutputFolder = v_model.set.OutputFolder + '_ablation_edge_' + str(num_ablation)	×
281
282	if not os.path.exists(v_model.set.OutputFolder):	×
283	os.mkdir(v_model.set.OutputFolder)	×
284
285	neighbour_to_ablate = [neighbours[num_ablation]]	×
286
287	# Calculate if the cell is neighbour on both sides
288	scutoid_face = None	×
289	neighbours_other_side = []	×
290	if location_filter == 0:	×
291	neighbours_other_side = cell_to_ablate[0].compute_neighbours(location_filter=2)	×
292	scutoid_face = np.nan	×
293	elif location_filter == 2:	×
294	neighbours_other_side = cell_to_ablate[0].compute_neighbours(location_filter=0)	×
295	scutoid_face = np.nan	×
296
297	if scutoid_face is not None:	×
298	if neighbour_to_ablate[0] in neighbours_other_side:	×
299	scutoid_face = True	×
300	else:
301	scutoid_face = False	×
302
303	# Get the centre of the tissue
304	centre_of_tissue = v_model.geo.compute_centre_of_tissue()	×
305	neighbour_to_ablate_cell = [cell for cell in v_model.geo.Cells if cell.ID == neighbour_to_ablate[0]][0]	×
306	distance_to_centre = np.mean([cell_to_ablate[0].compute_distance_to_centre(centre_of_tissue),	×
307	neighbour_to_ablate_cell.compute_distance_to_centre(centre_of_tissue)])
308
309	# Pick the neighbour and put it in the list
310	cells_to_ablate = [cell_to_ablate[0].ID, neighbour_to_ablate[0]]	×
311
312	# Get the edge that share both cells
313	edge_length_init = v_model.geo.get_edge_length(cells_to_ablate, location_filter)	×
314
315	# Ablate the edge
316	v_model.set.ablation = True	×
317	v_model.geo.cellsToAblate = cells_to_ablate	×
318	v_model.set.TInitAblation = v_model.t	×
319	if type_of_ablation == 'recoil_info_apical':	×
320	v_model.geo.ablate_cells(v_model.set, v_model.t, combine_cells=False)	×
321	v_model.geo.y_ablated = []	×
322	elif type_of_ablation == 'recoil_edge_info_apical':	×
323	v_model.geo.y_ablated = v_model.geo.ablate_edge(v_model.set, v_model.t, domain=location_filter,	×
324	adjacent_surface=False)
325
326	# Relax the system
327	initial_time = v_model.t	×
328	v_model.set.tend = v_model.t + t_end	×
329	if type_of_ablation == 'recoil_info_apical':	×
330	v_model.set.dt = 0.005	×
331	elif type_of_ablation == 'recoil_edge_info_apical':	×
332	v_model.set.dt = 0.005	×
333
334	v_model.set.Remodelling = False	×
335
336	v_model.set.dt0 = v_model.set.dt	×
337	if type_of_ablation == 'recoil_edge_info_apical':	×
338	v_model.set.RemodelingFrequency = 0.05	×
339	else:
340	v_model.set.RemodelingFrequency = 100	×
341	v_model.set.ablation = False	×
342	v_model.set.export_images = True	×
343	if v_model.set.export_images and not os.path.exists(v_model.set.OutputFolder + '/images'):	×
344	os.mkdir(v_model.set.OutputFolder + '/images')	×
345	edge_length_final_normalized = []	×
346	edge_length_final = []	×
347	recoil_speed = []	×
348	time_steps = []	×
349
350	# if os.path.exists(v_model.set.OutputFolder):
351	# list_of_files = os.listdir(v_model.set.OutputFolder)
352	# # Get file modification times and sort files by date
353	# files_with_dates = [(file, os.path.getmtime(os.path.join(v_model.set.OutputFolder, file))) for file in
354	# list_of_files]
355	# files_with_dates.sort(key=lambda x: x[1])
356	# for file in files_with_dates:
357	# load_state(v_model, os.path.join(v_model.set.OutputFolder, file[0]))
358	# compute_edge_length_v_model(cells_to_ablate, edge_length_final, edge_length_final_normalized,
359	# edge_length_init, initial_time, location_filter, recoil_speed,
360	# time_steps,
361	# v_model)
362
363	while v_model.t <= v_model.set.tend and not v_model.didNotConverge:	×
364	gr = v_model.single_iteration()	×
365
366	compute_edge_length_v_model(cells_to_ablate, edge_length_final, edge_length_final_normalized,	×
367	edge_length_init, initial_time, location_filter, recoil_speed, time_steps,
368	v_model)
369
370	if np.isnan(gr):	×
371	break	×
372
373	cell_to_ablate_ID = cell_to_ablate[0].ID	×
374	neighbour_to_ablate_ID = neighbour_to_ablate[0]	×
375
376	K, initial_recoil, error_bars = fit_ablation_equation(edge_length_final, time_steps)	×
377
378	# Generate a plot with the edge length final and the fit for each ablation
379	plt.figure()	×
380	plt.plot(time_steps, edge_length_final_normalized, 'o')	×
381	# Plot fit line of the Kelvin-Voigt model
382	plt.plot(time_steps, recoil_model(np.array(time_steps), initial_recoil, K), 'r')	×
383	plt.xlabel('Time (s)')	×
384	plt.ylabel('Edge length final')	×
385	plt.title('Ablation fit - ' + str(cell_to_ablate_ID) + ' ' + str(neighbour_to_ablate_ID))	×
386
387	# Save plot
388	if type_of_ablation == 'recoil_info_apical':	×
389	plt.savefig(	×
390	os.path.join(file_name_v_model.replace('before_ablation.pkl', 'ablation_fit_' + str(num_ablation) + '.png'))
391	)
392	elif type_of_ablation == 'recoil_edge_info_apical':	×
393	plt.savefig(	×
394	os.path.join(file_name_v_model.replace('before_ablation.pkl', 'ablation_edge_fit_' + str(num_ablation) + '.png'))
395	)
396	plt.close()	×
397
398	# Save the results
399	dict_to_save = {	×
400	'cell_to_ablate': cell_to_ablate_ID,
401	'neighbour_to_ablate': neighbour_to_ablate_ID,
402	'edge_length_init': edge_length_init,
403	'edge_length_final': edge_length_final,
404	'edge_length_final_normalized': edge_length_final_normalized,
405	'initial_recoil_in_s': initial_recoil,
406	'K': K,
407	'scutoid_face': scutoid_face,
408	'location_filter': location_filter,
409	'distance_to_centre': distance_to_centre,
410	'time_steps': time_steps,
411	}
412	list_of_dicts_to_save.append(dict_to_save)	×
413
414	recoiling_info_df_apical = pd.DataFrame(list_of_dicts_to_save)	×
415	recoiling_info_df_apical.to_excel(file_name_v_model.replace('before_ablation.pkl', type_of_ablation+'.xlsx'))	×
416	save_variables({'recoiling_info_df_apical': recoiling_info_df_apical},	×
417	file_name_v_model.replace('before_ablation.pkl', type_of_ablation+'.pkl'))
418
419	return list_of_dicts_to_save	×
420
421
422	def recoil_model(x, initial_recoil, K):	×
423	"""
424	Model of the recoil based on a Kelvin-Voigt model
425	:param x:
426	:param initial_recoil:
427	:param K:
428	:return: Recoil
429	"""
430	return (initial_recoil / K) * (1 - np.exp(-K * x))	×
431
432
433	def fit_ablation_equation(edge_length_final_normalized, time_steps):	×
434	"""
435	Fit the ablation equation. Thanks to Veronika Lachina.
436	:param edge_length_final_normalized:
437	:param time_steps:
438	:return: K, initial_recoil
439	"""
440
441	# Normalize the edge length
442	edge_length_init = edge_length_final_normalized[0]	×
443	edge_length_final_normalized = (edge_length_final_normalized - edge_length_init) / edge_length_init	×
444
445	# Fit the model to the data
446	[params, covariance] = curve_fit(recoil_model, time_steps, edge_length_final_normalized,	×
447	p0=[0.00001, 3], bounds=(0, np.inf))
448
449	# Get the error
450	error_bars = np.sqrt(np.diag(covariance))	×
451
452	initial_recoil, K = params	×
453	return K, initial_recoil, error_bars	×
454
455
456	def compute_edge_length_v_model(cells_to_ablate, edge_length_final, edge_length_final_normalized, edge_length_init,	×
457	initial_time, location_filter, recoil_speed, time_steps, v_model):
458	"""
459	Compute the edge length of the edge that share the cells_to_ablate
460	:param cells_to_ablate:
461	:param edge_length_final:
462	:param edge_length_final_normalized:
463	:param edge_length_init:
464	:param initial_time:
465	:param location_filter:
466	:param recoil_speed:
467	:param time_steps:
468	:param v_model:
469	:return:
470	"""
471	if v_model.t == initial_time:	×
472	return	×
473	# Get the edge length
474	edge_length_final.append(v_model.geo.get_edge_length(cells_to_ablate, location_filter))	×
475	edge_length_final_normalized.append((edge_length_final[-1] - edge_length_init) / edge_length_init)	×
476	print('Edge length final: ', edge_length_final[-1])	×
477	# In seconds. 1 t = 1 minute = 60 seconds
478	time_steps.append((v_model.t - initial_time) * 60)	×
479	# Calculate the recoil
480	recoil_speed.append(edge_length_final_normalized[-1] / time_steps[-1])	×
481
482	def create_video(folder, name_containing_images='top_'):	×
483	"""
484	Create a video of the images in the folder
485	:param folder:
486	:return:
487	"""
488
489	# Get the images
490	images = [img for img in os.listdir(folder) if img.endswith(".png") and name_containing_images in img]	×
491	images.sort(key=lambda x: int(x.split('_')[-1].split('.')[0]))	×
492
493	# Determine the width and height from the first image
494	image_path = os.path.join(folder, images[0])	×
495	frame = cv2.imread(image_path)	×
496	height, width, layers = frame.shape	×
497
498	# Define the codec and create VideoWriter object
499	fourcc = cv2.VideoWriter_fourcc(*'mp4v') # Use 'mp4v' for MP4	×
500	video = cv2.VideoWriter(os.path.join(folder, name_containing_images + 'video.mp4'), fourcc, 7, (width, height))	×
501
502	for image in images:	×
503	video.write(cv2.imread(os.path.join(folder, image)))	×
504
505	cv2.destroyAllWindows()	×
506	video.release()	×

Pablo1990 / pyVertexModel / 11858895593

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous