Add files via upload
This commit is contained in:
GitHub
parent
2e5879cd62
commit
6599c9c81d
9 changed files with 772 additions and 188 deletions
|
|
@ -5,6 +5,8 @@
|
||||||
"raw_path": "raw",
|
"raw_path": "raw",
|
||||||
"export_path": "output_M3879M_2021-10-05",
|
"export_path": "output_M3879M_2021-10-05",
|
||||||
"ref_image_path": "ref_images_M3879M_2021-10-05",
|
"ref_image_path": "ref_images_M3879M_2021-10-05",
|
||||||
|
"heartbeat_remove": true, // if gevi must be true; geci: who knows...
|
||||||
|
"gevi": true, // true => gevi, false => geci
|
||||||
// Ratio Sequence
|
// Ratio Sequence
|
||||||
"classical_ratio_mode": true, // true: a/d false: 1+a-d
|
"classical_ratio_mode": true, // true: a/d false: 1+a-d
|
||||||
// Regression
|
// Regression
|
||||||
|
|
|
||||||
|
|
@ -5,6 +5,8 @@
|
||||||
"raw_path": "raw",
|
"raw_path": "raw",
|
||||||
"export_path": "output_M_Sert_Cre_41",
|
"export_path": "output_M_Sert_Cre_41",
|
||||||
"ref_image_path": "ref_images_M_Sert_Cre_41",
|
"ref_image_path": "ref_images_M_Sert_Cre_41",
|
||||||
|
"heartbeat_remove": false,
|
||||||
|
"gevi": false, // true => gevi, false => geci
|
||||||
// Ratio Sequence
|
// Ratio Sequence
|
||||||
"classical_ratio_mode": true, // true: a/d false: 1+a-d
|
"classical_ratio_mode": true, // true: a/d false: 1+a-d
|
||||||
// Regression
|
// Regression
|
||||||
|
|
|
||||||
62
config_M_Sert_Cre_42.json
Normal file
62
config_M_Sert_Cre_42.json
Normal file
|
|
@ -0,0 +1,62 @@
|
||||||
|
{
|
||||||
|
"basic_path": "/data_1/hendrik",
|
||||||
|
"recoding_data": "2023-07-18",
|
||||||
|
"mouse_identifier": "M_Sert_Cre_42",
|
||||||
|
"raw_path": "raw",
|
||||||
|
"export_path": "output_M_Sert_Cre_42",
|
||||||
|
"ref_image_path": "ref_images_M_Sert_Cre_42",
|
||||||
|
"heartbeat_remove": false,
|
||||||
|
"gevi": false, // true => gevi, false => geci
|
||||||
|
// Ratio Sequence
|
||||||
|
"classical_ratio_mode": true, // true: a/d false: 1+a-d
|
||||||
|
// Regression
|
||||||
|
//"target_camera_acceptor": "acceptor",
|
||||||
|
"target_camera_acceptor": "",
|
||||||
|
"regressor_cameras_acceptor": [
|
||||||
|
"oxygenation",
|
||||||
|
"volume"
|
||||||
|
],
|
||||||
|
"target_camera_donor": "donor",
|
||||||
|
"regressor_cameras_donor": [
|
||||||
|
// "oxygenation",
|
||||||
|
"volume"
|
||||||
|
],
|
||||||
|
// binning
|
||||||
|
"binning_enable": true,
|
||||||
|
"binning_at_the_end": false,
|
||||||
|
"binning_kernel_size": 4,
|
||||||
|
"binning_stride": 4,
|
||||||
|
"binning_divisor_override": 1,
|
||||||
|
// alignment
|
||||||
|
"alignment_batch_size": 200,
|
||||||
|
"rotation_stabilization_threshold_factor": 3.0, // >= 1.0
|
||||||
|
"rotation_stabilization_threshold_border": 0.9, // <= 1.0
|
||||||
|
// Heart beat detection
|
||||||
|
"lower_freqency_bandpass": 5.0, // Hz
|
||||||
|
"upper_freqency_bandpass": 14.0, // Hz
|
||||||
|
"heartbeat_filtfilt_chuck_size": 10,
|
||||||
|
// Gauss smear
|
||||||
|
"gauss_smear_spatial_width": 8,
|
||||||
|
"gauss_smear_temporal_width": 0.1,
|
||||||
|
"gauss_smear_use_matlab_mask": false,
|
||||||
|
// LED Ramp on
|
||||||
|
"skip_frames_in_the_beginning": 100, // Frames
|
||||||
|
// PyTorch
|
||||||
|
"dtype": "float32",
|
||||||
|
"force_to_cpu": false,
|
||||||
|
// Save
|
||||||
|
"save_as_python": true, // produces .npz files (compressed)
|
||||||
|
"save_as_matlab": false, // produces .hd5 file (compressed)
|
||||||
|
// Save extra information
|
||||||
|
"save_alignment": false,
|
||||||
|
"save_heartbeat": false,
|
||||||
|
"save_factors": false,
|
||||||
|
"save_regression_coefficients": false,
|
||||||
|
// Not important parameter
|
||||||
|
"required_order": [
|
||||||
|
"acceptor",
|
||||||
|
"donor",
|
||||||
|
"oxygenation",
|
||||||
|
"volume"
|
||||||
|
]
|
||||||
|
}
|
||||||
62
config_M_Sert_Cre_45.json
Normal file
62
config_M_Sert_Cre_45.json
Normal file
|
|
@ -0,0 +1,62 @@
|
||||||
|
{
|
||||||
|
"basic_path": "/data_1/hendrik",
|
||||||
|
"recoding_data": "2023-07-18",
|
||||||
|
"mouse_identifier": "M_Sert_Cre_45",
|
||||||
|
"raw_path": "raw",
|
||||||
|
"export_path": "output_M_Sert_Cre_45",
|
||||||
|
"ref_image_path": "ref_images_M_Sert_Cre_45",
|
||||||
|
"heartbeat_remove": false,
|
||||||
|
"gevi": false, // true => gevi, false => geci
|
||||||
|
// Ratio Sequence
|
||||||
|
"classical_ratio_mode": true, // true: a/d false: 1+a-d
|
||||||
|
// Regression
|
||||||
|
//"target_camera_acceptor": "acceptor",
|
||||||
|
"target_camera_acceptor": "",
|
||||||
|
"regressor_cameras_acceptor": [
|
||||||
|
"oxygenation",
|
||||||
|
"volume"
|
||||||
|
],
|
||||||
|
"target_camera_donor": "donor",
|
||||||
|
"regressor_cameras_donor": [
|
||||||
|
// "oxygenation",
|
||||||
|
"volume"
|
||||||
|
],
|
||||||
|
// binning
|
||||||
|
"binning_enable": true,
|
||||||
|
"binning_at_the_end": false,
|
||||||
|
"binning_kernel_size": 4,
|
||||||
|
"binning_stride": 4,
|
||||||
|
"binning_divisor_override": 1,
|
||||||
|
// alignment
|
||||||
|
"alignment_batch_size": 200,
|
||||||
|
"rotation_stabilization_threshold_factor": 3.0, // >= 1.0
|
||||||
|
"rotation_stabilization_threshold_border": 0.9, // <= 1.0
|
||||||
|
// Heart beat detection
|
||||||
|
"lower_freqency_bandpass": 5.0, // Hz
|
||||||
|
"upper_freqency_bandpass": 14.0, // Hz
|
||||||
|
"heartbeat_filtfilt_chuck_size": 10,
|
||||||
|
// Gauss smear
|
||||||
|
"gauss_smear_spatial_width": 8,
|
||||||
|
"gauss_smear_temporal_width": 0.1,
|
||||||
|
"gauss_smear_use_matlab_mask": false,
|
||||||
|
// LED Ramp on
|
||||||
|
"skip_frames_in_the_beginning": 100, // Frames
|
||||||
|
// PyTorch
|
||||||
|
"dtype": "float32",
|
||||||
|
"force_to_cpu": false,
|
||||||
|
// Save
|
||||||
|
"save_as_python": true, // produces .npz files (compressed)
|
||||||
|
"save_as_matlab": false, // produces .hd5 file (compressed)
|
||||||
|
// Save extra information
|
||||||
|
"save_alignment": false,
|
||||||
|
"save_heartbeat": false,
|
||||||
|
"save_factors": false,
|
||||||
|
"save_regression_coefficients": false,
|
||||||
|
// Not important parameter
|
||||||
|
"required_order": [
|
||||||
|
"acceptor",
|
||||||
|
"donor",
|
||||||
|
"oxygenation",
|
||||||
|
"volume"
|
||||||
|
]
|
||||||
|
}
|
||||||
62
config_M_Sert_Cre_46.json
Normal file
62
config_M_Sert_Cre_46.json
Normal file
|
|
@ -0,0 +1,62 @@
|
||||||
|
{
|
||||||
|
"basic_path": "/data_1/hendrik",
|
||||||
|
"recoding_data": "2023-03-16",
|
||||||
|
"mouse_identifier": "M_Sert_Cre_46",
|
||||||
|
"raw_path": "raw",
|
||||||
|
"export_path": "output_M_Sert_Cre_46",
|
||||||
|
"ref_image_path": "ref_images_M_Sert_Cre_46",
|
||||||
|
"heartbeat_remove": false,
|
||||||
|
"gevi": false, // true => gevi, false => geci
|
||||||
|
// Ratio Sequence
|
||||||
|
"classical_ratio_mode": true, // true: a/d false: 1+a-d
|
||||||
|
// Regression
|
||||||
|
//"target_camera_acceptor": "acceptor",
|
||||||
|
"target_camera_acceptor": "",
|
||||||
|
"regressor_cameras_acceptor": [
|
||||||
|
"oxygenation",
|
||||||
|
"volume"
|
||||||
|
],
|
||||||
|
"target_camera_donor": "donor",
|
||||||
|
"regressor_cameras_donor": [
|
||||||
|
// "oxygenation",
|
||||||
|
"volume"
|
||||||
|
],
|
||||||
|
// binning
|
||||||
|
"binning_enable": true,
|
||||||
|
"binning_at_the_end": false,
|
||||||
|
"binning_kernel_size": 4,
|
||||||
|
"binning_stride": 4,
|
||||||
|
"binning_divisor_override": 1,
|
||||||
|
// alignment
|
||||||
|
"alignment_batch_size": 200,
|
||||||
|
"rotation_stabilization_threshold_factor": 3.0, // >= 1.0
|
||||||
|
"rotation_stabilization_threshold_border": 0.9, // <= 1.0
|
||||||
|
// Heart beat detection
|
||||||
|
"lower_freqency_bandpass": 5.0, // Hz
|
||||||
|
"upper_freqency_bandpass": 14.0, // Hz
|
||||||
|
"heartbeat_filtfilt_chuck_size": 10,
|
||||||
|
// Gauss smear
|
||||||
|
"gauss_smear_spatial_width": 8,
|
||||||
|
"gauss_smear_temporal_width": 0.1,
|
||||||
|
"gauss_smear_use_matlab_mask": false,
|
||||||
|
// LED Ramp on
|
||||||
|
"skip_frames_in_the_beginning": 100, // Frames
|
||||||
|
// PyTorch
|
||||||
|
"dtype": "float32",
|
||||||
|
"force_to_cpu": false,
|
||||||
|
// Save
|
||||||
|
"save_as_python": true, // produces .npz files (compressed)
|
||||||
|
"save_as_matlab": false, // produces .hd5 file (compressed)
|
||||||
|
// Save extra information
|
||||||
|
"save_alignment": false,
|
||||||
|
"save_heartbeat": false,
|
||||||
|
"save_factors": false,
|
||||||
|
"save_regression_coefficients": false,
|
||||||
|
// Not important parameter
|
||||||
|
"required_order": [
|
||||||
|
"acceptor",
|
||||||
|
"donor",
|
||||||
|
"oxygenation",
|
||||||
|
"volume"
|
||||||
|
]
|
||||||
|
}
|
||||||
|
|
@ -5,6 +5,8 @@
|
||||||
"raw_path": "raw",
|
"raw_path": "raw",
|
||||||
"export_path": "output_M_Sert_Cre_49",
|
"export_path": "output_M_Sert_Cre_49",
|
||||||
"ref_image_path": "ref_images_M_Sert_Cre_49",
|
"ref_image_path": "ref_images_M_Sert_Cre_49",
|
||||||
|
"heartbeat_remove": false,
|
||||||
|
"gevi": false, // true => gevi, false => geci
|
||||||
// Ratio Sequence
|
// Ratio Sequence
|
||||||
"classical_ratio_mode": true, // true: a/d false: 1+a-d
|
"classical_ratio_mode": true, // true: a/d false: 1+a-d
|
||||||
// Regression
|
// Regression
|
||||||
|
|
|
||||||
71
geci_loader.py
Normal file
71
geci_loader.py
Normal file
|
|
@ -0,0 +1,71 @@
|
||||||
|
import numpy as np
|
||||||
|
import os
|
||||||
|
import argh
|
||||||
|
|
||||||
|
|
||||||
|
# mouse:int = 0, 1, 2, 3, 4
|
||||||
|
def loader(mouse:int = 0, fpath:str = "/data_1/hendrik/gevi") -> None:
|
||||||
|
|
||||||
|
mouse_name = [
|
||||||
|
'M_Sert_Cre_41',
|
||||||
|
'M_Sert_Cre_42',
|
||||||
|
'M_Sert_Cre_45',
|
||||||
|
'M_Sert_Cre_46',
|
||||||
|
'M_Sert_Cre_49'
|
||||||
|
]
|
||||||
|
|
||||||
|
n_tris = [
|
||||||
|
[15, 15, 30, 30, 30, 30,], # 0 in cond 7
|
||||||
|
[15, 15, 30, 30, 30, 30,], # 0 in cond 7
|
||||||
|
[15, 15, 30, 30, 30, 30,], # 0 in cond 7
|
||||||
|
[20, 40, 20, 20,], # 0, 0, 0 in cond 5-7
|
||||||
|
[20, 40, 20, 20,], # 0, 0, 0 in cond 5-7
|
||||||
|
]
|
||||||
|
|
||||||
|
# 41, 42, 45, 46, 49:
|
||||||
|
# "1": "control",
|
||||||
|
# "2": "visual control grating 100 1s",
|
||||||
|
# "3": "optical Stimulation 20Hz 50% 5 Intervals",
|
||||||
|
# "4": "optical Stimulation 20Hz 100% 5 Intervals",
|
||||||
|
# "5": "optical Stimulation 20Hz 50% and grating 100",
|
||||||
|
# "6": "optical Stimulation 20Hz 100% and grating 100",
|
||||||
|
# "7": "grating 3s"
|
||||||
|
|
||||||
|
lbs = [
|
||||||
|
['control', 'visual control', 'op20 50 5', 'op20 100 5', 'op20 50 grat', 'op20 100 grat'],
|
||||||
|
['control', 'visual control', 'op20 50 5', 'op20 100 5', 'op20 50 grat', 'op20 100 grat'],
|
||||||
|
['control', 'visual control', 'op20 50 5', 'op20 100 5', 'op20 50 grat', 'op20 100 grat'],
|
||||||
|
['control', 'visual control', 'op20 50 5', 'op20 100 5'],
|
||||||
|
['control', 'visual control', 'op20 50 5', 'op20 100 5']
|
||||||
|
]
|
||||||
|
|
||||||
|
n_exp = len(n_tris[mouse])
|
||||||
|
|
||||||
|
for i_exp in range(n_exp):
|
||||||
|
n_tri = n_tris[mouse][i_exp]
|
||||||
|
for i_tri in range(n_tri):
|
||||||
|
|
||||||
|
experiment_name: str = f"Exp{i_exp+1:03d}_Trial{i_tri+1:03d}"
|
||||||
|
tmp_fname = os.path.join(
|
||||||
|
fpath, "output_" + mouse_name[mouse], experiment_name + "_acceptor_donor.npz"
|
||||||
|
)
|
||||||
|
print(f'Processing file "{tmp_fname}"...')
|
||||||
|
tmp = np.load(tmp_fname)
|
||||||
|
|
||||||
|
tmp_data_sequence = tmp["data_donor"]
|
||||||
|
tmp_light_signal = tmp["data_acceptor"]
|
||||||
|
|
||||||
|
if (i_exp == 0) and (i_tri == 0):
|
||||||
|
mask = tmp["mask"]
|
||||||
|
new_shape = [n_exp, *tmp_data_sequence.shape]
|
||||||
|
data_sequence = np.zeros(new_shape)
|
||||||
|
light_signal = np.zeros(new_shape)
|
||||||
|
data_sequence[i_exp] += tmp_data_sequence / n_tri
|
||||||
|
light_signal[i_exp] += tmp_light_signal / n_tri
|
||||||
|
|
||||||
|
np.save("dsq_" + mouse_name[mouse], data_sequence)
|
||||||
|
np.save("lsq_" + mouse_name[mouse], light_signal)
|
||||||
|
np.save("msq_" + mouse_name[mouse], mask)
|
||||||
|
|
||||||
|
if __name__ == "__main__":
|
||||||
|
argh.dispatch_command(loader)
|
||||||
209
geci_plot.py
Normal file
209
geci_plot.py
Normal file
|
|
@ -0,0 +1,209 @@
|
||||||
|
import numpy as np
|
||||||
|
import matplotlib.pyplot as plt
|
||||||
|
import h5py # type: ignore
|
||||||
|
import argh
|
||||||
|
import scipy # type: ignore
|
||||||
|
import json
|
||||||
|
import os
|
||||||
|
from jsmin import jsmin # type:ignore
|
||||||
|
from functions.get_trials import get_trials
|
||||||
|
|
||||||
|
|
||||||
|
def func_pow(x, a, b, c):
|
||||||
|
return -a * x**b + c
|
||||||
|
|
||||||
|
|
||||||
|
def func_exp(x, a, b, c):
|
||||||
|
return a * np.exp(-x / b) + c
|
||||||
|
|
||||||
|
|
||||||
|
# mouse: int = 0, 1, 2, 3, 4
|
||||||
|
def plot(
|
||||||
|
filename: str = "config_M_Sert_Cre_49.json",
|
||||||
|
experiment: int = 4,
|
||||||
|
skip_timesteps: int = 100,
|
||||||
|
# If there is no special ROI... Get one! This is just a backup
|
||||||
|
roi_control_path: str = "/data_1/hendrik/2023-03-15/ROI_control.mat",
|
||||||
|
roi_sdraken_path: str = "/data_1/hendrik/2023-03-15/ROI_sDarken.mat",
|
||||||
|
remove_fit: bool = True,
|
||||||
|
fit_power: bool = False, # True => -ax^b ; False => exp(-b)
|
||||||
|
) -> None:
|
||||||
|
|
||||||
|
# lbs = [
|
||||||
|
# [
|
||||||
|
# "control",
|
||||||
|
# "visual control",
|
||||||
|
# "op20 50 5",
|
||||||
|
# "op20 100 5",
|
||||||
|
# "op20 50 grat",
|
||||||
|
# "op20 100 grat",
|
||||||
|
# ],
|
||||||
|
# [
|
||||||
|
# "control",
|
||||||
|
# "visual control",
|
||||||
|
# "op20 50 5",
|
||||||
|
# "op20 100 5",
|
||||||
|
# "op20 50 grat",
|
||||||
|
# "op20 100 grat",
|
||||||
|
# ],
|
||||||
|
# [
|
||||||
|
# "control",
|
||||||
|
# "visual control",
|
||||||
|
# "op20 50 5",
|
||||||
|
# "op20 100 5",
|
||||||
|
# "op20 50 grat",
|
||||||
|
# "op20 100 grat",
|
||||||
|
# ],
|
||||||
|
# ["control", "visual control", "op20 50 5", "op20 100 5"],
|
||||||
|
# ["control", "visual control", "op20 50 5", "op20 100 5"],
|
||||||
|
# ]
|
||||||
|
|
||||||
|
if os.path.isfile(filename) is False:
|
||||||
|
print(f"{filename} is missing")
|
||||||
|
exit()
|
||||||
|
|
||||||
|
with open(filename, "r") as file:
|
||||||
|
config = json.loads(jsmin(file.read()))
|
||||||
|
|
||||||
|
raw_data_path: str = os.path.join(
|
||||||
|
config["basic_path"],
|
||||||
|
config["recoding_data"],
|
||||||
|
config["mouse_identifier"],
|
||||||
|
config["raw_path"],
|
||||||
|
)
|
||||||
|
|
||||||
|
if os.path.isdir(raw_data_path) is False:
|
||||||
|
print(f"ERROR: could not find raw directory {raw_data_path}!!!!")
|
||||||
|
exit()
|
||||||
|
|
||||||
|
trials = get_trials(raw_data_path, experiment).numpy()
|
||||||
|
assert trials.shape[0] > 0
|
||||||
|
|
||||||
|
with open(os.path.join(raw_data_path, f"Exp{experiment:03d}_Trial{trials[0]:03d}_Part001_meta.txt"), "r") as file:
|
||||||
|
metadata = json.loads(jsmin(file.read()))
|
||||||
|
|
||||||
|
experiment_names = metadata['sessionMetaData']['experimentNames'][str(experiment)]
|
||||||
|
|
||||||
|
roi_path: str = os.path.join(config["basic_path"], config["recoding_data"])
|
||||||
|
roi_control_mat: str = os.path.join(roi_path, "ROI_control.mat")
|
||||||
|
roi_sdarken_mat: str = os.path.join(roi_path, "ROI_sDarken.mat")
|
||||||
|
|
||||||
|
if os.path.isdir(roi_control_mat):
|
||||||
|
roi_control_path = roi_control_mat
|
||||||
|
|
||||||
|
if os.path.isdir(roi_sdarken_mat):
|
||||||
|
roi_sdraken_path = roi_sdarken_mat
|
||||||
|
|
||||||
|
print("Load data...")
|
||||||
|
data = np.load("dsq_" + config["mouse_identifier"] + ".npy", mmap_mode="r")
|
||||||
|
|
||||||
|
print("Load light signal...")
|
||||||
|
light = np.load("lsq_" + config["mouse_identifier"] + ".npy", mmap_mode="r")
|
||||||
|
|
||||||
|
print("Load mask...")
|
||||||
|
mask = np.load("msq_" + config["mouse_identifier"] + ".npy")
|
||||||
|
|
||||||
|
hf = h5py.File(roi_control_path, "r")
|
||||||
|
roi_lighten = np.array(hf["roi"]).T
|
||||||
|
roi_lighten *= mask
|
||||||
|
|
||||||
|
hf = h5py.File(roi_sdraken_path, "r")
|
||||||
|
roi_darken = np.array(hf["roi"]).T
|
||||||
|
roi_darken *= mask
|
||||||
|
|
||||||
|
plt.figure(1)
|
||||||
|
a_show = data[experiment - 1, :, :, 1000].copy()
|
||||||
|
a_show[(roi_darken + roi_lighten) < 0.5] = np.nan
|
||||||
|
plt.imshow(a_show)
|
||||||
|
plt.title(f"{config["mouse_identifier"]} -- Experiment: {experiment}")
|
||||||
|
plt.show(block=False)
|
||||||
|
|
||||||
|
plt.figure(2)
|
||||||
|
a_dontshow = data[experiment - 1, :, :, 1000].copy()
|
||||||
|
a_dontshow[(roi_darken + roi_lighten) > 0.5] = np.nan
|
||||||
|
plt.imshow(a_dontshow)
|
||||||
|
plt.title(f"{config["mouse_identifier"]} -- Experiment: {experiment}")
|
||||||
|
plt.show(block=False)
|
||||||
|
|
||||||
|
plt.figure(3)
|
||||||
|
light_exp = light[experiment - 1, :, :, skip_timesteps:].copy()
|
||||||
|
light_exp[(roi_darken + roi_lighten) < 0.5, :] = 0.0
|
||||||
|
light_signal = light_exp.mean(axis=(0, 1))
|
||||||
|
light_signal -= light_signal.min()
|
||||||
|
light_signal /= light_signal.max()
|
||||||
|
|
||||||
|
a_exp = data[experiment - 1, :, :, skip_timesteps:].copy()
|
||||||
|
|
||||||
|
if remove_fit:
|
||||||
|
combined_matrix = (roi_darken + roi_lighten) > 0
|
||||||
|
idx = np.where(combined_matrix)
|
||||||
|
for idx_pos in range(0, idx[0].shape[0]):
|
||||||
|
temp = a_exp[idx[0][idx_pos], idx[1][idx_pos], :]
|
||||||
|
temp -= temp.mean()
|
||||||
|
|
||||||
|
data_time = np.arange(0, temp.shape[0], dtype=np.float32) + skip_timesteps
|
||||||
|
data_time /= 100.0
|
||||||
|
|
||||||
|
data_min = temp.min()
|
||||||
|
data_max = temp.max()
|
||||||
|
data_delta = data_max - data_min
|
||||||
|
a_min = data_min - data_delta
|
||||||
|
b_min = 0.01
|
||||||
|
a_max = data_max + data_delta
|
||||||
|
if fit_power:
|
||||||
|
b_max = 10.0
|
||||||
|
else:
|
||||||
|
b_max = 100.0
|
||||||
|
c_min = data_min - data_delta
|
||||||
|
c_max = data_max + data_delta
|
||||||
|
|
||||||
|
try:
|
||||||
|
if fit_power:
|
||||||
|
popt, _ = scipy.optimize.curve_fit(
|
||||||
|
f=func_pow,
|
||||||
|
xdata=data_time,
|
||||||
|
ydata=np.nan_to_num(temp),
|
||||||
|
bounds=([a_min, b_min, c_min], [a_max, b_max, c_max]),
|
||||||
|
)
|
||||||
|
pattern: np.ndarray | None = func_pow(data_time, *popt)
|
||||||
|
else:
|
||||||
|
popt, _ = scipy.optimize.curve_fit(
|
||||||
|
f=func_exp,
|
||||||
|
xdata=data_time,
|
||||||
|
ydata=np.nan_to_num(temp),
|
||||||
|
bounds=([a_min, b_min, c_min], [a_max, b_max, c_max]),
|
||||||
|
)
|
||||||
|
pattern = func_exp(data_time, *popt)
|
||||||
|
|
||||||
|
assert pattern is not None
|
||||||
|
pattern -= pattern.mean()
|
||||||
|
|
||||||
|
scale = (temp * pattern).sum() / (pattern**2).sum()
|
||||||
|
pattern *= scale
|
||||||
|
|
||||||
|
except ValueError:
|
||||||
|
print(f"Fit failed: Position ({idx[0][idx_pos]}, {idx[1][idx_pos]}")
|
||||||
|
pattern = None
|
||||||
|
|
||||||
|
if pattern is not None:
|
||||||
|
temp -= pattern
|
||||||
|
|
||||||
|
darken = a_exp[roi_darken > 0.5, :].sum(axis=0) / (roi_darken > 0.5).sum()
|
||||||
|
lighten = a_exp[roi_lighten > 0.5, :].sum(axis=0) / (roi_lighten > 0.5).sum()
|
||||||
|
|
||||||
|
light_signal *= darken.max() - darken.min()
|
||||||
|
light_signal += darken.min()
|
||||||
|
|
||||||
|
time_axis = np.arange(0, lighten.shape[-1], dtype=np.float32) + skip_timesteps
|
||||||
|
time_axis /= 100.0
|
||||||
|
|
||||||
|
plt.plot(time_axis, light_signal, c="k", label="light")
|
||||||
|
plt.plot(time_axis, darken, label="sDarken")
|
||||||
|
plt.plot(time_axis, lighten, label="control")
|
||||||
|
plt.title(f"{config["mouse_identifier"]} -- Experiment: {experiment} ({experiment_names})")
|
||||||
|
plt.legend()
|
||||||
|
plt.show()
|
||||||
|
|
||||||
|
|
||||||
|
if __name__ == "__main__":
|
||||||
|
argh.dispatch_command(plot)
|
||||||
|
|
@ -99,7 +99,7 @@ def process_trial(
|
||||||
free_mem = cuda_total_memory - max(
|
free_mem = cuda_total_memory - max(
|
||||||
[torch.cuda.memory_reserved(device), torch.cuda.memory_allocated(device)]
|
[torch.cuda.memory_reserved(device), torch.cuda.memory_allocated(device)]
|
||||||
)
|
)
|
||||||
mylogger.info(f"CUDA memory: {free_mem//1024} MByte")
|
mylogger.info(f"CUDA memory: {free_mem // 1024} MByte")
|
||||||
|
|
||||||
mylogger.info(f"Data shape: {data.shape}")
|
mylogger.info(f"Data shape: {data.shape}")
|
||||||
mylogger.info("-==- Done -==-")
|
mylogger.info("-==- Done -==-")
|
||||||
|
|
@ -266,9 +266,9 @@ def process_trial(
|
||||||
batch_size=config["alignment_batch_size"],
|
batch_size=config["alignment_batch_size"],
|
||||||
fill_value=-100.0,
|
fill_value=-100.0,
|
||||||
)
|
)
|
||||||
mylogger.info(f"Rotation: {round(float(angle_refref[0]),2)} degree")
|
mylogger.info(f"Rotation: {round(float(angle_refref[0]), 2)} degree")
|
||||||
mylogger.info(
|
mylogger.info(
|
||||||
f"Translation: {round(float(tvec_refref[0]),1)} x {round(float(tvec_refref[1]),1)} pixel"
|
f"Translation: {round(float(tvec_refref[0]), 1)} x {round(float(tvec_refref[1]), 1)} pixel"
|
||||||
)
|
)
|
||||||
|
|
||||||
if config["save_alignment"]:
|
if config["save_alignment"]:
|
||||||
|
|
@ -285,7 +285,7 @@ def process_trial(
|
||||||
np.save(temp_path, tvec_refref.cpu())
|
np.save(temp_path, tvec_refref.cpu())
|
||||||
|
|
||||||
mylogger.info("Moving & rotating the oxygenation ref image")
|
mylogger.info("Moving & rotating the oxygenation ref image")
|
||||||
ref_image_oxygenation = tv.transforms.functional.affine(
|
ref_image_oxygenation = tv.transforms.functional.affine( # type: ignore
|
||||||
img=ref_image_oxygenation.unsqueeze(0),
|
img=ref_image_oxygenation.unsqueeze(0),
|
||||||
angle=-float(angle_refref),
|
angle=-float(angle_refref),
|
||||||
translate=[0, 0],
|
translate=[0, 0],
|
||||||
|
|
@ -295,7 +295,7 @@ def process_trial(
|
||||||
fill=-100.0,
|
fill=-100.0,
|
||||||
)
|
)
|
||||||
|
|
||||||
ref_image_oxygenation = tv.transforms.functional.affine(
|
ref_image_oxygenation = tv.transforms.functional.affine( # type: ignore
|
||||||
img=ref_image_oxygenation,
|
img=ref_image_oxygenation,
|
||||||
angle=0,
|
angle=0,
|
||||||
translate=[tvec_refref[1], tvec_refref[0]],
|
translate=[tvec_refref[1], tvec_refref[0]],
|
||||||
|
|
@ -313,8 +313,8 @@ def process_trial(
|
||||||
volume_index: int = config["required_order"].index("volume")
|
volume_index: int = config["required_order"].index("volume")
|
||||||
|
|
||||||
mylogger.info("Rotate acceptor")
|
mylogger.info("Rotate acceptor")
|
||||||
data[acceptor_index, ...] = tv.transforms.functional.affine(
|
data[acceptor_index, ...] = tv.transforms.functional.affine( # type: ignore
|
||||||
img=data[acceptor_index, ...],
|
img=data[acceptor_index, ...], # type: ignore
|
||||||
angle=-float(angle_refref),
|
angle=-float(angle_refref),
|
||||||
translate=[0, 0],
|
translate=[0, 0],
|
||||||
scale=1.0,
|
scale=1.0,
|
||||||
|
|
@ -324,7 +324,7 @@ def process_trial(
|
||||||
)
|
)
|
||||||
|
|
||||||
mylogger.info("Translate acceptor")
|
mylogger.info("Translate acceptor")
|
||||||
data[acceptor_index, ...] = tv.transforms.functional.affine(
|
data[acceptor_index, ...] = tv.transforms.functional.affine( # type: ignore
|
||||||
img=data[acceptor_index, ...],
|
img=data[acceptor_index, ...],
|
||||||
angle=0,
|
angle=0,
|
||||||
translate=[tvec_refref[1], tvec_refref[0]],
|
translate=[tvec_refref[1], tvec_refref[0]],
|
||||||
|
|
@ -335,7 +335,7 @@ def process_trial(
|
||||||
)
|
)
|
||||||
|
|
||||||
mylogger.info("Rotate oxygenation")
|
mylogger.info("Rotate oxygenation")
|
||||||
data[oxygenation_index, ...] = tv.transforms.functional.affine(
|
data[oxygenation_index, ...] = tv.transforms.functional.affine( # type: ignore
|
||||||
img=data[oxygenation_index, ...],
|
img=data[oxygenation_index, ...],
|
||||||
angle=-float(angle_refref),
|
angle=-float(angle_refref),
|
||||||
translate=[0, 0],
|
translate=[0, 0],
|
||||||
|
|
@ -346,7 +346,7 @@ def process_trial(
|
||||||
)
|
)
|
||||||
|
|
||||||
mylogger.info("Translate oxygenation")
|
mylogger.info("Translate oxygenation")
|
||||||
data[oxygenation_index, ...] = tv.transforms.functional.affine(
|
data[oxygenation_index, ...] = tv.transforms.functional.affine( # type: ignore
|
||||||
img=data[oxygenation_index, ...],
|
img=data[oxygenation_index, ...],
|
||||||
angle=0,
|
angle=0,
|
||||||
translate=[tvec_refref[1], tvec_refref[0]],
|
translate=[tvec_refref[1], tvec_refref[0]],
|
||||||
|
|
@ -359,7 +359,7 @@ def process_trial(
|
||||||
|
|
||||||
mylogger.info("Perform rotation between donor and volume and its ref images")
|
mylogger.info("Perform rotation between donor and volume and its ref images")
|
||||||
mylogger.info("for all frames and then rotate all the data accordingly")
|
mylogger.info("for all frames and then rotate all the data accordingly")
|
||||||
perform_donor_volume_rotation
|
|
||||||
(
|
(
|
||||||
data[acceptor_index, ...],
|
data[acceptor_index, ...],
|
||||||
data[donor_index, ...],
|
data[donor_index, ...],
|
||||||
|
|
@ -381,9 +381,9 @@ def process_trial(
|
||||||
|
|
||||||
mylogger.info(
|
mylogger.info(
|
||||||
f"angles: "
|
f"angles: "
|
||||||
f"min {round(float(angle_donor_volume.min()),2)} "
|
f"min {round(float(angle_donor_volume.min()), 2)} "
|
||||||
f"max {round(float(angle_donor_volume.max()),2)} "
|
f"max {round(float(angle_donor_volume.max()), 2)} "
|
||||||
f"mean {round(float(angle_donor_volume.mean()),2)} "
|
f"mean {round(float(angle_donor_volume.mean()), 2)} "
|
||||||
)
|
)
|
||||||
|
|
||||||
if config["save_alignment"]:
|
if config["save_alignment"]:
|
||||||
|
|
@ -417,15 +417,15 @@ def process_trial(
|
||||||
|
|
||||||
mylogger.info(
|
mylogger.info(
|
||||||
f"translation dim 0: "
|
f"translation dim 0: "
|
||||||
f"min {round(float(tvec_donor_volume[:,0].min()),1)} "
|
f"min {round(float(tvec_donor_volume[:, 0].min()), 1)} "
|
||||||
f"max {round(float(tvec_donor_volume[:,0].max()),1)} "
|
f"max {round(float(tvec_donor_volume[:, 0].max()), 1)} "
|
||||||
f"mean {round(float(tvec_donor_volume[:,0].mean()),1)} "
|
f"mean {round(float(tvec_donor_volume[:, 0].mean()), 1)} "
|
||||||
)
|
)
|
||||||
mylogger.info(
|
mylogger.info(
|
||||||
f"translation dim 1: "
|
f"translation dim 1: "
|
||||||
f"min {round(float(tvec_donor_volume[:,1].min()),1)} "
|
f"min {round(float(tvec_donor_volume[:, 1].min()), 1)} "
|
||||||
f"max {round(float(tvec_donor_volume[:,1].max()),1)} "
|
f"max {round(float(tvec_donor_volume[:, 1].max()), 1)} "
|
||||||
f"mean {round(float(tvec_donor_volume[:,1].mean()),1)} "
|
f"mean {round(float(tvec_donor_volume[:, 1].mean()), 1)} "
|
||||||
)
|
)
|
||||||
|
|
||||||
if config["save_alignment"]:
|
if config["save_alignment"]:
|
||||||
|
|
@ -471,172 +471,183 @@ def process_trial(
|
||||||
|
|
||||||
sample_frequency: float = 1.0 / meta_frame_time
|
sample_frequency: float = 1.0 / meta_frame_time
|
||||||
|
|
||||||
mylogger.info("Extract heartbeat from volume signal")
|
if config["gevi"]:
|
||||||
heartbeat_ts: torch.Tensor = bandpass(
|
assert config["heartbeat_remove"]
|
||||||
data=data[volume_index, ...].movedim(0, -1).clone(),
|
|
||||||
low_frequency=config["lower_freqency_bandpass"],
|
|
||||||
high_frequency=config["upper_freqency_bandpass"],
|
|
||||||
fs=sample_frequency,
|
|
||||||
filtfilt_chuck_size=config["heartbeat_filtfilt_chuck_size"],
|
|
||||||
)
|
|
||||||
heartbeat_ts = heartbeat_ts.flatten(start_dim=0, end_dim=-2)
|
|
||||||
mask_flatten: torch.Tensor = mask_positve.flatten(start_dim=0, end_dim=-1)
|
|
||||||
|
|
||||||
heartbeat_ts = heartbeat_ts[mask_flatten, :]
|
if config["heartbeat_remove"]:
|
||||||
|
mylogger.info("Extract heartbeat from volume signal")
|
||||||
heartbeat_ts = heartbeat_ts.movedim(0, -1)
|
heartbeat_ts: torch.Tensor = bandpass(
|
||||||
heartbeat_ts -= heartbeat_ts.mean(dim=0, keepdim=True)
|
data=data[volume_index, ...].movedim(0, -1).clone(),
|
||||||
|
|
||||||
try:
|
|
||||||
volume_heartbeat, _, _ = torch.linalg.svd(heartbeat_ts, full_matrices=False)
|
|
||||||
except torch.cuda.OutOfMemoryError:
|
|
||||||
mylogger.info("torch.cuda.OutOfMemoryError: Fallback to cpu")
|
|
||||||
volume_heartbeat_cpu, _, _ = torch.linalg.svd(
|
|
||||||
heartbeat_ts.cpu(), full_matrices=False
|
|
||||||
)
|
|
||||||
volume_heartbeat = volume_heartbeat_cpu.to(heartbeat_ts.data, copy=True)
|
|
||||||
del volume_heartbeat_cpu
|
|
||||||
|
|
||||||
volume_heartbeat = volume_heartbeat[:, 0]
|
|
||||||
volume_heartbeat -= volume_heartbeat[
|
|
||||||
config["skip_frames_in_the_beginning"] :
|
|
||||||
].mean()
|
|
||||||
|
|
||||||
del heartbeat_ts
|
|
||||||
|
|
||||||
if device != torch.device("cpu"):
|
|
||||||
torch.cuda.empty_cache()
|
|
||||||
mylogger.info("Empty CUDA cache")
|
|
||||||
free_mem = cuda_total_memory - max(
|
|
||||||
[torch.cuda.memory_reserved(device), torch.cuda.memory_allocated(device)]
|
|
||||||
)
|
|
||||||
mylogger.info(f"CUDA memory: {free_mem//1024} MByte")
|
|
||||||
|
|
||||||
if config["save_heartbeat"]:
|
|
||||||
temp_path = os.path.join(
|
|
||||||
config["export_path"], experiment_name + "_volume_heartbeat.npy"
|
|
||||||
)
|
|
||||||
mylogger.info(f"Save volume heartbeat to {temp_path}")
|
|
||||||
np.save(temp_path, volume_heartbeat.cpu())
|
|
||||||
mylogger.info("-==- Done -==-")
|
|
||||||
|
|
||||||
volume_heartbeat = volume_heartbeat.unsqueeze(0).unsqueeze(0)
|
|
||||||
norm_volume_heartbeat = (
|
|
||||||
volume_heartbeat[..., config["skip_frames_in_the_beginning"] :] ** 2
|
|
||||||
).sum(dim=-1)
|
|
||||||
|
|
||||||
heartbeat_coefficients: torch.Tensor = torch.zeros(
|
|
||||||
(data.shape[0], data.shape[-2], data.shape[-1]),
|
|
||||||
dtype=data.dtype,
|
|
||||||
device=data.device,
|
|
||||||
)
|
|
||||||
for i in range(0, data.shape[0]):
|
|
||||||
y = bandpass(
|
|
||||||
data=data[i, ...].movedim(0, -1).clone(),
|
|
||||||
low_frequency=config["lower_freqency_bandpass"],
|
low_frequency=config["lower_freqency_bandpass"],
|
||||||
high_frequency=config["upper_freqency_bandpass"],
|
high_frequency=config["upper_freqency_bandpass"],
|
||||||
fs=sample_frequency,
|
fs=sample_frequency,
|
||||||
filtfilt_chuck_size=config["heartbeat_filtfilt_chuck_size"],
|
filtfilt_chuck_size=config["heartbeat_filtfilt_chuck_size"],
|
||||||
)[..., config["skip_frames_in_the_beginning"] :]
|
|
||||||
y -= y.mean(dim=-1, keepdim=True)
|
|
||||||
|
|
||||||
heartbeat_coefficients[i, ...] = (
|
|
||||||
volume_heartbeat[..., config["skip_frames_in_the_beginning"] :] * y
|
|
||||||
).sum(dim=-1) / norm_volume_heartbeat
|
|
||||||
|
|
||||||
heartbeat_coefficients[i, ...] *= mask_positve.type(
|
|
||||||
dtype=heartbeat_coefficients.dtype
|
|
||||||
)
|
)
|
||||||
del y
|
heartbeat_ts = heartbeat_ts.flatten(start_dim=0, end_dim=-2)
|
||||||
|
mask_flatten: torch.Tensor = mask_positve.flatten(start_dim=0, end_dim=-1)
|
||||||
|
|
||||||
if config["save_heartbeat"]:
|
heartbeat_ts = heartbeat_ts[mask_flatten, :]
|
||||||
temp_path = os.path.join(
|
|
||||||
config["export_path"], experiment_name + "_heartbeat_coefficients.npy"
|
heartbeat_ts = heartbeat_ts.movedim(0, -1)
|
||||||
|
heartbeat_ts -= heartbeat_ts.mean(dim=0, keepdim=True)
|
||||||
|
|
||||||
|
try:
|
||||||
|
volume_heartbeat, _, _ = torch.linalg.svd(heartbeat_ts, full_matrices=False)
|
||||||
|
except torch.cuda.OutOfMemoryError:
|
||||||
|
mylogger.info("torch.cuda.OutOfMemoryError: Fallback to cpu")
|
||||||
|
volume_heartbeat_cpu, _, _ = torch.linalg.svd(
|
||||||
|
heartbeat_ts.cpu(), full_matrices=False
|
||||||
|
)
|
||||||
|
volume_heartbeat = volume_heartbeat_cpu.to(heartbeat_ts.data, copy=True)
|
||||||
|
del volume_heartbeat_cpu
|
||||||
|
|
||||||
|
volume_heartbeat = volume_heartbeat[:, 0]
|
||||||
|
volume_heartbeat -= volume_heartbeat[
|
||||||
|
config["skip_frames_in_the_beginning"] :
|
||||||
|
].mean()
|
||||||
|
|
||||||
|
del heartbeat_ts
|
||||||
|
|
||||||
|
if device != torch.device("cpu"):
|
||||||
|
torch.cuda.empty_cache()
|
||||||
|
mylogger.info("Empty CUDA cache")
|
||||||
|
free_mem = cuda_total_memory - max(
|
||||||
|
[
|
||||||
|
torch.cuda.memory_reserved(device),
|
||||||
|
torch.cuda.memory_allocated(device),
|
||||||
|
]
|
||||||
|
)
|
||||||
|
mylogger.info(f"CUDA memory: {free_mem // 1024} MByte")
|
||||||
|
|
||||||
|
if config["save_heartbeat"]:
|
||||||
|
temp_path = os.path.join(
|
||||||
|
config["export_path"], experiment_name + "_volume_heartbeat.npy"
|
||||||
|
)
|
||||||
|
mylogger.info(f"Save volume heartbeat to {temp_path}")
|
||||||
|
np.save(temp_path, volume_heartbeat.cpu())
|
||||||
|
mylogger.info("-==- Done -==-")
|
||||||
|
|
||||||
|
volume_heartbeat = volume_heartbeat.unsqueeze(0).unsqueeze(0)
|
||||||
|
norm_volume_heartbeat = (
|
||||||
|
volume_heartbeat[..., config["skip_frames_in_the_beginning"] :] ** 2
|
||||||
|
).sum(dim=-1)
|
||||||
|
|
||||||
|
heartbeat_coefficients: torch.Tensor = torch.zeros(
|
||||||
|
(data.shape[0], data.shape[-2], data.shape[-1]),
|
||||||
|
dtype=data.dtype,
|
||||||
|
device=data.device,
|
||||||
)
|
)
|
||||||
mylogger.info(f"Save heartbeat coefficients to {temp_path}")
|
for i in range(0, data.shape[0]):
|
||||||
np.save(temp_path, heartbeat_coefficients.cpu())
|
y = bandpass(
|
||||||
mylogger.info("-==- Done -==-")
|
data=data[i, ...].movedim(0, -1).clone(),
|
||||||
|
low_frequency=config["lower_freqency_bandpass"],
|
||||||
|
high_frequency=config["upper_freqency_bandpass"],
|
||||||
|
fs=sample_frequency,
|
||||||
|
filtfilt_chuck_size=config["heartbeat_filtfilt_chuck_size"],
|
||||||
|
)[..., config["skip_frames_in_the_beginning"] :]
|
||||||
|
y -= y.mean(dim=-1, keepdim=True)
|
||||||
|
|
||||||
mylogger.info("Remove heart beat from data")
|
heartbeat_coefficients[i, ...] = (
|
||||||
data -= heartbeat_coefficients.unsqueeze(1) * volume_heartbeat.unsqueeze(0).movedim(
|
volume_heartbeat[..., config["skip_frames_in_the_beginning"] :] * y
|
||||||
-1, 1
|
).sum(dim=-1) / norm_volume_heartbeat
|
||||||
)
|
|
||||||
mylogger.info("-==- Done -==-")
|
|
||||||
|
|
||||||
donor_heartbeat_factor = heartbeat_coefficients[donor_index, ...].clone()
|
heartbeat_coefficients[i, ...] *= mask_positve.type(
|
||||||
acceptor_heartbeat_factor = heartbeat_coefficients[acceptor_index, ...].clone()
|
dtype=heartbeat_coefficients.dtype
|
||||||
del heartbeat_coefficients
|
)
|
||||||
|
del y
|
||||||
|
|
||||||
if device != torch.device("cpu"):
|
if config["save_heartbeat"]:
|
||||||
torch.cuda.empty_cache()
|
temp_path = os.path.join(
|
||||||
mylogger.info("Empty CUDA cache")
|
config["export_path"], experiment_name + "_heartbeat_coefficients.npy"
|
||||||
free_mem = cuda_total_memory - max(
|
)
|
||||||
[torch.cuda.memory_reserved(device), torch.cuda.memory_allocated(device)]
|
mylogger.info(f"Save heartbeat coefficients to {temp_path}")
|
||||||
|
np.save(temp_path, heartbeat_coefficients.cpu())
|
||||||
|
mylogger.info("-==- Done -==-")
|
||||||
|
|
||||||
|
mylogger.info("Remove heart beat from data")
|
||||||
|
data -= heartbeat_coefficients.unsqueeze(1) * volume_heartbeat.unsqueeze(
|
||||||
|
0
|
||||||
|
).movedim(-1, 1)
|
||||||
|
mylogger.info("-==- Done -==-")
|
||||||
|
|
||||||
|
donor_heartbeat_factor = heartbeat_coefficients[donor_index, ...].clone()
|
||||||
|
acceptor_heartbeat_factor = heartbeat_coefficients[acceptor_index, ...].clone()
|
||||||
|
del heartbeat_coefficients
|
||||||
|
|
||||||
|
if device != torch.device("cpu"):
|
||||||
|
torch.cuda.empty_cache()
|
||||||
|
mylogger.info("Empty CUDA cache")
|
||||||
|
free_mem = cuda_total_memory - max(
|
||||||
|
[
|
||||||
|
torch.cuda.memory_reserved(device),
|
||||||
|
torch.cuda.memory_allocated(device),
|
||||||
|
]
|
||||||
|
)
|
||||||
|
mylogger.info(f"CUDA memory: {free_mem // 1024} MByte")
|
||||||
|
|
||||||
|
mylogger.info("Calculate scaling factor for donor and acceptor")
|
||||||
|
donor_factor: torch.Tensor = (
|
||||||
|
donor_heartbeat_factor + acceptor_heartbeat_factor
|
||||||
|
) / (2 * donor_heartbeat_factor)
|
||||||
|
acceptor_factor: torch.Tensor = (
|
||||||
|
donor_heartbeat_factor + acceptor_heartbeat_factor
|
||||||
|
) / (2 * acceptor_heartbeat_factor)
|
||||||
|
|
||||||
|
del donor_heartbeat_factor
|
||||||
|
del acceptor_heartbeat_factor
|
||||||
|
|
||||||
|
if config["save_factors"]:
|
||||||
|
temp_path = os.path.join(
|
||||||
|
config["export_path"], experiment_name + "_donor_factor.npy"
|
||||||
|
)
|
||||||
|
mylogger.info(f"Save donor factor to {temp_path}")
|
||||||
|
np.save(temp_path, donor_factor.cpu())
|
||||||
|
|
||||||
|
temp_path = os.path.join(
|
||||||
|
config["export_path"], experiment_name + "_acceptor_factor.npy"
|
||||||
|
)
|
||||||
|
mylogger.info(f"Save acceptor factor to {temp_path}")
|
||||||
|
np.save(temp_path, acceptor_factor.cpu())
|
||||||
|
mylogger.info("-==- Done -==-")
|
||||||
|
|
||||||
|
mylogger.info("Scale acceptor to heart beat amplitude")
|
||||||
|
mylogger.info("Calculate mean")
|
||||||
|
mean_values_acceptor = data[
|
||||||
|
acceptor_index, config["skip_frames_in_the_beginning"] :, ...
|
||||||
|
].nanmean(dim=0, keepdim=True)
|
||||||
|
|
||||||
|
mylogger.info("Remove mean")
|
||||||
|
data[acceptor_index, ...] -= mean_values_acceptor
|
||||||
|
mylogger.info("Apply acceptor_factor and mask")
|
||||||
|
data[acceptor_index, ...] *= acceptor_factor.unsqueeze(
|
||||||
|
0
|
||||||
|
) * mask_positve.unsqueeze(0)
|
||||||
|
mylogger.info("Add mean")
|
||||||
|
data[acceptor_index, ...] += mean_values_acceptor
|
||||||
|
mylogger.info("-==- Done -==-")
|
||||||
|
|
||||||
|
mylogger.info("Scale donor to heart beat amplitude")
|
||||||
|
mylogger.info("Calculate mean")
|
||||||
|
mean_values_donor = data[
|
||||||
|
donor_index, config["skip_frames_in_the_beginning"] :, ...
|
||||||
|
].nanmean(dim=0, keepdim=True)
|
||||||
|
mylogger.info("Remove mean")
|
||||||
|
data[donor_index, ...] -= mean_values_donor
|
||||||
|
mylogger.info("Apply donor_factor and mask")
|
||||||
|
data[donor_index, ...] *= donor_factor.unsqueeze(0) * mask_positve.unsqueeze(0)
|
||||||
|
mylogger.info("Add mean")
|
||||||
|
data[donor_index, ...] += mean_values_donor
|
||||||
|
mylogger.info("-==- Done -==-")
|
||||||
|
|
||||||
|
mylogger.info("Divide by mean over time")
|
||||||
|
data /= data[:, config["skip_frames_in_the_beginning"] :, ...].nanmean(
|
||||||
|
dim=1,
|
||||||
|
keepdim=True,
|
||||||
)
|
)
|
||||||
mylogger.info(f"CUDA memory: {free_mem//1024} MByte")
|
|
||||||
|
|
||||||
mylogger.info("Calculate scaling factor for donor and acceptor")
|
mylogger.info("-==- Done -==-")
|
||||||
donor_factor: torch.Tensor = (
|
|
||||||
donor_heartbeat_factor + acceptor_heartbeat_factor
|
|
||||||
) / (2 * donor_heartbeat_factor)
|
|
||||||
acceptor_factor: torch.Tensor = (
|
|
||||||
donor_heartbeat_factor + acceptor_heartbeat_factor
|
|
||||||
) / (2 * acceptor_heartbeat_factor)
|
|
||||||
|
|
||||||
del donor_heartbeat_factor
|
|
||||||
del acceptor_heartbeat_factor
|
|
||||||
|
|
||||||
if config["save_factors"]:
|
|
||||||
temp_path = os.path.join(
|
|
||||||
config["export_path"], experiment_name + "_donor_factor.npy"
|
|
||||||
)
|
|
||||||
mylogger.info(f"Save donor factor to {temp_path}")
|
|
||||||
np.save(temp_path, donor_factor.cpu())
|
|
||||||
|
|
||||||
temp_path = os.path.join(
|
|
||||||
config["export_path"], experiment_name + "_acceptor_factor.npy"
|
|
||||||
)
|
|
||||||
mylogger.info(f"Save acceptor factor to {temp_path}")
|
|
||||||
np.save(temp_path, acceptor_factor.cpu())
|
|
||||||
mylogger.info("-==- Done -==-")
|
|
||||||
|
|
||||||
mylogger.info("Scale acceptor to heart beat amplitude")
|
|
||||||
mylogger.info("Calculate mean")
|
|
||||||
mean_values_acceptor = data[
|
|
||||||
acceptor_index, config["skip_frames_in_the_beginning"] :, ...
|
|
||||||
].nanmean(dim=0, keepdim=True)
|
|
||||||
|
|
||||||
mylogger.info("Remove mean")
|
|
||||||
data[acceptor_index, ...] -= mean_values_acceptor
|
|
||||||
mylogger.info("Apply acceptor_factor and mask")
|
|
||||||
data[acceptor_index, ...] *= acceptor_factor.unsqueeze(0) * mask_positve.unsqueeze(
|
|
||||||
0
|
|
||||||
)
|
|
||||||
mylogger.info("Add mean")
|
|
||||||
data[acceptor_index, ...] += mean_values_acceptor
|
|
||||||
mylogger.info("-==- Done -==-")
|
|
||||||
|
|
||||||
mylogger.info("Scale donor to heart beat amplitude")
|
|
||||||
mylogger.info("Calculate mean")
|
|
||||||
mean_values_donor = data[
|
|
||||||
donor_index, config["skip_frames_in_the_beginning"] :, ...
|
|
||||||
].nanmean(dim=0, keepdim=True)
|
|
||||||
mylogger.info("Remove mean")
|
|
||||||
data[donor_index, ...] -= mean_values_donor
|
|
||||||
mylogger.info("Apply donor_factor and mask")
|
|
||||||
data[donor_index, ...] *= donor_factor.unsqueeze(0) * mask_positve.unsqueeze(0)
|
|
||||||
mylogger.info("Add mean")
|
|
||||||
data[donor_index, ...] += mean_values_donor
|
|
||||||
mylogger.info("-==- Done -==-")
|
|
||||||
|
|
||||||
mylogger.info("Divide by mean over time")
|
|
||||||
data /= data[:, config["skip_frames_in_the_beginning"] :, ...].nanmean(
|
|
||||||
dim=1,
|
|
||||||
keepdim=True,
|
|
||||||
)
|
|
||||||
data = data.nan_to_num(nan=0.0)
|
data = data.nan_to_num(nan=0.0)
|
||||||
mylogger.info("-==- Done -==-")
|
|
||||||
|
|
||||||
mylogger.info("Preparation for regression -- Gauss smear")
|
mylogger.info("Preparation for regression -- Gauss smear")
|
||||||
spatial_width = float(config["gauss_smear_spatial_width"])
|
spatial_width = float(config["gauss_smear_spatial_width"])
|
||||||
|
|
||||||
|
|
@ -669,7 +680,7 @@ def process_trial(
|
||||||
free_mem = cuda_total_memory - max(
|
free_mem = cuda_total_memory - max(
|
||||||
[torch.cuda.memory_reserved(device), torch.cuda.memory_allocated(device)]
|
[torch.cuda.memory_reserved(device), torch.cuda.memory_allocated(device)]
|
||||||
)
|
)
|
||||||
mylogger.info(f"CUDA memory: {free_mem//1024} MByte")
|
mylogger.info(f"CUDA memory: {free_mem // 1024} MByte")
|
||||||
|
|
||||||
overwrite_fft_gauss: None | torch.Tensor = None
|
overwrite_fft_gauss: None | torch.Tensor = None
|
||||||
for i in range(0, data_filtered.shape[0]):
|
for i in range(0, data_filtered.shape[0]):
|
||||||
|
|
@ -703,7 +714,7 @@ def process_trial(
|
||||||
free_mem = cuda_total_memory - max(
|
free_mem = cuda_total_memory - max(
|
||||||
[torch.cuda.memory_reserved(device), torch.cuda.memory_allocated(device)]
|
[torch.cuda.memory_reserved(device), torch.cuda.memory_allocated(device)]
|
||||||
)
|
)
|
||||||
mylogger.info(f"CUDA memory: {free_mem//1024} MByte")
|
mylogger.info(f"CUDA memory: {free_mem // 1024} MByte")
|
||||||
mylogger.info("-==- Done -==-")
|
mylogger.info("-==- Done -==-")
|
||||||
|
|
||||||
mylogger.info("Preperation for Regression")
|
mylogger.info("Preperation for Regression")
|
||||||
|
|
@ -747,6 +758,9 @@ def process_trial(
|
||||||
mylogger.info("-==- Done -==-")
|
mylogger.info("-==- Done -==-")
|
||||||
else:
|
else:
|
||||||
dual_signal_mode = False
|
dual_signal_mode = False
|
||||||
|
target_id = config["required_order"].index("acceptor")
|
||||||
|
data_acceptor = data[target_id, ...].clone()
|
||||||
|
data_acceptor[mask_negative, :] = 0.0
|
||||||
|
|
||||||
if len(config["target_camera_donor"]) > 0:
|
if len(config["target_camera_donor"]) > 0:
|
||||||
mylogger.info("Regression Donor")
|
mylogger.info("Regression Donor")
|
||||||
|
|
@ -781,6 +795,9 @@ def process_trial(
|
||||||
mylogger.info("-==- Done -==-")
|
mylogger.info("-==- Done -==-")
|
||||||
else:
|
else:
|
||||||
dual_signal_mode = False
|
dual_signal_mode = False
|
||||||
|
target_id = config["required_order"].index("donor")
|
||||||
|
data_donor = data[target_id, ...].clone()
|
||||||
|
data_donor[mask_negative, :] = 0.0
|
||||||
|
|
||||||
del data
|
del data
|
||||||
del data_filtered
|
del data_filtered
|
||||||
|
|
@ -791,24 +808,119 @@ def process_trial(
|
||||||
free_mem = cuda_total_memory - max(
|
free_mem = cuda_total_memory - max(
|
||||||
[torch.cuda.memory_reserved(device), torch.cuda.memory_allocated(device)]
|
[torch.cuda.memory_reserved(device), torch.cuda.memory_allocated(device)]
|
||||||
)
|
)
|
||||||
mylogger.info(f"CUDA memory: {free_mem//1024} MByte")
|
mylogger.info(f"CUDA memory: {free_mem // 1024} MByte")
|
||||||
|
|
||||||
|
# #####################
|
||||||
|
|
||||||
|
if config["gevi"]:
|
||||||
|
assert dual_signal_mode
|
||||||
|
else:
|
||||||
|
assert dual_signal_mode is False
|
||||||
|
|
||||||
|
if dual_signal_mode is False:
|
||||||
|
|
||||||
|
mylogger.info("mono signal model")
|
||||||
|
|
||||||
|
mylogger.info("Remove nan")
|
||||||
|
data_acceptor = torch.nan_to_num(data_acceptor, nan=0.0)
|
||||||
|
data_donor = torch.nan_to_num(data_donor, nan=0.0)
|
||||||
|
mylogger.info("-==- Done -==-")
|
||||||
|
|
||||||
|
if config["binning_enable"] and config["binning_at_the_end"]:
|
||||||
|
mylogger.info("Binning of data")
|
||||||
|
mylogger.info(
|
||||||
|
(
|
||||||
|
f"kernel_size={int(config['binning_kernel_size'])}, "
|
||||||
|
f"stride={int(config['binning_stride'])}, "
|
||||||
|
"divisor_override=None"
|
||||||
|
)
|
||||||
|
)
|
||||||
|
|
||||||
|
data_acceptor = binning(
|
||||||
|
data_acceptor.unsqueeze(-1),
|
||||||
|
kernel_size=int(config["binning_kernel_size"]),
|
||||||
|
stride=int(config["binning_stride"]),
|
||||||
|
divisor_override=None,
|
||||||
|
).squeeze(-1)
|
||||||
|
|
||||||
|
data_donor = binning(
|
||||||
|
data_donor.unsqueeze(-1),
|
||||||
|
kernel_size=int(config["binning_kernel_size"]),
|
||||||
|
stride=int(config["binning_stride"]),
|
||||||
|
divisor_override=None,
|
||||||
|
).squeeze(-1)
|
||||||
|
|
||||||
|
mask_positve = (
|
||||||
|
binning(
|
||||||
|
mask_positve.unsqueeze(-1).unsqueeze(-1).type(dtype=dtype),
|
||||||
|
kernel_size=int(config["binning_kernel_size"]),
|
||||||
|
stride=int(config["binning_stride"]),
|
||||||
|
divisor_override=None,
|
||||||
|
)
|
||||||
|
.squeeze(-1)
|
||||||
|
.squeeze(-1)
|
||||||
|
)
|
||||||
|
mask_positve = (mask_positve > 0).type(torch.bool)
|
||||||
|
|
||||||
|
if config["save_as_python"]:
|
||||||
|
|
||||||
|
temp_path = os.path.join(
|
||||||
|
config["export_path"], experiment_name + "_acceptor_donor.npz"
|
||||||
|
)
|
||||||
|
mylogger.info(f"Save data donor and acceptor and mask to {temp_path}")
|
||||||
|
np.savez_compressed(
|
||||||
|
temp_path,
|
||||||
|
data_acceptor=data_acceptor.cpu(),
|
||||||
|
data_donor=data_donor.cpu(),
|
||||||
|
mask=mask_positve.cpu(),
|
||||||
|
)
|
||||||
|
|
||||||
|
if config["save_as_matlab"]:
|
||||||
|
temp_path = os.path.join(
|
||||||
|
config["export_path"], experiment_name + "_acceptor_donor.hd5"
|
||||||
|
)
|
||||||
|
mylogger.info(f"Save data donor and acceptor and mask to {temp_path}")
|
||||||
|
file_handle = h5py.File(temp_path, "w")
|
||||||
|
|
||||||
|
mask_positve = mask_positve.movedim(0, -1)
|
||||||
|
data_acceptor = data_acceptor.movedim(1, -1).movedim(0, -1)
|
||||||
|
data_donor = data_donor.movedim(1, -1).movedim(0, -1)
|
||||||
|
_ = file_handle.create_dataset(
|
||||||
|
"mask",
|
||||||
|
data=mask_positve.type(torch.uint8).cpu(),
|
||||||
|
compression="gzip",
|
||||||
|
compression_opts=9,
|
||||||
|
)
|
||||||
|
_ = file_handle.create_dataset(
|
||||||
|
"data_acceptor",
|
||||||
|
data=data_acceptor.cpu(),
|
||||||
|
compression="gzip",
|
||||||
|
compression_opts=9,
|
||||||
|
)
|
||||||
|
_ = file_handle.create_dataset(
|
||||||
|
"data_donor",
|
||||||
|
data=data_donor.cpu(),
|
||||||
|
compression="gzip",
|
||||||
|
compression_opts=9,
|
||||||
|
)
|
||||||
|
mylogger.info("Reminder: How to read with matlab:")
|
||||||
|
mylogger.info(f"mask = h5read('{temp_path}','/mask');")
|
||||||
|
mylogger.info(f"data_acceptor = h5read('{temp_path}','/data_acceptor');")
|
||||||
|
mylogger.info(f"data_donor = h5read('{temp_path}','/data_donor');")
|
||||||
|
file_handle.close()
|
||||||
|
return
|
||||||
|
# #####################
|
||||||
|
|
||||||
mylogger.info("Calculate ratio sequence")
|
mylogger.info("Calculate ratio sequence")
|
||||||
if dual_signal_mode:
|
|
||||||
if config["classical_ratio_mode"]:
|
if config["classical_ratio_mode"]:
|
||||||
mylogger.info("via acceptor / donor")
|
mylogger.info("via acceptor / donor")
|
||||||
ratio_sequence: torch.Tensor = data_acceptor / data_donor
|
ratio_sequence: torch.Tensor = data_acceptor / data_donor
|
||||||
mylogger.info("via / mean over time")
|
mylogger.info("via / mean over time")
|
||||||
ratio_sequence /= ratio_sequence.mean(dim=-1, keepdim=True)
|
ratio_sequence /= ratio_sequence.mean(dim=-1, keepdim=True)
|
||||||
else:
|
|
||||||
mylogger.info("via 1.0 + acceptor - donor")
|
|
||||||
ratio_sequence = 1.0 + data_acceptor - data_donor
|
|
||||||
else:
|
else:
|
||||||
mylogger.info("mono signal model")
|
mylogger.info("via 1.0 + acceptor - donor")
|
||||||
if len(config["target_camera_donor"]) > 0:
|
ratio_sequence = 1.0 + data_acceptor - data_donor
|
||||||
ratio_sequence = data_donor.clone()
|
|
||||||
else:
|
|
||||||
ratio_sequence = data_acceptor.clone()
|
|
||||||
|
|
||||||
mylogger.info("Remove nan")
|
mylogger.info("Remove nan")
|
||||||
ratio_sequence = torch.nan_to_num(ratio_sequence, nan=0.0)
|
ratio_sequence = torch.nan_to_num(ratio_sequence, nan=0.0)
|
||||||
|
|
|
||||||
Loading…
Reference in a new issue