diff --git a/DeepLabStream.py b/DeepLabStream.py index 264355b..16d3d42 100755 --- a/DeepLabStream.py +++ b/DeepLabStream.py @@ -340,7 +340,7 @@ def get_pose_mp(input_q, output_q): scmap, locref, ANIMALS_NUMBER, config ) # Use the line below to use raw DLC output rather then DLStream optimization - # peaks = pose + #peaks = pose if MODEL_ORIGIN == "MADLC": peaks = get_ma_pose(frame, config, sess, inputs, outputs) analysis_time = time.time() - start_time diff --git a/Readme.md b/Readme.md index 00896d5..10dca0e 100644 --- a/Readme.md +++ b/Readme.md @@ -13,7 +13,7 @@ DeepLabStream is a python based multi-purpose tool that enables the realtime tra Our toolbox was orginally adapted from the previously published [DeepLabCut](https://github.com/AlexEMG/DeepLabCut) ([Mathis et al., 2018](https://www.nature.com/articles/s41593-018-0209-y)) and expanded on its core capabilities, but is now able to utilize a variety of different network architectures for online pose estimation ([SLEAP](https://github.com/murthylab/sleap), [DLC-Live](https://github.com/DeepLabCut/DeepLabCut-live), [DeepPosekit's](https://github.com/jgraving/DeepPoseKit) StackedDenseNet, StackedHourGlass and [LEAP](https://github.com/murthylab/sleap)). -DeepLabStreams core feature is the utilization of real-time tracking to orchestrate closed-loop experiments. This can be achieved using any type of camera-based video stream (incl. multiple streams). It enables running experimental protocols that are dependent on a constant stream of bodypart positions and feedback activation of several input/output devices. It's capabilities range from simple region of interest (ROI) based triggers to headdirection or behavior dependent stimulation. +DeepLabStreams core feature is the utilization of real-time tracking to orchestrate closed-loop experiments. This can be achieved using any type of camera-based video stream (incl. multiple streams). It enables running experimental protocols that are dependent on a constant stream of bodypart positions and feedback activation of several input/output devices. It's capabilities range from simple region of interest (ROI) based triggers to headdirection or behavior dependent stimulation, including online classification ([SiMBA](https://www.biorxiv.org/content/10.1101/2020.04.19.049452v2), [B-SOID](https://www.biorxiv.org/content/10.1101/770271v2)). ![DLS_Stim](docs/DLSSTim_example.gif) @@ -25,6 +25,14 @@ DeepLabStreams core feature is the utilization of real-time tracking to orchestr ## New features: +#### 03/2021: Online Behavior Classification using SiMBA and B-SOID: + +- full integration of online classification of user-defined behavior using [SiMBA](https://github.com/sgoldenlab/simba) and [B-SOID](https://github.com/YttriLab/B-SOID). +- SOCIAL CLASSIFICATION with SiMBA 14bp two animal classification (more to come!) +- Unsupervised Classification with B-SOID +- New wiki guide and example experiment to get started with online classification: [Advanced Behavior Classification](https://github.com/SchwarzNeuroconLab/DeepLabStream/wiki/Advanced-Behavior-Classification) +- this version has new requirements (numba, pure, scikit-learn), so be sure to install them (e.g. `pip install -r requirements.txt`). + #### 02/2021: Multiple Animal Experiments (Pre-release): Full [SLEAP](https://github.com/murthylab/sleap) integration (Full release coming soon!) - Updated [Installation](https://github.com/SchwarzNeuroconLab/DeepLabStream/wiki/Installation-&-Testing) (for SLEAP support) @@ -33,7 +41,8 @@ DeepLabStreams core feature is the utilization of real-time tracking to orchestr #### 01/2021: DLStream was published in [Communications Biology](https://www.nature.com/articles/s42003-021-01654-9) -#### 12/2021: New pose estimation model integration ([DLC-Live](https://github.com/DeepLabCut/DeepLabCut-live)) and pre-release of further integration ([DeepPosekit's](https://github.com/jgraving/DeepPoseKit) StackedDenseNet, StackedHourGlass and [LEAP](https://github.com/murthylab/sleap)) +#### 12/2021: New pose estimation model integration + - ([DLC-Live](https://github.com/DeepLabCut/DeepLabCut-live)) and pre-release of further integration ([DeepPosekit's](https://github.com/jgraving/DeepPoseKit) StackedDenseNet, StackedHourGlass and [LEAP](https://github.com/murthylab/sleap)) ## Quick Reference: @@ -131,7 +140,6 @@ If you encounter any issues or errors, you can check out the wiki article ([Help If you use this code or data please cite: - Schweihoff, J.F., Loshakov, M., Pavlova, I. et al. DeepLabStream enables closed-loop behavioral experiments using deep learning-based markerless, real-time posture detection. Commun Biol 4, 130 (2021). https://doi.org/10.1038/s42003-021-01654-9 @@ -147,3 +155,23 @@ Developed by: - Matvey Loshakov, matveyloshakov@gmail.com Corresponding Author: Martin Schwarz, Martin.Schwarz@ukbonn.de + +## Other References + +If you are using any of the following open-source code please cite them accordingly: + +> Simple Behavioral Analysis (SimBA) – an open source toolkit for computer classification of complex social behaviors in experimental animals; + Simon RO Nilsson, Nastacia L. Goodwin, Jia Jie Choong, Sophia Hwang, Hayden R Wright, Zane C Norville, Xiaoyu Tong, Dayu Lin, Brandon S. Bentzley, Neir Eshel, Ryan J McLaughlin, Sam A. Golden + bioRxiv 2020.04.19.049452; doi: https://doi.org/10.1101/2020.04.19.049452 + +> B-SOiD: An Open Source Unsupervised Algorithm for Discovery of Spontaneous Behaviors; + Alexander I. Hsu, Eric A. Yttri + bioRxiv 770271; doi: https://doi.org/10.1101/770271 + +> SLEAP: Multi-animal pose tracking; +Talmo D. Pereira, Nathaniel Tabris, Junyu Li, Shruthi Ravindranath, Eleni S. Papadoyannis, Z. Yan Wang, David M. Turner, Grace McKenzie-Smith, Sarah D. Kocher, Annegret L. Falkner, Joshua W. Shaevitz, Mala Murthy +bioRxiv 2020.08.31.276246; doi: https://doi.org/10.1101/2020.08.31.276246 + +>Real-time, low-latency closed-loop feedback using markerless posture tracking; + Gary A Kane, Gonçalo Lopes, Jonny L Saunders, Alexander Mathis, Mackenzie W Mathis; + eLife 2020;9:e61909 doi: 10.7554/eLife.61909 diff --git a/convert_classifier.py b/convert_classifier.py new file mode 100644 index 0000000..3f8097a --- /dev/null +++ b/convert_classifier.py @@ -0,0 +1,29 @@ +import pickle +import os +from pure_sklearn.map import convert_estimator + + +def load_classifier(path_to_sav): + """Load saved classifier""" + file = open(path_to_sav, "rb") + classifier = pickle.load(file) + file.close() + return classifier + + +def convert_classifier(path): + # convert to pure python estimator + print("Loading classifier...") + clf = load_classifier(path) + dir_path = os.path.dirname(path) + filename = os.path.basename(path) + filename, _ = filename.split(".") + clf_pure_predict = convert_estimator(clf) + with open(dir_path + "/" + filename + "_pure.sav", "wb") as f: + pickle.dump(clf_pure_predict, f) + print(f"Converted Classifier {filename}") + + +if __name__ == "__main__": + path_to_classifier = "PATH_TO_CLASSIFIER" + convert_classifier(path_to_classifier) diff --git a/experiments/custom/classifier.py b/experiments/custom/classifier.py new file mode 100644 index 0000000..517c72b --- /dev/null +++ b/experiments/custom/classifier.py @@ -0,0 +1,665 @@ +""" +DeepLabStream +© J.Schweihoff, M. Loshakov +University Bonn Medical Faculty, Germany +https://github.com/SchwarzNeuroconLab/DeepLabStream +Licensed under GNU General Public License v3.0 +""" + +import multiprocessing as mp +import time +import pickle + +from utils.configloader import PATH_TO_CLASSIFIER, TIME_WINDOW +from experiments.custom.featureextraction import ( + SimbaFeatureExtractor, + SimbaFeatureExtractorStandard14bp, + BsoidFeatureExtractor, +) + + +class Classifier: + """Empty base class for classification trigger. Loads pretrained classifier, extracts features from skeleton sequence + and passes it to the classifier. Returns found motif and result if used as trigger.""" + + def __init__(self, win_len: int = 1): + self._classifier = self.load_classifier(PATH_TO_CLASSIFIER) + self._win_len = win_len + self.last_result = None + + @staticmethod + def load_classifier(path_to_sav): + """Load saved classifier""" + # import pickle + file = open(path_to_sav, "rb") + classifier = pickle.load(file) + file.close() + return classifier + + def classify(self, features): + """predicts motif from features""" + prediction = self._classifier.predict(features) + self.last_result = prediction + return prediction + + def get_last_result(self, skeleton_window: list): + """Returns predicted last prediction""" + return self.last_result + + def get_win_len(self): + return self._win_len + + +class SiMBAClassifier: + """SiMBA base class for simple behavior classification trigger. Loads pretrained classifier, gets passed features + from SimbaFeatureExtractor. Returns probability of prediction that can be incorporated into triggers.""" + + def __init__(self): + self._classifier = self.load_classifier(PATH_TO_CLASSIFIER) + self.last_result = 0.0 + self._pure = self._check_pure() + + @staticmethod + def load_classifier(path_to_sav): + """Load saved classifier""" + # import pickle + file = open(path_to_sav, "rb") + classifier = pickle.load(file) + file.close() + return classifier + + def _check_pure(self): + if "pure" in str(self._classifier): + return True + else: + return False + + def classify(self, features): + """predicts motif probability from features""" + if self._pure: + # pure-predict needs a list instead of a numpy array + prediction = self._classifier.predict_proba(list(features)) + # pure-predict returns a nested list + probability = prediction[0][1] + else: + prediction = self._classifier.predict_proba(features) + probability = prediction.item(1) + self.last_result = probability + return probability + + def get_last_result(self): + """Returns predicted last prediction""" + return self.last_result + + +class BsoidClassifier: + """BSOID base class for multiple behavior classification trigger. Loads pretrained classifier, gets passed features + from SimbaFeatureExtractor. Returns probability of prediction that can be incorporated into triggers.""" + + def __init__(self): + self._classifier = self.load_classifier(PATH_TO_CLASSIFIER) + self.last_result = 0.0 + + @staticmethod + def load_classifier(path_to_sav): + """Load saved classifier""" + import joblib + + file = open(path_to_sav, "rb") + [_, _, _, clf, _, predictions] = joblib.load(file) + file.close() + return clf + + def classify(self, features): + """predicts motif probability from features :param feats: list, multiple feats (original feature space) + :param clf: Obj, MLP classifier + :return nonfs_labels: list, label/100ms + Adapted from BSOID; https://github.com/YttriLab/B-SOID + """ + labels_fslow = [] + for i in range(0, len(features)): + labels = self._classifier.predict(features[i].T) + labels_fslow.append(labels) + self.last_result = labels_fslow + + return labels_fslow + + def get_last_result(self): + """Returns predicted last prediction""" + return self.last_result + + +"""Feature Extraction and Classification in the same pool""" + + +def example_feat_classifier_pool_run(input_q: mp.Queue, output_q: mp.Queue): + feature_extractor = SimbaFeatureExtractor(TIME_WINDOW) + classifier = Classifier() # initialize classifier + while True: + skel_time_window = None + feature_id = 0 + if input_q.full(): + skel_time_window, feature_id = input_q.get() + if skel_time_window is not None: + start_time = time.time() + features = feature_extractor.extract_features(skel_time_window) + last_prob = classifier.classify(features) + output_q.put((last_prob, feature_id)) + end_time = time.time() + # print("Classification time: {:.2f} msec".format((end_time-start_time)*1000)) + else: + pass + + +def simba_feat_classifier_pool_run(input_q: mp.Queue, output_q: mp.Queue): + feature_extractor = SimbaFeatureExtractorStandard14bp(TIME_WINDOW) + # feature_extractor = SimbaFeatureExtractor(TIME_WINDOW) + classifier = SiMBAClassifier() # initialize classifier + while True: + skel_time_window = None + feature_id = 0 + if input_q.full(): + skel_time_window, feature_id = input_q.get() + if skel_time_window is not None: + start_time = time.time() + features = feature_extractor.extract_features(skel_time_window) + # end_time = time.time() + # print( + # "Feature extraction time: {:.2f} msec".format( + # (end_time - start_time) * 1000 + # ) + # ) + + last_prob = classifier.classify(features) + output_q.put((last_prob, feature_id)) + # end_time2 = time.time() + # print( + # "Classification time: {:.2f} msec".format( + # (end_time2 - end_time) * 1000 + # ) + # ) + else: + pass + + +def bsoid_feat_classifier_pool_run(input_q: mp.Queue, output_q: mp.Queue): + feature_extractor = BsoidFeatureExtractor(TIME_WINDOW) + classifier = BsoidClassifier() # initialize classifier + while True: + skel_time_window = None + feature_id = 0 + if input_q.full(): + skel_time_window, feature_id = input_q.get() + if skel_time_window is not None: + start_time = time.time() + features = feature_extractor.extract_features(skel_time_window) + last_prob = classifier.classify(features) + output_q.put((last_prob, feature_id)) + end_time = time.time() + # print("Classification time: {:.2f} msec".format((end_time-start_time)*1000)) + # print("Feature ID: "+ feature_id) + else: + pass + + +class FeatureExtractionClassifierProcessPool: + """ + Class to help work with protocol function in multiprocessing + spawns a pool of processes that tackle the frame-by-frame issue. + """ + + def __init__(self, pool_size: int): + """ + Setting up the three queues and the process itself + """ + self._running = False + self._pool_size = pool_size + self._process_pool = self.initiate_pool( + example_feat_classifier_pool_run, pool_size + ) + + @staticmethod + def initiate_pool(process_func, pool_size: int): + """creates list of process dictionaries that are used to classify features + :param process_func: function that will be passed to mp.Process object, should contain classification + :param pool_size: number of processes created by function, should be enough to enable constistent feature classification without skipped frames + :""" + process_pool = [] + + for i in range(pool_size): + input_queue = mp.Queue(1) + output_queue = mp.Queue(1) + classification_process = mp.Process( + target=process_func, args=(input_queue, output_queue) + ) + process_pool.append( + dict( + process=classification_process, + input=input_queue, + output=output_queue, + running=False, + ) + ) + + return process_pool + + def start(self): + """ + Starting all processes + """ + for process in self._process_pool: + process["process"].start() + + def end(self): + """ + Ending all processes + """ + for process in self._process_pool: + process["input"].close() + process["output"].close() + process["process"].terminate() + + def get_status(self): + """ + Getting current status of the running protocol + """ + return self._running + + def pass_time_window(self, skel_time_window: tuple, debug: bool = False): + """ + Passing the features to the process pool + First checks if processes got their first input yet + Checks which process is already done and then gives new input + breaks for loop if an idle process was found + :param features tuple: feature list from feature extractor and feature_id used to identify processing sequence + :param debug bool: reporting of process + feature id to identify discrepancies in processing sequence + """ + for process in self._process_pool: + if not process["running"]: + if process["input"].empty(): + process["input"].put(skel_time_window) + process["running"] = True + if debug: + print( + "First Input", + process["process"].name, + "ID: " + str(skel_time_window[1]), + ) + break + + elif process["input"].empty() and process["output"].full(): + process["input"].put(skel_time_window) + if debug: + print( + "Input", + process["process"].name, + "ID: " + str(skel_time_window[1]), + ) + break + + def get_result(self, debug: bool = False): + """ + Getting result from the process pool + takes result from first finished process in pool + :param debug bool: reporting of process + feature id to identify discrepancies in processing sequence + + """ + result = (None, 0) + for process in self._process_pool: + if process["output"].full(): + result = process["output"].get() + if debug: + print("Output", process["process"].name, "ID: " + str(result[1])) + break + return result + + +class FeatSimbaProcessPool(FeatureExtractionClassifierProcessPool): + """ + Class to help work with protocol function in multiprocessing + spawns a pool of processes that tackle the frame-by-frame issue. + """ + + def __init__(self, pool_size: int): + """ + Setting up the three queues and the process itself + """ + super().__init__(pool_size) + self._process_pool = super().initiate_pool( + simba_feat_classifier_pool_run, pool_size + ) + + +class FeatBsoidProcessPool(FeatureExtractionClassifierProcessPool): + """ + Class to help work with protocol function in multiprocessing + spawns a pool of processes that tackle the frame-by-frame issue. + """ + + def __init__(self, pool_size: int): + """ + Setting up the three queues and the process itself + """ + super().__init__(pool_size) + self._process_pool = super().initiate_pool( + bsoid_feat_classifier_pool_run, pool_size + ) + + +"""Simple process protocols """ + + +def example_classifier_run( + input_classification_q: mp.Queue, output_classification_q: mp.Queue +): + classifier = Classifier() # initialize classifier + while True: + features = None + if input_classification_q.full(): + features = input_classification_q.get() + if features is not None: + last_prob = classifier.classify(features) + output_classification_q.put(last_prob) + else: + pass + + +def simba_classifier_run(input_q: mp.Queue, output_q: mp.Queue): + classifier = SiMBAClassifier() # initialize classifier + while True: + features = None + if input_q.full(): + features = input_q.get() + if features is not None: + start_time = time.time() + last_prob = classifier.classify(features) + output_q.put((last_prob)) + end_time = time.time() + # print("Classification time: {:.2f} msec".format((end_time-start_time)*1000)) + else: + pass + + +def bsoid_classifier_run(input_q: mp.Queue, output_q: mp.Queue): + classifier = BsoidClassifier() # initialize classifier + while True: + features = None + if input_q.full(): + features = input_q.get() + if features is not None: + start_time = time.time() + last_prob = classifier.classify(features) + output_q.put((last_prob)) + end_time = time.time() + print( + "Classification time: {:.2f} msec".format( + (end_time - start_time) * 1000 + ) + ) + else: + pass + + +class ClassifierProcess: + """ + Class to help work with protocol function in multiprocessing + Modified from stimulus_process.py + """ + + def __init__(self): + """ + Setting up the three queues and the process itself + """ + self.input_queue = mp.Queue(1) + self.output_queue = mp.Queue(1) + self._classification_process = None + self._running = False + self._classification_process = mp.Process( + target=example_classifier_run, args=(self.input_queue, self.output_queue) + ) + + def start(self): + """ + Starting the process + """ + self._classification_process.start() + + def end(self): + """ + Ending the process + """ + self.input_queue.close() + self.output_queue.close() + self._classification_process.terminate() + + def get_status(self): + """ + Getting current status of the running protocol + """ + return self._running + + def pass_features(self, features): + """ + Passing the features to the process + """ + if self.input_queue.empty(): + self.input_queue.put(features) + self._running = True + + def get_result(self): + """ + Getting result from the process + """ + if self.output_queue.full(): + self._running = False + return self.output_queue.get() + + +class SimbaClassifier_Process(ClassifierProcess): + def __init__(self): + super().__init__() + self.input_queue = mp.Queue(1) + self.output_queue = mp.Queue(1) + self._classification_process = mp.Process( + target=simba_classifier_run, args=(self.input_queue, self.output_queue) + ) + + +class BsoidClassifier_Process(ClassifierProcess): + def __init__(self): + super().__init__() + self.input_queue = mp.Queue(1) + self.output_queue = mp.Queue(1) + self._classification_process = mp.Process( + target=bsoid_classifier_run, args=(self.input_queue, self.output_queue) + ) + + +"""Processing pool for classification""" + + +def example_classifier_pool_run(input_q: mp.Queue, output_q: mp.Queue): + classifier = Classifier() # initialize classifier + while True: + features = None + feature_id = 0 + if input_q.full(): + features, feature_id = input_q.get() + if features is not None: + start_time = time.time() + last_prob = classifier.classify(features) + output_q.put((last_prob, feature_id)) + end_time = time.time() + # print("Classification time: {:.2f} msec".format((end_time-start_time)*1000)) + else: + pass + + +def simba_classifier_pool_run(input_q: mp.Queue, output_q: mp.Queue): + classifier = SiMBAClassifier() # initialize classifier + while True: + features = None + feature_id = 0 + if input_q.full(): + features, feature_id = input_q.get() + if features is not None: + start_time = time.time() + last_prob = classifier.classify(features) + output_q.put((last_prob, feature_id)) + end_time = time.time() + # print("Classification time: {:.2f} msec".format((end_time-start_time)*1000)) + else: + pass + + +def bsoid_classifier_pool_run(input_q: mp.Queue, output_q: mp.Queue): + classifier = BsoidClassifier() # initialize classifier + while True: + features = None + feature_id = 0 + if input_q.full(): + features, feature_id = input_q.get() + if features is not None: + start_time = time.time() + last_prob = classifier.classify(features) + output_q.put((last_prob, feature_id)) + end_time = time.time() + # print("Classification time: {:.2f} msec".format((end_time-start_time)*1000)) + # print("Feature ID: "+ feature_id) + + else: + pass + + +class ClassifierProcessPool: + """ + Class to help work with protocol function in multiprocessing + spawns a pool of processes that tackle the frame-by-frame issue. + """ + + def __init__(self, pool_size: int): + """ + Setting up the three queues and the process itself + """ + self._running = False + self._pool_size = pool_size + self._process_pool = self.initiate_pool(example_classifier_pool_run, pool_size) + + @staticmethod + def initiate_pool(process_func, pool_size: int): + """creates list of process dictionaries that are used to classify features + :param process_func: function that will be passed to mp.Process object, should contain classification + :param pool_size: number of processes created by function, should be enough to enable constistent feature classification without skipped frames + :""" + process_pool = [] + + for i in range(pool_size): + input_queue = mp.Queue(1) + output_queue = mp.Queue(1) + classification_process = mp.Process( + target=process_func, args=(input_queue, output_queue) + ) + process_pool.append( + dict( + process=classification_process, + input=input_queue, + output=output_queue, + running=False, + ) + ) + + return process_pool + + def start(self): + """ + Starting all processes + """ + for process in self._process_pool: + process["process"].start() + + def end(self): + """ + Ending all processes + """ + for process in self._process_pool: + process["input"].close() + process["output"].close() + process["process"].terminate() + + def get_status(self): + """ + Getting current status of the running protocol + """ + return self._running + + def pass_features(self, features: tuple, debug: bool = False): + """ + Passing the features to the process pool + First checks if processes got their first input yet + Checks which process is already done and then gives new input + breaks for loop if an idle process was found + :param features tuple: feature list from feature extractor and feature_id used to identify processing sequence + :param debug bool: reporting of process + feature id to identify discrepancies in processing sequence + """ + for process in self._process_pool: + if not process["running"]: + if process["input"].empty(): + process["input"].put(features) + process["running"] = True + if debug: + print( + "First Input", + process["process"].name, + "ID: " + str(features[1]), + ) + break + + elif process["input"].empty() and process["output"].full(): + process["input"].put(features) + if debug: + print("Input", process["process"].name, "ID: " + str(features[1])) + break + + def get_result(self, debug: bool = False): + """ + Getting result from the process pool + takes result from first finished process in pool + :param debug bool: reporting of process + feature id to identify discrepancies in processing sequence + + """ + result = (None, 0) + for process in self._process_pool: + if process["output"].full(): + result = process["output"].get() + if debug: + print("Output", process["process"].name, "ID: " + str(result[1])) + break + return result + + +class SimbaProcessPool(ClassifierProcessPool): + """ + Class to help work with protocol function in multiprocessing + spawns a pool of processes that tackle the frame-by-frame issue. + """ + + def __init__(self, pool_size: int): + """ + Setting up the three queues and the process itself + """ + super().__init__(pool_size) + self._process_pool = super().initiate_pool(simba_classifier_pool_run, pool_size) + + +class BsoidProcessPool(ClassifierProcessPool): + """ + Class to help work with protocol function in multiprocessing + spawns a pool of processes that tackle the frame-by-frame issue. + """ + + def __init__(self, pool_size: int): + """ + Setting up the three queues and the process itself + """ + super().__init__(pool_size) + self._process_pool = super().initiate_pool(bsoid_classifier_pool_run, pool_size) diff --git a/experiments/custom/experiments.py b/experiments/custom/experiments.py index 34d3050..fcc7333 100644 --- a/experiments/custom/experiments.py +++ b/experiments/custom/experiments.py @@ -22,11 +22,262 @@ OutsideTrigger, DirectionTrigger, SpeedTrigger, + SimbaThresholdBehaviorPoolTrigger, + BsoidClassBehaviorPoolTrigger, SocialInteractionTrigger, ) + from utils.plotter import plot_triggers_response from utils.analysis import angle_between_vectors from experiments.custom.stimulation import show_visual_stim_img, laser_switch +from experiments.custom.classifier import FeatBsoidProcessPool, FeatSimbaProcessPool + + +from utils.configloader import THRESHOLD, POOL_SIZE, TRIGGER + + +""" experimental classification experiment using Simba trained classifiers in a pool which are converted using the pure-predict package""" + + +class SimbaBehaviorPoolExperiment: + """ + Test experiment for Simba classification + Simple class to contain all of the experiment properties and includes classification + Uses multiprocess to ensure the best possible performance and + to showcase that it is possible to work with any type of equipment, even timer-dependant + """ + + def __init__(self): + """Classifier process and initiation of behavior trigger""" + self.experiment_finished = False + self._process_experiment = ExampleProtocolProcess() + # this process has feature extraction and classification in one process + # simplifies everything if the feature extraction script is within the parallel process + self._process_pool = FeatSimbaProcessPool(POOL_SIZE) + # pass classifier to trigger, so that check_skeleton is the only function that passes skeleton + # initiate in experiment, so that process can be started with start_experiment + self._behaviortrigger = SimbaThresholdBehaviorPoolTrigger( + prob_threshold=THRESHOLD, class_process_pool=self._process_pool, debug=False + ) + self._event = None + # is not fully utilized in this experiment but is useful to keep for further adaptation + self._current_trial = None + self._max_reps = 999 + self._trial_count = {trial: 0 for trial in self._trials} + self._trial_timers = {trial: Timer(0) for trial in self._trials} + self._exp_timer = Timer(9999) + + def check_skeleton(self, frame, skeleton): + """ + Checking each passed animal skeleton for a pre-defined set of conditions + Outputting the visual representation, if exist + Advancing trials according to inherent logic of an experiment + :param frame: frame, on which animal skeleton was found + :param skeleton: skeleton, consisting of multiple joints of an animal + """ + self.check_exp_timer() # checking if experiment is still on + for trial in self._trial_count: + # checking if any trial hit a predefined cap + if self._trial_count[trial] >= self._max_reps: + self.stop_experiment() + + if not self.experiment_finished: + for trial in self._trials: + # check for all trials if condition is met + # this passes the skeleton to the trigger, where the feature extraction is done and the extracted features + # are passed to the classifier process + result, response = self._trials[trial]["trigger"]( + skeleton, target_prob=self._trials[trial]["target_prob"] + ) + plot_triggers_response(frame, response) + # if the trigger is reporting back that the behavior is found: do something + # currently nothing is done, expect counting the occurances + if result: + if self._current_trial is None: + if not self._trial_timers[trial].check_timer(): + self._current_trial = trial + self._trial_timers[trial].reset() + self._trial_count[trial] += 1 + print(trial, self._trial_count[trial]) + else: + if self._current_trial == trial: + self._current_trial = None + self._trial_timers[trial].start() + + self._process_experiment.set_trial(self._current_trial) + + @property + def _trials(self): + """ + Defining the trials + """ + trials = { + "DLStream_test": dict( + trigger=self._behaviortrigger.check_skeleton, target_prob=None, count=0 + ) + } + return trials + + def check_exp_timer(self): + """ + Checking the experiment timer + """ + if not self._exp_timer.check_timer(): + print("Experiment is finished") + print("Time ran out.") + self.stop_experiment() + + def start_experiment(self): + """ + Start the experiment + """ + self._process_experiment.start() + self._process_pool.start() + if not self.experiment_finished: + self._exp_timer.start() + + def stop_experiment(self): + """ + Stop the experiment and reset the timer + """ + self.experiment_finished = True + self._process_experiment.end() + self._process_pool.end() + print("Experiment completed!") + self._exp_timer.reset() + + def get_trial(self): + """ + Check which trial is going on right now + """ + return self._event + + def get_info(self): + """ returns optional info""" + info = self._behaviortrigger.get_last_prob() + return info + + +""" experimental classification experiment using BSOID trained classifiers in a pool""" + + +class BsoidBehaviorPoolExperiment: + """ + Test experiment for BSOID classification + Simple class to contain all of the experiment properties and includes classification + Uses multiprocess to ensure the best possible performance and + to showcase that it is possible to work with any type of equipment, even timer-dependant + """ + + def __init__(self): + """Classifier process and initiation of behavior trigger""" + self.experiment_finished = False + self._process_experiment = ExampleProtocolProcess() + self._process_pool = FeatBsoidProcessPool(POOL_SIZE) + # pass classifier to trigger, so that check_skeleton is the only function that passes skeleton + # initiate in experiment, so that process can be started with start_experiment + self._behaviortrigger = BsoidClassBehaviorPoolTrigger( + target_class=TRIGGER, class_process_pool=self._process_pool, debug=False + ) + self._event = None + # is not fully utilized in this experiment but is usefull to keep for further adaptation + self._current_trial = None + self._trial_count = {trial: 0 for trial in self._trials} + self._trial_timers = {trial: Timer(10) for trial in self._trials} + self._exp_timer = Timer(600) + + def check_skeleton(self, frame, skeleton): + """ + Checking each passed animal skeleton for a pre-defined set of conditions + Outputting the visual representation, if exist + Advancing trials according to inherent logic of an experiment + :param frame: frame, on which animal skeleton was found + :param skeleton: skeleton, consisting of multiple joints of an animal + """ + self.check_exp_timer() # checking if experiment is still on + for trial in self._trial_count: + # checking if any trial hit a predefined cap + if self._trial_count[trial] >= 10: + self.stop_experiment() + + if not self.experiment_finished: + for trial in self._trials: + # check for all trials if condition is met + # this passes the skeleton to the trigger, where the feature extraction is done and the extracted features + # are passed to the classifier process + result, response = self._trials[trial]["trigger"]( + skeleton, target_class=self._trials[trial]["target_class"] + ) + plot_triggers_response(frame, response) + # if the trigger is reporting back that the behavior is found: do something + # currently nothing is done, expect counting the occurances + if result: + if self._current_trial is None: + if not self._trial_timers[trial].check_timer(): + self._current_trial = trial + self._trial_timers[trial].reset() + self._trial_count[trial] += 1 + print(trial, self._trial_count[trial]) + else: + if self._current_trial == trial: + self._current_trial = None + self._trial_timers[trial].start() + self._process_experiment.set_trial(self._current_trial) + else: + pass + return result, response + + @property + def _trials(self): + """ + Defining the trials + Target class is the cluster of interest and can be changed with every call of check_skeleton + """ + trials = { + "DLStream_test": dict( + trigger=self._behaviortrigger.check_skeleton, target_class=None, count=0 + ) + } + return trials + + def check_exp_timer(self): + """ + Checking the experiment timer + """ + if not self._exp_timer.check_timer(): + print("Experiment is finished") + print("Time ran out.") + self.stop_experiment() + + def start_experiment(self): + """ + Start the experiment + """ + self._process_pool.start() + self._process_experiment.start() + if not self.experiment_finished: + self._exp_timer.start() + + def stop_experiment(self): + """ + Stop the experiment and reset the timer + """ + self.experiment_finished = True + self._process_experiment.end() + self._process_pool.end() + print("Experiment completed!") + self._exp_timer.reset() + + def get_trial(self): + """ + Check which trial is going on right now + """ + return self._event + + def get_info(self): + """ returns optional info""" + info = self._behaviortrigger.get_last_prob() + return info """Social or multiple animal experiments in combination with SLEAP or non-flattened maDLC pose estimation""" diff --git a/experiments/custom/featureextraction.py b/experiments/custom/featureextraction.py new file mode 100644 index 0000000..ecedd1e --- /dev/null +++ b/experiments/custom/featureextraction.py @@ -0,0 +1,1435 @@ +""" +DeepLabStream +© J.Schweihoff, M. Loshakov +University Bonn Medical Faculty, Germany +https://github.com/SchwarzNeuroconLab/DeepLabStream +Licensed under GNU General Public License v3.0 + +Simba Feature extraction functions were provided by Simon Nilsson from Golden Lab +Main developer of SiMBA https://github.com/sgoldenlab/simba +and integrated into the SimbaFeatureExtractor + +Bsoid Feature extraction functions were adpated from code orginally by Alexander Hsu from Yttri Lab +Main developer of B-Soid https://github.com/YttriLab/B-SOID +and integrated into the BsoidFeatureExtractor +""" +import scipy +from scipy.spatial import ConvexHull +from utils.configloader import PIXPERMM, FRAMERATE +import numpy as np +import pandas as pd +import math +import time +from numba import jit +import itertools + +"""For standard SIMBA feature extraction""" + + +def count_values_in_range(series, values_in_range_min, values_in_range_max): + return series.between(left=values_in_range_min, right=values_in_range_max).sum() + + +@jit(nopython=True, cache=True) +def euclidean_distance(x1, x2, y1, y2): + result = np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2) + return result + + +@jit(nopython=True, cache=True) +def angle3pt(ax, ay, bx, by, cx, cy): + ang = math.degrees(math.atan2(cy - by, cx - bx) - math.atan2(ay - by, ax - bx)) + return ang + 360 if ang < 0 else ang + + +def as_strided(*args): + return np.lib.stride_tricks.as_strided(*args) + + +class SimbaFeatureExtractorStandard14bp: + """Feature extraction module for integration in behavior trigger, takes list of postures as input + and calculates features to pass to classifier. Features and classifier have to match! + Designed to work with Simba https://github.com/sgoldenlab/simba""" + + def __init__(self, input_array_length, debug=False): + self._debug = debug + self.input_array_length = input_array_length + # TODO: Collect bodypart position in input array from skeleton automatically (this will make this much easier!) + self._currPixPerMM = PIXPERMM + self._fps = FRAMERATE + self._num_bodyparts = 14 + self._roll_window_values = [2, 5, 6, 7.5, 15] + self._roll_window_values_len = len(self._roll_window_values) + self._roll_windows = [int(self._fps / num) for num in self._roll_window_values] + self._columnheaders = [ + "Nose_1_x", + "Nose_1_y", + "Ear_left_1_x", + "Ear_left_1_y", + "Ear_right_1_x", + "Ear_right_1_y", + "Center_1_x", + "Center_1_y", + "Lat_left_1_x", + "Lat_left_1_y", + "Lat_right_1_x", + "Lat_right_1_y", + "Tail_base_1_x", + "Tail_base_1_y", + "Nose_2_x", + "Nose_2_y", + "Ear_left_2_x", + "Ear_left_2_y", + "Ear_right_2_x", + "Ear_right_2_y", + "Center_2_x", + "Center_2_y", + "Lat_left_2_x", + "Lat_left_2_y", + "Lat_right_2_x", + "Lat_right_2_y", + "Tail_base_2_x", + "Tail_base_2_y", + ] + self._pheaders = [ + "Nose_1_p", + "Ear_left_1_p", + "Ear_right_1_p", + "Center_1_p", + "Lat_left_1_p", + "Lat_right_1_p", + "Tail_base_1_p", + "Nose_2_p", + "Ear_left_2_p", + "Ear_right_2_p", + "Center_2_p", + "Lat_left_2_p", + "Lat_right_2_p", + "Tail_base_2_p", + ] + + def convert_pandas(self, input_array): + """This is a hardcoded version for the 14 bp SiMBA classification""" + input_array_np = np.array(input_array) + + input_df = pd.DataFrame(input_array_np, columns=self._columnheaders) + # add fake prediction column + for p_clm in self._pheaders: + input_df[p_clm] = 0.99 + + input_df = input_df.fillna(0) + input_df = input_df.drop(input_df.index[[0]]) + input_df = input_df.apply(pd.to_numeric) + input_df = input_df.reset_index() + input_df = input_df.reset_index(drop=True) + return input_df + + def extract_features_simba14bp(self, input_df): + # adapted from SIMBA: https://github.com/sgoldenlab/simba/blob/master/simba/features_scripts/extract_features_14bp.py + ########### CREATE PD FOR RAW DATA AND PD FOR MOVEMENT BETWEEN FRAMES ########### + M1_hull_large_euclidean_list = [] + M1_hull_small_euclidean_list = [] + M1_hull_mean_euclidean_list = [] + M1_hull_sum_euclidean_list = [] + M2_hull_large_euclidean_list = [] + M2_hull_small_euclidean_list = [] + M2_hull_mean_euclidean_list = [] + M2_hull_sum_euclidean_list = [] + + start_time = time.time() + ########### MOUSE AREAS ########################################### + input_df["Mouse_1_poly_area"] = input_df.apply( + lambda x: ConvexHull( + np.array( + [ + [x["Ear_left_1_x"], x["Ear_left_1_y"]], + [x["Ear_right_1_x"], x["Ear_right_1_y"]], + [x["Nose_1_x"], x["Nose_1_y"]], + [x["Lat_left_1_x"], x["Lat_left_1_y"]], + [x["Lat_right_1_x"], x["Lat_right_1_y"]], + [x["Tail_base_1_x"], x["Tail_base_1_y"]], + [x["Center_1_x"], x["Center_1_y"]], + ] + ) + ).area, + axis=1, + ) + input_df["Mouse_1_poly_area"] = ( + input_df["Mouse_1_poly_area"] / self._currPixPerMM + ) + + input_df["Mouse_2_poly_area"] = input_df.apply( + lambda x: ConvexHull( + np.array( + [ + [x["Ear_left_2_x"], x["Ear_left_2_y"]], + [x["Ear_right_2_x"], x["Ear_right_2_y"]], + [x["Nose_2_x"], x["Nose_2_y"]], + [x["Lat_left_2_x"], x["Lat_left_2_y"]], + [x["Lat_right_2_x"], x["Lat_right_2_y"]], + [x["Tail_base_2_x"], x["Tail_base_2_y"]], + [x["Center_2_x"], x["Center_2_y"]], + ] + ) + ).area, + axis=1, + ) + input_df["Mouse_2_poly_area"] = ( + input_df["Mouse_2_poly_area"] / self._currPixPerMM + ) + + if self._debug: + print("Evaluating convex hulls...") + print((time.time() - start_time) * 1000, "ms") + + ########### CREATE SHIFTED DATAFRAME FOR DISTANCE CALCULATIONS ########################################### + start_time = time.time() + input_df_shifted = input_df.shift(periods=1) + input_df_shifted = input_df_shifted.rename( + columns={ + "Ear_left_1_x": "Ear_left_1_x_shifted", + "Ear_left_1_y": "Ear_left_1_y_shifted", + "Ear_left_1_p": "Ear_left_1_p_shifted", + "Ear_right_1_x": "Ear_right_1_x_shifted", + "Ear_right_1_y": "Ear_right_1_y_shifted", + "Ear_right_1_p": "Ear_right_1_p_shifted", + "Nose_1_x": "Nose_1_x_shifted", + "Nose_1_y": "Nose_1_y_shifted", + "Nose_1_p": "Nose_1_p_shifted", + "Center_1_x": "Center_1_x_shifted", + "Center_1_y": "Center_1_y_shifted", + "Center_1_p": "Center_1_p_shifted", + "Lat_left_1_x": "Lat_left_1_x_shifted", + "Lat_left_1_y": "Lat_left_1_y_shifted", + "Lat_left_1_p": "Lat_left_1_p_shifted", + "Lat_right_1_x": "Lat_right_1_x_shifted", + "Lat_right_1_y": "Lat_right_1_y_shifted", + "Lat_right_1_p": "Lat_right_1_p_shifted", + "Tail_base_1_x": "Tail_base_1_x_shifted", + "Tail_base_1_y": "Tail_base_1_y_shifted", + "Tail_base_1_p": "Tail_base_1_p_shifted", + "Ear_left_2_x": "Ear_left_2_x_shifted", + "Ear_left_2_y": "Ear_left_2_y_shifted", + "Ear_left_2_p": "Ear_left_2_p_shifted", + "Ear_right_2_x": "Ear_right_2_x_shifted", + "Ear_right_2_y": "Ear_right_2_y_shifted", + "Ear_right_2_p": "Ear_right_2_p_shifted", + "Nose_2_x": "Nose_2_x_shifted", + "Nose_2_y": "Nose_2_y_shifted", + "Nose_2_p": "Nose_2_p_shifted", + "Center_2_x": "Center_2_x_shifted", + "Center_2_y": "Center_2_y_shifted", + "Center_2_p": "Center_2_p_shifted", + "Lat_left_2_x": "Lat_left_2_x_shifted", + "Lat_left_2_y": "Lat_left_2_y_shifted", + "Lat_left_2_p": "Lat_left_2_p_shifted", + "Lat_right_2_x": "Lat_right_2_x_shifted", + "Lat_right_2_y": "Lat_right_2_y_shifted", + "Lat_right_2_p": "Lat_right_2_p_shifted", + "Tail_base_2_x": "Tail_base_2_x_shifted", + "Tail_base_2_y": "Tail_base_2_y_shifted", + "Tail_base_2_p": "Tail_base_2_p_shifted", + "Mouse_1_poly_area": "Mouse_1_poly_area_shifted", + "Mouse_2_poly_area": "Mouse_2_poly_area_shifted", + } + ) + input_df_combined = pd.concat( + [input_df, input_df_shifted], axis=1, join="inner" + ) + input_df_combined = input_df_combined.fillna(0) + input_df_combined = input_df_combined.reset_index(drop=True) + + if self._debug: + print("Creating shifted dataframes for distance calculations") + print((time.time() - start_time) * 1000, "ms") + + ########### EUCLIDEAN DISTANCES ########################################### + start_time = time.time() + + # within mice + + eucl_distance_dict_wm = dict( + nose_to_tail=("Nose", "Tail_base"), + width=("Lat_left", "Lat_right"), + Ear_distance=("Ear_right", "Ear_left"), + Nose_to_centroid=("Nose", "Center"), + Nose_to_lateral_left=("Nose", "Lat_left"), + Nose_to_lateral_right=("Nose", "Lat_right"), + Centroid_to_lateral_left=("Center", "Lat_left"), + Centroid_to_lateral_right=("Center", "Lat_right"), + ) + + mice = [1, 2] + for mouse in mice: + for distance_measurement, bodyparts in eucl_distance_dict_wm.items(): + x1 = input_df[f"{bodyparts[0]}_{mouse}_x"].to_numpy() + y1 = input_df[f"{bodyparts[0]}_{mouse}_y"].to_numpy() + x2 = input_df[f"{bodyparts[1]}_{mouse}_x"].to_numpy() + y2 = input_df[f"{bodyparts[1]}_{mouse}_y"].to_numpy() + input_df[f"Mouse_{mouse}_{distance_measurement}"] = ( + euclidean_distance(x1, x2, y1, y2) / self._currPixPerMM + ) + + # between mice + + eucl_distance_dict_bm = dict( + Centroid_distance=("Center_1", "Center_2"), + Nose_to_nose_distance=("Nose_1", "Nose_2"), + M1_Nose_to_M2_lat_left=("Nose_1", "Lat_left_2"), + M1_Nose_to_M2_lat_right=("Nose_1", "Lat_right_2"), + M2_Nose_to_M1_lat_left=("Nose_2", "Lat_left_1"), + M2_Nose_to_M1_lat_right=("Nose_2", "Lat_right_1"), + M1_Nose_to_M2_tail_base=("Nose_1", "Tail_base_2"), + M2_Nose_to_M1_tail_base=("Nose_2", "Tail_base_1"), + ) + + for distance_measurement, bodyparts in eucl_distance_dict_bm.items(): + x1 = input_df[f"{bodyparts[0]}_x"].to_numpy() + y1 = input_df[f"{bodyparts[0]}_y"].to_numpy() + x2 = input_df[f"{bodyparts[1]}_x"].to_numpy() + y2 = input_df[f"{bodyparts[1]}_y"].to_numpy() + input_df[f"{distance_measurement}"] = ( + euclidean_distance(x1, x2, y1, y2) / self._currPixPerMM + ) + + # Movement + + bp_list = ( + "Center", + "Nose", + "Lat_left", + "Lat_right", + "Tail_base", + "Ear_left", + "Ear_right", + ) + + mice = [1, 2] + for mouse in mice: + for bp in bp_list: + x1 = input_df_combined[f"{bp}_{mouse}_x_shifted"].to_numpy() + y1 = input_df_combined[f"{bp}_{mouse}_y_shifted"].to_numpy() + x2 = input_df_combined[f"{bp}_{mouse}_x"].to_numpy() + y2 = input_df_combined[f"{bp}_{mouse}_y"].to_numpy() + "Movement_mouse_1_centroid" + if bp == "Center": + input_df[f"Movement_mouse_{mouse}_centroid"] = ( + euclidean_distance(x1, x2, y1, y2) / self._currPixPerMM + ) + elif bp == "Ear_left": + input_df[f"Movement_mouse_{mouse}_left_ear"] = ( + euclidean_distance(x1, x2, y1, y2) / self._currPixPerMM + ) + elif bp == "Ear_right": + input_df[f"Movement_mouse_{mouse}_right_ear"] = ( + euclidean_distance(x1, x2, y1, y2) / self._currPixPerMM + ) + elif bp == "Lat_left": + input_df[f"Movement_mouse_{mouse}_lateral_left"] = ( + euclidean_distance(x1, x2, y1, y2) / self._currPixPerMM + ) + elif bp == "Lat_right": + input_df[f"Movement_mouse_{mouse}_lateral_right"] = ( + euclidean_distance(x1, x2, y1, y2) / self._currPixPerMM + ) + else: + input_df[f"Movement_mouse_{mouse}_{bp.lower()}"] = ( + euclidean_distance(x1, x2, y1, y2) / self._currPixPerMM + ) + + input_df["Mouse_1_polygon_size_change"] = ( + input_df_combined["Mouse_1_poly_area_shifted"] + - input_df_combined["Mouse_1_poly_area"] + ) + input_df["Mouse_2_polygon_size_change"] = ( + input_df_combined["Mouse_2_poly_area_shifted"] + - input_df_combined["Mouse_2_poly_area"] + ) + + if self._debug: + print("Calculating euclidean distances...") + print((time.time() - start_time) * 1000, "ms") + + ########### HULL - EUCLIDEAN DISTANCES ########################################### + start_time = time.time() + + for index, row in input_df.iterrows(): + M1_np_array = np.array( + [ + [row["Ear_left_1_x"], row["Ear_left_1_y"]], + [row["Ear_right_1_x"], row["Ear_right_1_y"]], + [row["Nose_1_x"], row["Nose_1_y"]], + [row["Center_1_x"], row["Center_1_y"]], + [row["Lat_left_1_x"], row["Lat_left_1_y"]], + [row["Lat_right_1_x"], row["Lat_right_1_y"]], + [row["Tail_base_1_x"], row["Tail_base_1_y"]], + ] + ).astype(int) + M2_np_array = np.array( + [ + [row["Ear_left_2_x"], row["Ear_left_2_y"]], + [row["Ear_right_2_x"], row["Ear_right_2_y"]], + [row["Nose_2_x"], row["Nose_2_y"]], + [row["Center_2_x"], row["Center_2_y"]], + [row["Lat_left_2_x"], row["Lat_left_2_y"]], + [row["Lat_right_2_x"], row["Lat_right_2_y"]], + [row["Tail_base_2_x"], row["Tail_base_2_y"]], + ] + ).astype(int) + M1_dist_euclidean = scipy.spatial.distance.cdist( + M1_np_array, M1_np_array, metric="euclidean" + ) + M1_dist_euclidean = M1_dist_euclidean[M1_dist_euclidean != 0] + M1_hull_large_euclidean = np.amax(M1_dist_euclidean) + M1_hull_small_euclidean = np.min(M1_dist_euclidean) + M1_hull_mean_euclidean = np.mean(M1_dist_euclidean) + M1_hull_sum_euclidean = np.sum(M1_dist_euclidean) + M1_hull_large_euclidean_list.append(M1_hull_large_euclidean) + M1_hull_small_euclidean_list.append(M1_hull_small_euclidean) + M1_hull_mean_euclidean_list.append(M1_hull_mean_euclidean) + M1_hull_sum_euclidean_list.append(M1_hull_sum_euclidean) + M2_dist_euclidean = scipy.spatial.distance.cdist( + M2_np_array, M2_np_array, metric="euclidean" + ) + M2_dist_euclidean = M2_dist_euclidean[M2_dist_euclidean != 0] + M2_hull_large_euclidean = np.amax(M2_dist_euclidean) + M2_hull_small_euclidean = np.min(M2_dist_euclidean) + M2_hull_mean_euclidean = np.mean(M2_dist_euclidean) + M2_hull_sum_euclidean = np.sum(M2_dist_euclidean) + M2_hull_large_euclidean_list.append(M2_hull_large_euclidean) + M2_hull_small_euclidean_list.append(M2_hull_small_euclidean) + M2_hull_mean_euclidean_list.append(M2_hull_mean_euclidean) + M2_hull_sum_euclidean_list.append(M2_hull_sum_euclidean) + input_df["M1_largest_euclidean_distance_hull"] = list( + map(lambda x: x / self._currPixPerMM, M1_hull_large_euclidean_list) + ) + input_df["M1_smallest_euclidean_distance_hull"] = list( + map(lambda x: x / self._currPixPerMM, M1_hull_small_euclidean_list) + ) + input_df["M1_mean_euclidean_distance_hull"] = list( + map(lambda x: x / self._currPixPerMM, M1_hull_mean_euclidean_list) + ) + input_df["M1_sum_euclidean_distance_hull"] = list( + map(lambda x: x / self._currPixPerMM, M1_hull_sum_euclidean_list) + ) + input_df["M2_largest_euclidean_distance_hull"] = list( + map(lambda x: x / self._currPixPerMM, M2_hull_large_euclidean_list) + ) + input_df["M2_smallest_euclidean_distance_hull"] = list( + map(lambda x: x / self._currPixPerMM, M2_hull_small_euclidean_list) + ) + input_df["M2_mean_euclidean_distance_hull"] = list( + map(lambda x: x / self._currPixPerMM, M2_hull_mean_euclidean_list) + ) + input_df["M2_sum_euclidean_distance_hull"] = list( + map(lambda x: x / self._currPixPerMM, M2_hull_sum_euclidean_list) + ) + input_df["Sum_euclidean_distance_hull_M1_M2"] = ( + input_df["M1_sum_euclidean_distance_hull"] + + input_df["M2_sum_euclidean_distance_hull"] + ) + + if self._debug: + print("Calculating hull variables...") + print((time.time() - start_time) * 1000, "ms") + + ########### COLLAPSED MEASURES ########################################### + + start_time = time.time() + input_df["Total_movement_centroids"] = ( + input_df["Movement_mouse_1_centroid"] + + input_df["Movement_mouse_2_centroid"] + ) + input_df["Total_movement_all_bodyparts_M1"] = ( + input_df["Movement_mouse_1_centroid"] + + input_df["Movement_mouse_1_nose"] + + input_df["Movement_mouse_1_tail_base"] + + input_df["Movement_mouse_1_left_ear"] + + input_df["Movement_mouse_1_right_ear"] + + input_df["Movement_mouse_1_lateral_left"] + + input_df["Movement_mouse_1_lateral_right"] + ) + input_df["Total_movement_all_bodyparts_M2"] = ( + input_df["Movement_mouse_2_centroid"] + + input_df["Movement_mouse_2_nose"] + + input_df["Movement_mouse_2_tail_base"] + + input_df["Movement_mouse_2_left_ear"] + + input_df["Movement_mouse_2_right_ear"] + + input_df["Movement_mouse_2_lateral_left"] + + input_df["Movement_mouse_2_lateral_right"] + ) + input_df["Total_movement_all_bodyparts_both_mice"] = ( + input_df["Total_movement_all_bodyparts_M1"] + + input_df["Total_movement_all_bodyparts_M2"] + ) + + if self._debug: + print("Collapsed measures") + print((time.time() - start_time) * 1000, "ms") + + ########### CALC ROLLING WINDOWS MEDIANS AND MEANS ########################################### + # step simplification: remove pd.DataFrame.rolling() for fixed selection. In a limited time window a rolling approach is not neccessary + start_time = time.time() + + for roll_value in self._roll_windows: + + parameter_dict = dict( + Mouse1_width="Mouse_1_width", + Mouse2_width="Mouse_2_width", + Distance="Centroid_distance", + Movement="Total_movement_centroids", + Sum_euclid_distances_hull="Sum_euclidean_distance_hull_M1_M2", + ) + for key, clm_name in parameter_dict.items(): + + clm_array = input_df[clm_name][:roll_value].to_numpy() + + currentcolname = f"{key}_mean_" + str(roll_value) + input_df[currentcolname] = np.mean(clm_array) + + currentcolname = f"{key}_median_" + str(roll_value) + input_df[currentcolname] = np.median(clm_array) + + currentcolname = f"{key}_sum_" + str(roll_value) + input_df[currentcolname] = np.sum(clm_array) + + clm_name = "euclidean_distance_hull" + clm_name2 = "euclid_distances" + for mouse in mice: + for value in ("largest", "smallest", "mean"): + + clm_array = input_df[f"M{mouse}_{value}_{clm_name}"][ + :roll_value + ].to_numpy() + + currentcolname = f"Mouse{mouse}_{value}_{clm_name2}_mean_" + str( + roll_value + ) + input_df[currentcolname] = np.mean(clm_array) + + currentcolname = f"Mouse{mouse}_{value}_{clm_name2}_median_" + str( + roll_value + ) + input_df[currentcolname] = np.median(clm_array) + + currentcolname = f"Mouse{mouse}_{value}_{clm_name2}_sum_" + str( + roll_value + ) + input_df[currentcolname] = np.sum(clm_array) + + clm_list = [ + "Total_movement_all_bodyparts_both_mice", + "Total_movement_centroids", + ] + + for clm_name in clm_list: + clm_array = input_df[clm_name][:roll_value].to_numpy() + + currentcolname = clm_name + "_mean_" + str(roll_value) + input_df[currentcolname] = np.mean(clm_array) + + currentcolname = clm_name + "_median_" + str(roll_value) + input_df[currentcolname] = np.median(clm_array) + + currentcolname = clm_name + "_sum_" + str(roll_value) + input_df[currentcolname] = np.sum(clm_array) + + parameter_dict = dict( + Nose_movement="nose", + Centroid_movement="centroid", + Tail_base_movement="tail_base", + ) + for mouse in mice: + for key, bp in parameter_dict.items(): + clm_array = input_df[f"Movement_mouse_{mouse}_{bp.lower()}"][ + :roll_value + ].to_numpy() + + currentcolname = f"{key}_M{mouse}_mean_" + str(roll_value) + input_df[currentcolname] = np.mean(clm_array) + + currentcolname = f"{key}_M{mouse}_median_" + str(roll_value) + input_df[currentcolname] = np.median(clm_array) + + currentcolname = f"{key}_M{mouse}_sum_" + str(roll_value) + input_df[currentcolname] = np.sum(clm_array) + + if self._debug: + print("Calculating rolling windows: medians, medians, and sums...") + print((time.time() - start_time) * 1000, "ms") + + ########### BODY PARTS RELATIVE TO EACH OTHER ################## + + ################# EMPETY ######################################### + + ########### ANGLES ########################################### + start_time = time.time() + input_df["Mouse_1_angle"] = input_df.apply( + lambda x: angle3pt( + x["Nose_1_x"], + x["Nose_1_y"], + x["Center_1_x"], + x["Center_1_y"], + x["Tail_base_1_x"], + x["Tail_base_1_y"], + ), + axis=1, + ) + input_df["Mouse_2_angle"] = input_df.apply( + lambda x: angle3pt( + x["Nose_2_x"], + x["Nose_2_y"], + x["Center_2_x"], + x["Center_2_y"], + x["Tail_base_2_x"], + x["Tail_base_2_y"], + ), + axis=1, + ) + input_df["Total_angle_both_mice"] = ( + input_df["Mouse_1_angle"] + input_df["Mouse_2_angle"] + ) + for roll_value in self._roll_windows: + currentcolname = "Total_angle_both_mice_" + str(roll_value) + input_df[currentcolname] = np.sum( + input_df["Total_angle_both_mice"][:roll_value].to_numpy() + ) + + if self._debug: + print("Calculating angles...") + print((time.time() - start_time) * 1000, "ms") + + ########### DEVIATIONS ########################################### + start_time = time.time() + parameter_dict = dict( + Total_movement_all_bodyparts_both_mice_deviation="Total_movement_all_bodyparts_both_mice", + Sum_euclid_distances_hull_deviation="Sum_euclidean_distance_hull_M1_M2", + M1_smallest_euclid_distances_hull_deviation="M1_smallest_euclidean_distance_hull", + M1_largest_euclid_distances_hull_deviation="M1_largest_euclidean_distance_hull", + M1_mean_euclid_distances_hull_deviation="M1_mean_euclidean_distance_hull", + Centroid_distance_deviation="Centroid_distance", + Total_angle_both_mice_deviation="Total_angle_both_mice", + Movement_mouse_1_deviation_centroid="Movement_mouse_1_centroid", + Movement_mouse_2_deviation_centroid="Movement_mouse_2_centroid", + Mouse_1_polygon_deviation="Mouse_1_poly_area", + Mouse_2_polygon_deviation="Mouse_2_poly_area", + ) + for key, value in parameter_dict.items(): + currentClm = input_df[value].to_numpy() + input_df[key] = np.mean(currentClm) - currentClm + + for roll_value in self._roll_windows: + + clm_names = ( + "Total_movement_all_bodyparts_both_mice_mean_", + "Sum_euclid_distances_hull_mean_", + "Mouse1_smallest_euclid_distances_mean_", + "Mouse1_largest_euclid_distances_mean_", + "Movement_mean_", + "Distance_mean_", + "Total_angle_both_mice_", + ) + + for part in clm_names: + currentcolname = part + str(roll_value) + currentClm = input_df[currentcolname].to_numpy() + output_array = np.mean(currentClm) - currentClm + input_df[f"{currentcolname}_deviation"] = output_array + # same values will be applied to another clm ... weird but what should i do? + # Yes I checked. Twice, actually 3 times. + input_df[f"{currentcolname}_percentile_rank"] = output_array + + if self._debug: + print("Calculating deviations...") + print((time.time() - start_time) * 1000, "ms") + + ########### PERCENTILE RANK ########################################### + + start_time = time.time() + input_df["Movement_percentile_rank"] = input_df[ + "Total_movement_centroids" + ].rank(pct=True) + input_df["Distance_percentile_rank"] = input_df["Centroid_distance"].rank( + pct=True + ) + input_df["Movement_mouse_1_percentile_rank"] = input_df[ + "Movement_mouse_1_centroid" + ].rank(pct=True) + input_df["Movement_mouse_2_percentile_rank"] = input_df[ + "Movement_mouse_1_centroid" + ].rank(pct=True) + input_df["Movement_mouse_1_deviation_percentile_rank"] = input_df[ + "Movement_mouse_1_deviation_centroid" + ].rank(pct=True) + input_df["Movement_mouse_2_deviation_percentile_rank"] = input_df[ + "Movement_mouse_2_deviation_centroid" + ].rank(pct=True) + input_df["Centroid_distance_percentile_rank"] = input_df[ + "Centroid_distance" + ].rank(pct=True) + input_df["Centroid_distance_deviation_percentile_rank"] = input_df[ + "Centroid_distance_deviation" + ].rank(pct=True) + + if self._debug: + print("Calculating percentile ranks...") + print((time.time() - start_time) * 1000, "ms") + + ########### CALCULATE STRAIGHTNESS OF POLYLINE PATH: tortuosity ########################################### + # as_strided = np.lib.stride_tricks.as_strided + start_time = time.time() + + win_size = 3 + centroidList_Mouse1_x = as_strided( + input_df.Center_1_x, + (len(input_df) - (win_size - 1), win_size), + (input_df.Center_1_x.values.strides * 2), + ) + centroidList_Mouse1_y = as_strided( + input_df.Center_1_y, + (len(input_df) - (win_size - 1), win_size), + (input_df.Center_1_y.values.strides * 2), + ) + centroidList_Mouse2_x = as_strided( + input_df.Center_2_x, + (len(input_df) - (win_size - 1), win_size), + (input_df.Center_2_x.values.strides * 2), + ) + centroidList_Mouse2_y = as_strided( + input_df.Center_2_y, + (len(input_df) - (win_size - 1), win_size), + (input_df.Center_2_y.values.strides * 2), + ) + + for k in range(self._roll_window_values_len): + start = 0 + end = start + int(self._roll_window_values[k]) + tortuosity_M1 = [] + tortuosity_M2 = [] + for y in range(len(input_df)): + tortuosity_List_M1 = [] + tortuosity_List_M2 = [] + CurrCentroidList_Mouse1_x = centroidList_Mouse1_x[start:end] + CurrCentroidList_Mouse1_y = centroidList_Mouse1_y[start:end] + CurrCentroidList_Mouse2_x = centroidList_Mouse2_x[start:end] + CurrCentroidList_Mouse2_y = centroidList_Mouse2_y[start:end] + for i in range(len(CurrCentroidList_Mouse1_x)): + currMovementAngle_mouse1 = angle3pt( + CurrCentroidList_Mouse1_x[i][0], + CurrCentroidList_Mouse1_y[i][0], + CurrCentroidList_Mouse1_x[i][1], + CurrCentroidList_Mouse1_y[i][1], + CurrCentroidList_Mouse1_x[i][2], + CurrCentroidList_Mouse1_y[i][2], + ) + currMovementAngle_mouse2 = angle3pt( + CurrCentroidList_Mouse2_x[i][0], + CurrCentroidList_Mouse2_y[i][0], + CurrCentroidList_Mouse2_x[i][1], + CurrCentroidList_Mouse2_y[i][1], + CurrCentroidList_Mouse2_x[i][2], + CurrCentroidList_Mouse2_y[i][2], + ) + tortuosity_List_M1.append(currMovementAngle_mouse1) + tortuosity_List_M2.append(currMovementAngle_mouse2) + tortuosity_M1.append(sum(tortuosity_List_M1) / (2 * math.pi)) + tortuosity_M2.append(sum(tortuosity_List_M2) / (2 * math.pi)) + start += 1 + end += 1 + currentcolname1 = str("Tortuosity_Mouse1_") + str( + self._roll_window_values[k] + ) + # currentcolname2 = str('Tortuosity_Mouse2_') + str(self._roll_window_values[k]) + input_df[currentcolname1] = tortuosity_M1 + # input_df[currentcolname2] = tortuosity_M2 + if self._debug: + print("Calculating path tortuosities...") + print((time.time() - start_time) * 1000, "ms") + + ########### CALC THE NUMBER OF LOW PROBABILITY DETECTIONS & TOTAL PROBABILITY VALUE FOR ROW########################################### + start_time = time.time() + + # SKIPPING BECAUSE DLSTREAM DOES NOT USE PROBABILITY + input_df["Sum_probabilities"] = 0.99 * self._num_bodyparts + input_df["Sum_probabilities_deviation"] = 0 + # output of ranking with all 0 + input_df["Sum_probabilities_deviation_percentile_rank"] = 0.53571 + input_df["Sum_probabilities_percentile_rank"] = 0.53571 + + input_df["Low_prob_detections_0.1"] = 0.0 + input_df["Low_prob_detections_0.5"] = 0.0 + input_df["Low_prob_detections_0.75"] = 0.0 + if self._debug: + print("Calculating pose probability scores...") + print((time.time() - start_time) * 1000, "ms") + + ########### DROP COORDINATE COLUMNS ########################################### + + input_df = input_df.reset_index(drop=True) + input_df = input_df.fillna(0) + input_df = input_df.drop(columns=["index"]) + # drop coordinates and likelyhood + # input_df = input_df.drop(columns = self._columnheaders + self._pheaders) + + return input_df + + def extract_features(self, input_array): + """Takes bp coordinates of length input_list_length and extract features. + :return extracted feature list for input to classifier" + Adapted from code provided by Simon Nilsson from Golden Lab; Main developer of SiMBA https://github.com/sgoldenlab/simba""" + + if len(input_array) == self.input_array_length: + start_time = time.time() + input_df = self.convert_pandas(input_array) + features = self.extract_features_simba14bp(input_df) + + if self._debug: + print( + "Full feature extraction time: ", (time.time() - start_time) * 1000 + ) + return features.to_numpy() + + else: + return None + + def get_currPixPerMM(self): + return self.currPixPerMM + + def get_input_array_length(self): + return self.input_array_length + + def set_input_array_length(self, new_input_array_length): + self.input_array_length = new_input_array_length + + +"""For fast simba extraction""" + + +### EUCLIDIAN DISTANCES IN SINGLE FRAMES +@jit(nopython=True, cache=True) +def EuclidianDistCalc(inArr): + """provided by Simon Nilsson from Golden Lab; Main developer of SiMBA https://github.com/sgoldenlab/simba""" + Mouse_1_nose_to_tail = int( + np.sqrt( + ((inArr[-1][4] - inArr[-1][12]) ** 2 + (inArr[-1][5] - inArr[-1][13]) ** 2) + ) + ) + Mouse_2_nose_to_tail = int( + np.sqrt( + ( + (inArr[-1][18] - inArr[-1][26]) ** 2 + + (inArr[-1][19] - inArr[-1][27]) ** 2 + ) + ) + ) + Mouse_1_Ear_distance = int( + np.sqrt( + ((inArr[-1][0] - inArr[-1][2]) ** 2 + (inArr[-1][1] - inArr[-1][3]) ** 2) + ) + ) + Centroid_distance = int( + np.sqrt( + ((inArr[-1][6] - inArr[-1][20]) ** 2 + (inArr[-1][7] - inArr[-1][21]) ** 2) + ) + ) + Nose_to_nose_distance = int( + np.sqrt( + ((inArr[-1][4] - inArr[-1][18]) ** 2 + (inArr[-1][5] - inArr[-1][19]) ** 2) + ) + ) + M1_Nose_to_M2_lat_left = int( + np.sqrt( + ((inArr[-1][4] - inArr[-1][22]) ** 2 + (inArr[-1][5] - inArr[-1][23]) ** 2) + ) + ) + M1_Nose_to_M2_lat_right = int( + np.sqrt( + ((inArr[-1][4] - inArr[-1][24]) ** 2 + (inArr[-1][5] - inArr[-1][25]) ** 2) + ) + ) + M2_Nose_to_M1_lat_left = int( + np.sqrt( + ((inArr[-1][18] - inArr[-1][8]) ** 2 + (inArr[-1][19] - inArr[-1][9]) ** 2) + ) + ) + M2_Nose_to_M1_lat_right = int( + np.sqrt( + ( + (inArr[-1][18] - inArr[-1][10]) ** 2 + + (inArr[-1][19] - inArr[-1][11]) ** 2 + ) + ) + ) + M1_Nose_to_M2_tail_base = int( + np.sqrt( + ((inArr[-1][4] - inArr[-1][26]) ** 2 + (inArr[-1][5] - inArr[-1][27]) ** 2) + ) + ) + M2_Nose_to_M1_tail_base = int( + np.sqrt( + ( + (inArr[-1][18] - inArr[-1][12]) ** 2 + + (inArr[-1][19] - inArr[-1][13]) ** 2 + ) + ) + ) + + return np.array( + [ + Mouse_1_nose_to_tail, + Mouse_2_nose_to_tail, + Mouse_1_Ear_distance, + Centroid_distance, + Nose_to_nose_distance, + M1_Nose_to_M2_lat_left, + M1_Nose_to_M2_lat_right, + M2_Nose_to_M1_lat_left, + M2_Nose_to_M1_lat_right, + M1_Nose_to_M2_tail_base, + M2_Nose_to_M1_tail_base, + ] + ) + + +### EUCLIDEAN DISTANCES BETWEEN CENTROIDS IN ROLLING WINDOWS +@jit(nopython=True, cache=True) +def distancesBetweenBps(bpArray, bp): + # TODO: adapt to flexible window + """provided by Simon Nilsson from Golden Lab; Main developer of SiMBA https://github.com/sgoldenlab/simba""" + + frames2Process = bpArray.shape[0] + bps2process = int(bpArray.shape[1] / 2) + outputArray = np.zeros((frames2Process, bps2process - 1)) + for frame in range(frames2Process): + outputArray[frame] = np.sqrt( + (bpArray[frame][0] - bpArray[frame][2]) ** 2 + + (bpArray[frame][1] - bpArray[frame][3]) ** 2 + ) + Distance_median_2, Distance_mean_2, Distance_sum_2 = ( + int(np.median(outputArray)), + int(np.mean(outputArray)), + int(np.sum(outputArray)), + ) + if bp == "centroid": + # currently hardcoded to 15 frames + msArr200, msArr166, msArr133, msArr66 = ( + outputArray[8:15], + outputArray[10:15], + outputArray[11:15], + outputArray[13:15], + ) + Distance_median_5, Distance_mean_5, Distance_sum_5 = ( + int(np.median(msArr200)), + int(np.mean(msArr200)), + int(np.sum(msArr200)), + ) + Distance_median_6, Distance_mean_6, Distance_sum_6 = ( + int(np.median(msArr166)), + int(np.mean(msArr166)), + int(np.sum(msArr166)), + ) + Distance_median_7, Distance_mean_7, Distance_sum_7 = ( + int(np.median(msArr133)), + int(np.mean(msArr133)), + int(np.sum(msArr133)), + ) + Distance_median_15, Distance_mean_15, Distance_sum_15 = ( + int(np.median(msArr66)), + int(np.mean(msArr66)), + int(np.sum(msArr66)), + ) + + return np.array( + [ + Distance_median_2, + Distance_mean_2, + Distance_sum_2, + Distance_median_5, + Distance_mean_5, + Distance_sum_5, + Distance_median_6, + Distance_mean_6, + Distance_sum_6, + Distance_median_7, + Distance_mean_7, + Distance_sum_7, + Distance_median_15, + Distance_mean_15, + Distance_sum_15, + ] + ) + if bp == "width": + return np.array([Distance_median_2, Distance_mean_2, Distance_sum_2]) + + +### EUCLIDEAN DISTANCES BETWEEN BODY-PARTS WITHIN EACH ANIMALS HULL +@jit(nopython=True, cache=True) +def bpDistancesInHull(animal1arraySplit): + """provided by Simon Nilsson from Golden Lab; Main developer of SiMBA https://github.com/sgoldenlab/simba""" + + frames2process = animal1arraySplit.shape[0] + bps2process = animal1arraySplit[1].shape[0] + outputArray = np.zeros((frames2process, bps2process)) + for frame in range(frames2process): + for bodypart in range(bps2process): + for otherBps in [x for x in range(bps2process) if x != bodypart]: + outputArray[frame][bodypart] = np.sqrt( + ( + animal1arraySplit[frame][bodypart][0] + - animal1arraySplit[frame][otherBps][0] + ) + ** 2 + + ( + animal1arraySplit[frame][bodypart][1] + - animal1arraySplit[frame][otherBps][1] + ) + ** 2 + ) + flattenedOutput = outputArray.flatten() + Mean_euclid_distances_median_2 = np.median(flattenedOutput) + Mean_euclid_distances_mean_2 = np.mean(flattenedOutput) + return np.array([Mean_euclid_distances_median_2, Mean_euclid_distances_mean_2]) + + +### TOTAL MOVEMENT OF ALL ANIMALS IN ROLLING WINDOWS +@jit(nopython=True, cache=True) +def TotalMovementBodyparts(arrayConcat_Animal1, arrayConcat_Animal2): + """provided by Simon Nilsson from Golden Lab; Main developer of SiMBA https://github.com/sgoldenlab/simba""" + + frames2process = arrayConcat_Animal2.shape[0] + bps2process = int(arrayConcat_Animal2.shape[1] / 2) + outputArray_animal_1, outputArray_animal_2 = ( + np.zeros((frames2process, bps2process)), + np.zeros((frames2process, bps2process)), + ) + for frame in range(frames2process): + for bp_current, bp_shifted in zip(range(0, 7), range(7, 14)): + outputArray_animal_1[frame][bp_current] = np.sqrt( + ( + arrayConcat_Animal1[frame][bp_current][0] + - arrayConcat_Animal1[frame][bp_shifted][0] + ) + ** 2 + + ( + arrayConcat_Animal1[frame][bp_current][1] + - arrayConcat_Animal1[frame][bp_shifted][1] + ) + ** 2 + ) + outputArray_animal_2[frame][bp_current] = np.sqrt( + ( + arrayConcat_Animal2[frame][bp_current][0] + - arrayConcat_Animal2[frame][bp_shifted][0] + ) + ** 2 + + ( + arrayConcat_Animal2[frame][bp_current][1] + - arrayConcat_Animal2[frame][bp_shifted][1] + ) + ** 2 + ) + sumAnimal1, sumAnimal2 = np.sum(outputArray_animal_1, axis=1), np.sum( + outputArray_animal_2, axis=1 + ) + sumConcat = np.concatenate((sumAnimal1, sumAnimal2)) + Total_movement_all_bodyparts_both_mice_median_2 = int(np.median(sumConcat)) + Total_movement_all_bodyparts_both_mice_mean_2 = int(np.mean(sumConcat)) + Total_movement_all_bodyparts_both_mice_sum_2 = int(np.sum(sumConcat)) + last200msArrayAnimal1, last200msArrayAnimal2 = sumAnimal1[9:13], sumAnimal2[9:13] + sumConcat200ms = np.concatenate((last200msArrayAnimal1, last200msArrayAnimal2)) + Total_movement_all_bodyparts_both_mice_mean_5 = int(np.mean(sumConcat200ms)) + Total_movement_all_bodyparts_both_mice_sum_5 = int(np.sum(sumConcat200ms)) + + return ( + np.array( + [ + Total_movement_all_bodyparts_both_mice_median_2, + Total_movement_all_bodyparts_both_mice_mean_2, + Total_movement_all_bodyparts_both_mice_sum_2, + Total_movement_all_bodyparts_both_mice_mean_5, + Total_movement_all_bodyparts_both_mice_sum_5, + ] + ), + outputArray_animal_1, + outputArray_animal_2, + ) + + +### MOVEMENTS OF INDIVIDUAL BODY-PARTS +@jit(nopython=True, cache=True) +def singleAnimalBpMovements(tail_1, tail_2, center_1, center_2, nose_1, nose_2): + """provided by Simon Nilsson from Golden Lab; Main developer of SiMBA https://github.com/sgoldenlab/simba""" + + Tail_base_movement_M1_median_2, Tail_base_movement_M2_median_2 = ( + int(np.median(tail_1)), + int(np.median(tail_2)), + ) + Tail_base_movement_M2_mean_2, Tail_base_movement_M2_sum_2 = ( + int(np.mean(tail_2)), + int(np.sum(tail_2)), + ) + Centroid_movement_M1_mean_2, Centroid_movement_M1_sum_2 = ( + int(np.mean(center_1)), + int(np.sum(center_1)), + ) + Centroid_movement_M2_mean_2, Centroid_movement_M2_sum_2 = ( + int(np.mean(center_2)), + int(np.sum(center_2)), + ) + Nose_movement_M1_median_2, Nose_movement_M1_mean_2 = int(np.median(nose_1)), int( + np.mean(nose_1) + ) + Nose_movement_M1_sum_2, Nose_movement_M2_mean_2, Nose_movement_M2_sum_2 = ( + int(np.sum(nose_1)), + int(np.mean(nose_2)), + int(np.sum(nose_2)), + ) + return np.array( + [ + Tail_base_movement_M1_median_2, + Tail_base_movement_M2_median_2, + Tail_base_movement_M2_mean_2, + Tail_base_movement_M2_sum_2, + Centroid_movement_M1_mean_2, + Centroid_movement_M1_sum_2, + Centroid_movement_M2_mean_2, + Centroid_movement_M2_sum_2, + Nose_movement_M1_median_2, + Nose_movement_M1_mean_2, + Nose_movement_M1_sum_2, + Nose_movement_M2_mean_2, + Nose_movement_M2_sum_2, + ] + ) + + +class SimbaFeatureExtractor: + """Feature extraction module for integration in behavior trigger, takes list of postures as input + and calculates features to pass to classifier. Features and classifier have to match! + Designed to work with Simba https://github.com/sgoldenlab/simba""" + + def __init__(self, input_array_length): + self.currPixPerMM = PIXPERMM + self.input_array_length = input_array_length + # TODO: Collect bodypart position in input array from skeleton automatically (this will make this much easier!) + + def get_currPixPerMM(self): + return self.currPixPerMM + + def get_input_array_length(self): + return self.input_array_length + + def set_input_array_length(self, new_input_array_length): + self.input_array_length = new_input_array_length + + def extract_features(self, input_array): + """Takes bp coordinates of length input_list_length and extract features. + :return extracted feature list for input to classifier" + Adapted from code provided by Simon Nilsson from Golden Lab; Main developer of SiMBA https://github.com/sgoldenlab/simba""" + + def append2featureList(featureList, measures2add): + featureList.extend(measures2add) + return featureList + + featureList = [] + if len(input_array) == self.input_array_length: + input_array = np.array(input_array).astype(int) + """Start extracting features""" + ### EUCLIDIAN DISTANCES IN SINGLE FRAMES + distanceMeasures = EuclidianDistCalc(input_array) + featureList = append2featureList(featureList, list(distanceMeasures)) + + ### EUCLIDEAN DISTANCES BETWEEN CENTROIDS/width IN ROLLING WINDOWS + centroidM1x, centroidM1y, centroidM2x, centroidM2y = ( + input_array[:, [6]], + input_array[:, [7]], + input_array[:, [20]], + input_array[:, [21]], + ) + bpArray = np.concatenate( + (centroidM1x, centroidM1y, centroidM2x, centroidM2y), 1 + ) + distancefeatures = distancesBetweenBps(bpArray, "centroid") + featureList = append2featureList(featureList, list(distancefeatures)) + + latLeftM1x, latLeftM1y, latRightM1x, latRightM1y = ( + input_array[:, [8]], + input_array[:, [9]], + input_array[:, [10]], + input_array[:, [11]], + ) + bpArray = np.concatenate( + (latLeftM1x, latLeftM1y, latRightM1x, latRightM1y), 1 + ) + animal_1_Widths = distancesBetweenBps(bpArray, "width") + featureList = append2featureList(featureList, list(animal_1_Widths)) + + latLeftM2x, latLeftM2y, latRightM2x, latRightM2y = ( + input_array[:, [22]], + input_array[:, [23]], + input_array[:, [24]], + input_array[:, [25]], + ) + bpArray = np.concatenate( + (latLeftM2x, latLeftM2y, latRightM2x, latRightM2y), 1 + ) + animal_2_Widths = distancesBetweenBps(bpArray, "width") + featureList = append2featureList(featureList, list(animal_2_Widths)) + + ### EUCLIDEAN DISTANCES BETWEEN BODY-PARTS WITHIN EACH ANIMALS HULL + animal1array, animal2array = input_array[:, 0:14], input_array[:, 14:28] + animal1arraySplit, animal2arraySplit = ( + animal1array.reshape(-1, 7, 2), + animal2array.reshape(-1, 7, 2), + ) + hullVars1 = bpDistancesInHull(animal1arraySplit) + hullVars2 = bpDistancesInHull(animal2arraySplit) + featureList = append2featureList(featureList, list(hullVars1)) + featureList = append2featureList(featureList, list(hullVars2)) + + ### TOTAL MOVEMENT OF ALL ANIMALS IN ROLLING WINDOWS + animalSplitArray = input_array.reshape(-1, 7, 2) + animal_1_SplitArray = animalSplitArray[0::2] + animal_2_SplitArray = animalSplitArray[1::2] + shiftedSplitArray_animal_1 = np.roll(animal_1_SplitArray, -1, axis=0) + shiftedSplitArray_animal_2 = np.roll(animal_2_SplitArray, -1, axis=0) + shiftedSplitArray_animal_1, shiftedSplitArray_animal_2 = ( + shiftedSplitArray_animal_1[:-1].copy(), + shiftedSplitArray_animal_2[:-1].copy(), + ) + animal_1_SplitArray, animal_2_SplitArray = ( + animal_1_SplitArray[:-1].copy(), + animal_2_SplitArray[:-1].copy(), + ) + arrayConcat_Animal1 = np.concatenate( + (animal_1_SplitArray, shiftedSplitArray_animal_1), 1 + ) + arrayConcat_Animal2 = np.concatenate( + (animal_2_SplitArray, shiftedSplitArray_animal_2), 1 + ) + ( + totalMovementVars, + outputArray_animal_1, + outputArray_animal_2, + ) = TotalMovementBodyparts(arrayConcat_Animal1, arrayConcat_Animal2) + featureList = append2featureList(featureList, list(totalMovementVars)) + + tail1, tail2 = outputArray_animal_1[:, [6]], outputArray_animal_2[:, [6]] + center1, center2 = ( + outputArray_animal_1[:, [3]], + outputArray_animal_2[:, [3]], + ) + nose1, nose2 = outputArray_animal_1[:, [0]], outputArray_animal_2[:, [0]] + indBpMovs = singleAnimalBpMovements( + tail1, tail2, center1, center2, nose1, nose2 + ) + featureList = append2featureList(featureList, list(indBpMovs)) + + featureList = [x / self.currPixPerMM for x in featureList] + featureList = np.array(featureList) + featureList = featureList.reshape(1, -1) + + return featureList + + else: + return None + + +class BsoidFeatureExtractor: + """Feature extraction module for integration in behavior trigger, takes list of postures as input + and calculates features to pass to classifier. Features and classifier have to match! + Designed to work with BSOID; https://github.com/YttriLab/B-SOID""" + + def __init__(self, input_array_length): + self.input_array_length = input_array_length + self._fps = FRAMERATE + + def get_currPixPerMM(self): + return None + + def get_input_array_length(self): + return self.input_array_length + + def set_input_array_length(self, new_input_array_length): + self.input_array_length = new_input_array_length + + def extract_features(self, input_array): + """ + Extracts features based on (x,y) positions + :param input_array: list of poses + :return f_10fps: 2D array, extracted features + Adapted from BSOID; https://github.com/YttriLab/B-SOID + """ + + def boxcar_center(a, n): + a1 = pd.Series(a) + moving_avg = np.array( + a1.rolling(window=n, min_periods=1, center=True).mean() + ) + + return moving_avg + + def adp_filt_pose(pose_estimation): + """Adapted from adp_filt function in BSOID""" + currdf = np.array(pose_estimation) + datax = currdf[:, :, 0] + datay = currdf[:, :, 1] + # TODO: Adapt filter to work without workaround and skeleton + # data_lh = currdf[:, :, 2] + data_lh = np.full_like(datax, 0.9) + currdf_filt = np.zeros((datax.shape[0], (datax.shape[1]) * 2)) + perc_rect = [] + for i in range(data_lh.shape[1]): + perc_rect.append(0) + for x in range(data_lh.shape[1]): + a, b = np.histogram(data_lh[1:, x].astype(np.float)) + rise_a = np.where(np.diff(a) >= 0) + if rise_a[0][0] > 1: + llh = b[rise_a[0][0]] + else: + llh = b[rise_a[0][1]] + data_lh_float = data_lh[:, x].astype(np.float) + perc_rect[x] = np.sum(data_lh_float < llh) / data_lh.shape[0] + currdf_filt[0, (2 * x) : (2 * x + 2)] = np.hstack( + [datax[0, x], datay[0, x]] + ) + for i in range(1, data_lh.shape[0]): + if data_lh_float[i] < llh: + currdf_filt[i, (2 * x) : (2 * x + 2)] = currdf_filt[ + i - 1, (2 * x) : (2 * x + 2) + ] + else: + currdf_filt[i, (2 * x) : (2 * x + 2)] = np.hstack( + [datax[i, x], datay[i, x]] + ) + currdf_filt = np.array(currdf_filt) + currdf_filt = currdf_filt.astype(np.float) + return currdf_filt, perc_rect + + if len(input_array) == self.input_array_length: + data, p_sub_threshold = adp_filt_pose(input_array) + data = np.array([data]) + win_len = np.int(np.round(0.05 / (1 / self._fps)) * 2 - 1) + feats = [] + for m in range(len(data)): # 1 + dataRange = len(data[m]) # 5 + dxy_r = [] + dis_r = [] + for r in range(dataRange): # 0-4 + if r < dataRange - 1: + dis = [] + for c in range(0, data[m].shape[1], 2): # 0-17, 2 + dis.append( + np.linalg.norm( + data[m][r + 1, c : c + 2] - data[m][r, c : c + 2] + ) + ) + dis_r.append(dis) + dxy = [] + for i, j in itertools.combinations( + range(0, data[m].shape[1], 2), 2 + ): # 0-17, 2 + dxy.append(data[m][r, i : i + 2] - data[m][r, j : j + 2]) + dxy_r.append(dxy) + dis_r = np.array(dis_r) + dxy_r = np.array(dxy_r) + dis_smth = [] + dxy_eu = np.zeros([dataRange, dxy_r.shape[1]]) # 5,36 + ang = np.zeros([dataRange - 1, dxy_r.shape[1]]) + dxy_smth = [] + ang_smth = [] + for l in range(dis_r.shape[1]): # 0-8 + dis_smth.append(boxcar_center(dis_r[:, l], win_len)) + for k in range(dxy_r.shape[1]): # 0-35 + for kk in range(dataRange): # 0 + dxy_eu[kk, k] = np.linalg.norm(dxy_r[kk, k, :]) + if kk < dataRange - 1: + b_3d = np.hstack([dxy_r[kk + 1, k, :], 0]) + a_3d = np.hstack([dxy_r[kk, k, :], 0]) + c = np.cross(b_3d, a_3d) + ang[kk, k] = np.dot( + np.dot(np.sign(c[2]), 180) / np.pi, + math.atan2( + np.linalg.norm(c), + np.dot(dxy_r[kk, k, :], dxy_r[kk + 1, k, :]), + ), + ) + dxy_smth.append(boxcar_center(dxy_eu[:, k], win_len)) + ang_smth.append(boxcar_center(ang[:, k], win_len)) + dis_smth = np.array(dis_smth) + dxy_smth = np.array(dxy_smth) + ang_smth = np.array(ang_smth) + feats.append(np.vstack((dxy_smth[:, 1:], ang_smth, dis_smth))) + + f_10fps = [] + for n in range(0, len(feats)): + feats1 = np.zeros(len(data[n])) + for s in range(math.floor(self._fps / 10)): + for k in range( + round(self._fps / 10) + s, + len(feats[n][0]), + round(self._fps / 10), + ): + if k > round(self._fps / 10) + s: + feats1 = np.concatenate( + ( + feats1.reshape(feats1.shape[0], feats1.shape[1]), + np.hstack( + ( + np.mean( + ( + feats[n][ + 0 : dxy_smth.shape[0], + range( + k - round(self._fps / 10), k + ), + ] + ), + axis=1, + ), + np.sum( + ( + feats[n][ + dxy_smth.shape[0] : feats[ + n + ].shape[0], + range( + k - round(self._fps / 10), k + ), + ] + ), + axis=1, + ), + ) + ).reshape(len(feats[0]), 1), + ), + axis=1, + ) + else: + feats1 = np.hstack( + ( + np.mean( + ( + feats[n][ + 0 : dxy_smth.shape[0], + range(k - round(self._fps / 10), k), + ] + ), + axis=1, + ), + np.sum( + ( + feats[n][ + dxy_smth.shape[0] : feats[n].shape[0], + range(k - round(self._fps / 10), k), + ] + ), + axis=1, + ), + ) + ).reshape(len(feats[0]), 1) + f_10fps.append(feats1) + return f_10fps + + else: + return None diff --git a/experiments/custom/stimulus_process.py b/experiments/custom/stimulus_process.py index 5e8ae09..8e02de7 100644 --- a/experiments/custom/stimulus_process.py +++ b/experiments/custom/stimulus_process.py @@ -16,6 +16,7 @@ withdraw_liqreward, DigitalModDevice, ) +from experiments.utils.gpio_control import DigitalArduinoDevice import random @@ -76,15 +77,19 @@ def get_start_time(self): def example_protocol_run(condition_q: mp.Queue): current_trial = None # dmod_device = DigitalModDevice('Dev1/PFI0') + # led_machine = DigitalArduinoDevice("COM5") while True: + # if no protocol is selected, running default picture (background) if condition_q.full(): current_trial = condition_q.get() if current_trial is not None: show_visual_stim_img(type=current_trial, name="DlStream") # dmod_device.toggle() + # led_machine.turn_on() else: show_visual_stim_img(name="DlStream") # dmod_device.turn_off() + # led_machine.turn_off() if cv2.waitKey(1) & 0xFF == ord("q"): break diff --git a/experiments/custom/triggers.py b/experiments/custom/triggers.py index f853232..c179fd8 100644 --- a/experiments/custom/triggers.py +++ b/experiments/custom/triggers.py @@ -6,15 +6,19 @@ Licensed under GNU General Public License v3.0 """ - +from utils.poser import transform_2pose from utils.analysis import ( angle_between_vectors, calculate_distance, EllipseROI, RectangleROI, ) -from utils.configloader import RESOLUTION +from utils.configloader import RESOLUTION, TIME_WINDOW, FRAMERATE from collections import deque +from experiments.custom.featureextraction import ( + SimbaFeatureExtractor, + BsoidFeatureExtractor, +) import numpy as np import time @@ -664,3 +668,523 @@ def check_skeleton(self, skeleton: dict): response = (result, response_body) return response + + +"""Behavior classifier trigger""" + + +class SimbaThresholdBehaviorPoolTrigger: + """ + Trigger to check if animal's behavior is classified as specific motif above threshold probability. + Uses processing pool with feature extraction and classification in the same process + + """ + + def __init__(self, prob_threshold: float, class_process_pool, debug: bool = False): + """ + Initialising trigger with following parameters: + :param float prob_threshold: threshold probability of prediction that is returned by classifier and should be used as trigger. + If you plan to use the classifier for multiple trial triggers in the same experiment with different thresholds. We recommend setting up the + trigger_probability during check_skeleton + :param class_process_pool: list of dictionaries with keys process: mp.Process, input: mp.queue, output: mp.queue; + used for lossless frame-by-frame classification + + """ + self._trigger_threshold = prob_threshold + self._process_pool = class_process_pool + self._last_prob = 0.0 + self._feature_id = 0 + self._center = None + self._debug = debug + self._skeleton = None + self._time_window_len = TIME_WINDOW + # feature extraction will happen in the classification process + # self.feat_extractor = None + self._time_window = deque(maxlen=self._time_window_len) + + def fill_time_window(self, skeleton: dict): + """Transforms skeleton input into flat numpy array of coordinates to pass to feature extraction""" + # todo: remove bodyparts that are not used automatically + key_selection = {"0_tail_tip", "1_tail_tip"} + skeleton_selection = { + k: skeleton[k] for k in skeleton.keys() if k not in key_selection + } + flat_values = transform_2pose(skeleton_selection).flatten() + + # if not enough animals are present, padd the rest with default value "0,0" + # TODO: Padd automatically; remove hardcoding + if flat_values.shape[0] < 28: + flat_values = np.pad( + flat_values, + (0, 28 - flat_values.shape[0]), + "constant", + constant_values=0, + ) + # this appends the new row to the deque time_window, which will drop the "oldest" entry due to a maximum + # length of time_window_len + self._time_window.append(flat_values) + + def check_skeleton(self, skeleton, target_prob: float = None): + """ + Checking skeleton for trigger, will pass skeleton window to classifier if window length is reached and + collect skeletons otherwise + :param skeleton: a skeleton dictionary, returned by calculate_skeletons() from poser file + :param target_prob: optional, overwrites self._trigger_prob with target probability. this is supposed to enable the + set up of different trials (with different motif thresholds) in the experiment without the necessaty to init to + classifiers: default None + :return: response, a tuple of result (bool) and response body + Response body is used for plotting and outputting results to trials dataframes + """ + self.fill_time_window(skeleton) + + """Checks if necessary time window was collected and passes it to classifier""" + # if enough postures where collected + + if len(self._time_window) == self._time_window_len: + # if the last classification is done and was taken + self._feature_id += 1 + self._process_pool.pass_time_window( + (self._time_window, self._feature_id), debug=self._debug + ) + # check if a process from the pool is done with the result + result, feature_id = self._process_pool.get_result(debug=self._debug) + if result is not None: + self._last_prob = result + # else: + # self._last_prob = 0.0 + + if target_prob is not None: + self._trigger_threshold = target_prob + # choosing a point to draw near the skeleton + self._center = (50, 50) + result = False + text = "Current probability: {:.2f}".format(self._last_prob) + + if self._trigger_threshold <= self._last_prob: + result = True + text = "Motif matched: {:.2f}".format(self._last_prob) + + color = (0, 255, 0) if result else (0, 0, 255) + response_body = { + "plot": {"text": dict(text=text, org=self._center, color=color)} + } + response = (result, response_body) + return response + + def get_trigger_threshold(self): + return self._trigger_threshold + + def get_last_prob(self): + return self._last_prob + + def get_time_window_len(self): + return self._time_window_len + + +class BsoidClassBehaviorPoolTrigger: + """ + Trigger to check if animal's behavior is classified as specific motif with BSOID trained classifier. + Uses processing pool with feature extraction and classification in the same process + """ + + def __init__(self, target_class: int, class_process_pool, debug: bool = False): + """ + Initialising trigger with following parameters: + :param int target_class: target classification category that should be used as trigger. Must match "Group" number of cluster in BSOID. + If you plan to use the classifier for multiple trial triggers in the same experiment with different thresholds. We recommend setting up the + target_class during check_skeleton + :param class_process_pool: list of dictionaries with keys process: mp.Process, input: mp.queue, output: mp.queue; + used for lossless frame-by-frame classification + + """ + self._trigger = target_class + self._process_pool = class_process_pool + self._last_result = [0] + self._feature_id = 0 + self._center = None + self._debug = debug + self._skeleton = None + self._time_window_len = TIME_WINDOW + # not used in this version + self.feat_extractor = None + self._time_window = deque(maxlen=self._time_window_len) + + def fill_time_window(self, skeleton): + from utils.poser import transform_2pose + + pose = transform_2pose(skeleton) + self._time_window.appendleft(pose) + + def check_skeleton(self, skeleton, target_class: int = None): + """ + Checking skeleton for trigger, will pass skeleton window to classifier if window length is reached and + collect skeletons otherwise + :param skeleton: a skeleton dictionary, returned by calculate_skeletons() from poser file + :param target_class: optional, overwrites self._trigger with target probability. this is supposed to enable the + set up of different trials (with different motifs/categories) in the experiment without the necessity to init to + classifiers: default None + :return: response, a tuple of result (bool) and response body + Response body is used for plotting and outputting results to trials dataframes + """ + self.fill_time_window(skeleton) + + """Checks if necessary time window was collected and passes it to classifier""" + if len(self._time_window) == self._time_window_len: + self._feature_id += 1 + self._process_pool.pass_time_window( + (self._time_window, self._feature_id), debug=self._debug + ) + # check if a process from the pool is done with the result + clf_result, feature_id = self._process_pool.get_result(debug=self._debug) + if clf_result is not None: + self._last_result = clf_result[0] + if target_class is not None: + self._trigger = target_class + # choosing a point to draw near the skeleton + self._center = skeleton[list(skeleton.keys())[0]] + # self._center = (50,50) + result = False + # text = 'Current probability: {:.2f}'.format(self._last_prob) + text = "Current Class: {}".format(self._last_result) + + if self._last_result[0] == self._trigger: + result = True + text = "Motif matched: {}".format(self._last_result) + + color = (0, 255, 0) if result else (0, 0, 255) + response_body = { + "plot": {"text": dict(text=text, org=self._center, color=color)} + } + response = (result, response_body) + return response + + def get_trigger_threshold(self): + return self._trigger + + def get_last_prob(self): + return self._last_prob + + def get_time_window_len(self): + return self._time_window_len + + +class SimbaThresholdBehaviorPoolTrigger_old: + """ + THIS VERSION HAS A SEPERATE FEATURE EXTRACTION WHICH WILL AFFECT OVERALL PERFORMANCE + Trigger to check if animal's behavior is classified as specific motif above threshold probability. + """ + + def __init__(self, prob_threshold: float, class_process_pool, debug: bool = False): + """ + Initialising trigger with following parameters: + :param float prob_threshold: threshold probability of prediction that is returned by classifier and should be used as trigger. + If you plan to use the classifier for multiple trial triggers in the same experiment with different thresholds. We recommend setting up the + trigger_probability during check_skeleton + :param class_process_pool: list of dictionaries with keys process: mp.Process, input: mp.queue, output: mp.queue; + used for lossless frame-by-frame classification + + """ + self._trigger_threshold = prob_threshold + self._process_pool = class_process_pool + self._last_prob = 0.0 + self._feature_id = 0 + self._center = None + self._debug = debug + self._skeleton = None + self._time_window_len = TIME_WINDOW + self.feat_extractor = SimbaFeatureExtractor( + input_array_length=self._time_window_len + ) + self._time_window = deque(maxlen=self._time_window_len) + + def fill_time_window(self, skeleton: dict): + """Transforms skeleton input into flat numpy array of coordinates to pass to feature extraction""" + # todo: remove bodyparts that are not used automatically + key_selection = {"0_tail_tip", "1_tail_tip"} + skeleton_selection = { + k: skeleton[k] for k in skeleton.keys() if k not in key_selection + } + flat_values = transform_2pose(skeleton_selection).flatten() + + # if not enough animals are present, padd the rest with default value "0,0" + # TODO: Padd automatically; remove hardcoding + if flat_values.shape[0] < 28: + flat_values = np.pad( + flat_values, + (0, 28 - flat_values.shape[0]), + "constant", + constant_values=0, + ) + # this appends the new row to the deque time_window, which will drop the "oldest" entry due to a maximum + # length of time_window_len + self._time_window.append(flat_values) + + def check_skeleton(self, skeleton, target_prob: float = None): + """ + Checking skeleton for trigger, will pass skeleton window to classifier if window length is reached and + collect skeletons otherwise + :param skeleton: a skeleton dictionary, returned by calculate_skeletons() from poser file + :param target_prob: optional, overwrites self._trigger_prob with target probability. this is supposed to enable the + set up of different trials (with different motif thresholds) in the experiment without the necessaty to init to + classifiers: default None + :return: response, a tuple of result (bool) and response body + Response body is used for plotting and outputting results to trials dataframes + """ + self.fill_time_window(skeleton) + f_extract_output = None + """Checks if necessary time window was collected and passes it to classifier""" + if len(self._time_window) == self._time_window_len: + start_time = time.time() + f_extract_output = self.feat_extractor.extract_features(self._time_window) + if self._debug: + end_time = time.time() + print( + "Feature extraction time: {:.2f} msec".format( + (end_time - start_time) * 1000 + ) + ) + # if enough postures where collected and their features extracted + if f_extract_output is not None: + # if the last classification is done and was taken + self._feature_id += 1 + self._process_pool.pass_features( + (f_extract_output, self._feature_id), debug=self._debug + ) + # check if a process from the pool is done with the result + result, feature_id = self._process_pool.get_result(debug=self._debug) + if result is not None: + self._last_prob = result + # else: + # self._last_prob = 0.0 + + if target_prob is not None: + self._trigger_threshold = target_prob + # choosing a point to draw near the skeleton + self._center = (50, 50) + result = False + text = "Current probability: {:.2f}".format(self._last_prob) + + if self._trigger_threshold <= self._last_prob: + result = True + text = "Motif matched: {:.2f}".format(self._last_prob) + + color = (0, 255, 0) if result else (0, 0, 255) + response_body = { + "plot": {"text": dict(text=text, org=self._center, color=color)} + } + response = (result, response_body) + return response + + def get_trigger_threshold(self): + return self._trigger_threshold + + def get_last_prob(self): + return self._last_prob + + def get_time_window_len(self): + return self._time_window_len + + +class BsoidClassBehaviorTrigger: + """ + THIS VERSION HAS A SEPERATE FEATURE EXTRACTION WHICH CAN AFFECT OVERALL PERFORMANCE AND USES ONLY 1 CLASSIFIER AT A TIME + WE RECOMMEND USING THE POOL VERSION. + If you only want/need to run 1 classifier, just set pool_size to 1 + + Trigger to check if animal's behavior is classified as specific motif with BSOID trained classifier. + """ + + def __init__(self, target_class: int, path_to_sav: str, debug: bool = False): + """ + Initialising trigger with following parameters: + :param int target_class: target classification category that should be used as trigger. Must match "Group" number of cluster in BSOID. + If you plan to use the classifier for multiple trial triggers in the same experiment with different thresholds. We recommend setting up the + target_class during check_skeleton + :param str path_to_sav: path to saved classifier, will be passed to classifier module + + """ + self._trigger = target_class + self._last_result = [0] + self._center = None + self._debug = debug + self._skeleton = None + self._classifier, self._time_window_len = self._init_classifier( + path_to_sav + ) # initialize classifier + self.feat_extractor = BsoidFeatureExtractor( + self._time_window_len, fps=FRAMERATE + ) + self._time_window = deque(maxlen=self._time_window_len) + + @staticmethod + def _init_classifier(path_to_sav): + from experiments.custom.classifier import BsoidClassifier + + """Put your classifier of choice in here""" + classifier = BsoidClassifier(path_to_clf=path_to_sav) + win_len = classifier.get_win_len() + return classifier, win_len + + def fill_time_window(self, skeleton): + from utils.poser import transform_2pose + + pose = transform_2pose(skeleton) + self._time_window.appendleft(pose) + + def check_skeleton(self, skeleton, trigger: float = None): + """ + Checking skeleton for trigger, will pass skeleton window to classifier if window length is reached and + collect skeletons otherwise + :param skeleton: a skeleton dictionary, returned by calculate_skeletons() from poser file + :param trigger: optional, overwrites self._trigger with target probability. this is supposed to enable the + set up of different trials (with different motifs/categories) in the experiment without the necessity to init to + classifiers: default None + :return: response, a tuple of result (bool) and response body + Response body is used for plotting and outputting results to trials dataframes + """ + self.fill_time_window(skeleton) + # self._time_window.append(temp_feature) + # self._time_window = temp_feature + f_extract_output = None + """Checks if necessary time window was collected and passes it to classifier""" + if len(self._time_window) == self._time_window_len: + start_time = time.time() + f_extract_output = self.feat_extractor.extract_features(self._time_window) + end_time = time.time() + if self._debug: + print( + "Feature extraction time: {:.2f} msec".format( + (end_time - start_time) * 1000 + ) + ) + if f_extract_output is not None: + self._last_result, _, _ = self._classifier.classify(f_extract_output) + else: + self._last_result = [0] + if trigger is not None: + self._trigger = trigger + # choosing a point to draw near the skeleton + self._center = skeleton[list(skeleton.keys())[0]] + # self._center = (50,50) + result = False + # text = 'Current probability: {:.2f}'.format(self._last_prob) + text = "Current Class: {}".format(self._last_result) + + if self._last_result[0] == self._trigger: + result = True + text = "Motif matched: {}".format(self._last_result) + + color = (0, 255, 0) if result else (0, 0, 255) + response_body = { + "plot": {"text": dict(text=text, org=self._center, color=color)} + } + response = (result, response_body) + return response + + def get_trigger_threshold(self): + return self._trigger + + def get_last_prob(self): + return self._last_prob + + def get_time_window_len(self): + return self._time_window_len + + +class BsoidClassBehaviorPoolTrigger_old: + """ + THIS VERSION HAS A SEPERATE FEATURE EXTRACTION WHICH CAN AFFECT OVERALL PERFORMANCE + + Trigger to check if animal's behavior is classified as specific motif with BSOID trained classifier. + """ + + def __init__(self, target_class: int, class_process_pool, debug: bool = False): + """ + Initialising trigger with following parameters: + :param int target_class: target classification category that should be used as trigger. Must match "Group" number of cluster in BSOID. + If you plan to use the classifier for multiple trial triggers in the same experiment with different thresholds. We recommend setting up the + target_class during check_skeleton + :param class_process_pool: list of dictionaries with keys process: mp.Process, input: mp.queue, output: mp.queue; + used for lossless frame-by-frame classification + + """ + self._trigger = target_class + self._process_pool = class_process_pool + self._last_result = [0] + self._feature_id = 0 + self._center = None + self._debug = debug + self._skeleton = None + self._time_window_len = TIME_WINDOW + self.feat_extractor = BsoidFeatureExtractor(self._time_window_len) + self._time_window = deque(maxlen=self._time_window_len) + + def fill_time_window(self, skeleton): + from utils.poser import transform_2pose + + pose = transform_2pose(skeleton) + self._time_window.appendleft(pose) + + def check_skeleton(self, skeleton, target_class: int = None): + """ + Checking skeleton for trigger, will pass skeleton window to classifier if window length is reached and + collect skeletons otherwise + :param skeleton: a skeleton dictionary, returned by calculate_skeletons() from poser file + :param target_class: optional, overwrites self._trigger with target probability. this is supposed to enable the + set up of different trials (with different motifs/categories) in the experiment without the necessity to init to + classifiers: default None + :return: response, a tuple of result (bool) and response body + Response body is used for plotting and outputting results to trials dataframes + """ + self.fill_time_window(skeleton) + # self._time_window.append(temp_feature) + # self._time_window = temp_feature + f_extract_output = None + """Checks if necessary time window was collected and passes it to classifier""" + if len(self._time_window) == self._time_window_len: + start_time = time.time() + f_extract_output = self.feat_extractor.extract_features(self._time_window) + end_time = time.time() + if self._debug: + print( + "Feature extraction time: {:.2f} msec".format( + (end_time - start_time) * 1000 + ) + ) + if f_extract_output is not None: + self._feature_id += 1 + self._process_pool.pass_features( + (f_extract_output, self._feature_id), debug=self._debug + ) + # check if a process from the pool is done with the result + clf_result, feature_id = self._process_pool.get_result(debug=self._debug) + if clf_result is not None: + self._last_result = clf_result[0] + if target_class is not None: + self._trigger = target_class + # choosing a point to draw near the skeleton + self._center = skeleton[list(skeleton.keys())[0]] + # self._center = (50,50) + result = False + # text = 'Current probability: {:.2f}'.format(self._last_prob) + text = "Current Class: {}".format(self._last_result) + + if self._last_result[0] == self._trigger: + result = True + text = "Motif matched: {}".format(self._last_result) + + color = (0, 255, 0) if result else (0, 0, 255) + response_body = { + "plot": {"text": dict(text=text, org=self._center, color=color)} + } + response = (result, response_body) + return response + + def get_trigger_threshold(self): + return self._trigger + + def get_last_prob(self): + return self._last_prob + + def get_time_window_len(self): + return self._time_window_len diff --git a/requirements.txt b/requirements.txt index 9ade97f..fb8c932 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,4 @@ + PySide2==5.15.2 gpiozero==1.5.1 pigpio==1.78 @@ -7,6 +8,9 @@ click==7.1.2 opencv-python==3.4.5.20 opencv-contrib-python==4.4.0.46 numpy>=1.14.5 -pandas==1.2.2 +pandas==1.1.4 +scikit-learn== 0.24.1 scikit-image==0.17.2 scipy==1.4.1 +pure-predict==0.0.4 +numba==0.51.1 diff --git a/settings.ini b/settings.ini index 0532ce3..b77fd75 100644 --- a/settings.ini +++ b/settings.ini @@ -1,20 +1,20 @@ [Streaming] -RESOLUTION = 1920, 1080 +RESOLUTION = 960, 540 FRAMERATE = 30 OUTPUT_DIRECTORY = /Output #if you have connected multiple cameras (USB), you will need to select the number OpenCV has given them. #Default is "0", which takes the first available camera. CAMERA_SOURCE = 0 -#you can use camera, ipwebcam or video to select your input source -STREAMING_SOURCE = video +#you can use "camera", "ipwebcam" or "video" to select your input source +STREAMING_SOURCE = camera [Pose Estimation] #possible origins are: SLEAP, DLC, DLC-LIVE,MADLC, DEEPPOSEKIT -MODEL_ORIGIN = ORIGIN +MODEL_ORIGIN = MODEL_ORIGIN #takes path to model or models (in case of SLEAP topdown, bottom up) in style "string" or "string , string", without "" # E.g.: MODEL_PATH = D:\SLEAP\models\baseline_model.centroids , D:\SLEAP\models\baseline_model.topdown MODEL_PATH = PATH_TO_MODEL -MODEL_NAME = MODEL_NAME +MODEL_NAME = NAME_OF_MODEL ; only used in DLC-LIVE and DeepPoseKit for now; if left empty or to short, auto-naming will be enabled in style bp1, bp2 ... ALL_BODYPARTS = bp1, bp2, bp3, bp4 @@ -26,6 +26,19 @@ EXP_NAME = ExampleExperiment #if you want the experiment to be recorded as a raw video set this to "True". RECORD_EXP = True +[Classification] +PATH_TO_CLASSIFIER = PATH_TO_CLASSIFIER +#time window used for feature extraction (currently only works with 15) +TIME_WINDOW = 15 +#number of parallel classifiers to run, this is dependent on your performance time. You need at least 1 more classifier then your average classification time. +POOL_SIZE = 1 +#threshold to accept a classification probability as positive detection (SIMBA + ) +THRESHOLD = 0.9 +# class/category of identified behavior to use as trigger (only used for B-SOID) +TRIGGER = NUMBER_OF_CLUSTER +#feature extraction currently works with millimeter not px, so be sure to enter the factor (as in simba). +PIXPERMM = 1 + [Video] #Full path to video that you want to use as input. Needs "STREAMING_SOURCE" set to "video"! VIDEO_SOURCE = PATH_TO_VIDEO diff --git a/utils/configloader.py b/utils/configloader.py index 11ca3c2..03ec8b6 100644 --- a/utils/configloader.py +++ b/utils/configloader.py @@ -39,7 +39,8 @@ def get_script_path(): MODEL_PATH = model_path_string[0] if len(model_path_string) <= 1 else model_path_string MODEL_NAME = dsc_config["Pose Estimation"].get("MODEL_NAME") ALL_BODYPARTS = tuple( - part for part in dsc_config["Pose Estimation"].get("ALL_BODYPARTS").split(",") + part.strip() + for part in dsc_config["Pose Estimation"].get("ALL_BODYPARTS").split(",") ) # Streaming items @@ -73,6 +74,15 @@ def get_script_path(): START_TIME = time.time() +# Classification +PATH_TO_CLASSIFIER = dsc_config["Classification"].get("PATH_TO_CLASSIFIER") +PIXPERMM = dsc_config["Classification"].getfloat("PIXPERMM") +THRESHOLD = dsc_config["Classification"].getfloat("THRESHOLD") +TRIGGER = dsc_config["Classification"].getfloat("TRIGGER") +POOL_SIZE = dsc_config["Classification"].getint("POOL_SIZE") +TIME_WINDOW = dsc_config["Classification"].getint("TIME_WINDOW") + + """advanced settings""" STREAMS = [ str(part).strip() for part in adv_dsc_config["Streaming"].get("STREAMS").split(",") diff --git a/utils/poser.py b/utils/poser.py index bf1b86d..21fe9d7 100644 --- a/utils/poser.py +++ b/utils/poser.py @@ -176,6 +176,28 @@ def find_local_peaks_new( return all_peaks +def filter_pose_by_likelihood(pose, threshold: float = 0.1): + """ + !!!FOR NOW THIS FUNCTION SETS THEM TO (0,0) due to missing float dtype in coordinate handling. Will be updated later on. + filters pose estimation by likelihood threshold. Estimates below threshold are set to NaN and handled downstream + of this function in calculate skeletons. + :param pose: pose estimation (e.g., from DLC) + :param threshold: likelihood threshold to filter by + :return filtered_pose: pose estimation filtered by likelihood (may contain NaN) + """ + # TODO: Update to NaN + + filtered_pose = pose.copy() + + for num, bp in enumerate(filtered_pose): + if bp[2] < threshold: + # set new threshold to "2" (number outside of normal range to signify filter + # TODO: This should be np.NaN not 0,0 + filtered_pose[num] = np.array([0, 0, 2]) + + return filtered_pose + + def calculate_dlstream_skeletons(peaks: dict, animals_number: int) -> list: """ Creating skeletons from given peaks @@ -272,8 +294,8 @@ def get_ma_pose(image, config, session, inputs, outputs): inputs, outputs, outall=True, - nms_radius=config["nmsradius"], - det_min_score=config["minconfidence"], + nms_radius=5.0, + det_min_score=0.1, c_engine=False, ) @@ -459,6 +481,38 @@ def handle_missing_bp(animal_skeletons: list): return animal_skeletons +def arrange_flatskeleton(skeleton, n_animals, n_bp_animal, switch_dict): + """changes sequence of bodypart sets (skeletons) in multi animal tracking with flat skeleton output (multiple animals in single skeleton) by switching position of pairs. + E.g. in pose estimation with different fur colors. Note: When switching muliple animals the new position of the previous switches will be used. + :param skeleton: flat skeleton of pose estimation in style {bp1: (x,y), bp2: (x2,y2) ...} + :param n_animals: number of animals in total represented by skeleton + :param n_bp_animal: number of bodyparts per animal in skeleton + :param switch_dict: dictionary containing position of bodypart set (animal) in flat skeleton as key and bp set to exchange with as value. + e.g. switch_dict = dict(1 = 2, 3 = 4) + :return: skeleton with new order + """ + flat_pose = transform_2pose(skeleton) + ra_dict = {} + # slicing the animals out + for num_animal in range(n_animals): + ra_dict[num_animal] = flat_pose[ + num_animal * n_bp_animal : num_animal * n_bp_animal + n_bp_animal + ] + # switching positions + for orig_pos, switch_pos in switch_dict.items(): + # extract old + orig = ra_dict[orig_pos] + switch = ra_dict[switch_pos] + # set to new position + ra_dict[orig_pos] = switch + ra_dict[switch_pos] = orig + # extracting pose + arranged_pose = np.array([*ra_dict.values()]).reshape(flat_pose.shape) + # transforming it to skeleton + flat_skeleton = transform_2skeleton(arranged_pose) + return flat_skeleton + + def calculate_skeletons_dlc_live(pose) -> list: """ Creating skeletons from given pose