Python files

00_clip_video.py

@@ -0,0 +1,251 @@
+
+# This code clips the raw videos based on claps detected by the user during the execution and stores them in a new file
+
+import subprocess
+import glob
+import sounddevice as sd
+import argparse
+import os
+
+
+from scipy.signal import butter, filtfilt
+from scipy.io import wavfile
+from pylab import *
+
+def butter_lowpass(cutoff_low, cutoff_high, fs, order=5):
+    nyq = 0.5 * fs
+    cutoff_low = cutoff_low / nyq
+    cutoff_high = cutoff_high / nyq
+    b, a = butter(order, [cutoff_low, cutoff_high], btype='bandpass',
+                  analog=False)
+    return b, a
+
+
+def butter_lowpass_filter(data, cutoff_low, cutoff_high, fs, order=5):
+    b, a = butter_lowpass(cutoff_low, cutoff_high, fs, order=order)
+    y = filtfilt(b, a, data)
+    return y
+
+
+def find_start_frame_by_clap(path_to_audio_file, startt=0, stopt=500,
+                             where=None):
+    """
+    Detect clap in audio file. Sound is played and detection must be confirmed
+    by pressing [d] or discarded by pressing [e]. [c] to restart, [a] to
+    listen again,
+
+    :param
+    path_to_audio_file: path to .wav file
+    startt: default 0, time in seconds to start search
+    stopt: default 100, time in seconds to end search. By default only within
+    the first minute and 30 seconds the clap is searched.
+    where: valuesallowed are 'start' or 'end' if provided first half or last
+    half resp. is set to zero
+
+    :return:
+    second_of_clap: Second of .wav file in which clap takes place
+    """
+    audioFreq, audio = wavfile.read(path_to_audio_file)
+    audio.setflags(write=1)
+
+    # set part of sequence to zero if searched for start / end clap
+    if where != None:
+        if where == 'end':
+            audio[0:int(len(audio) / 2)] = 0
+        elif where == 'start':
+            audio[int(len(audio) / 2):len(audio)] = 0
+        else:
+            print('where must be \'start\' oder \'end\'')
+
+    cutoff_low = 1000  # desired cutoff frequency of the filter, Hz
+    cutoff_high = 1400
+    fs = 3000
+    order = 5
+
+    foundStartClap = False
+    k = 0
+    while not foundStartClap:
+        k += 10
+        s = int(audioFreq * startt)
+        e = int(audioFreq * (stopt + k))
+        # Filter .wav file within given range of time and find argmax amplitude
+        y = butter_lowpass_filter(audio[s:e, 0], cutoff_low,
+                                  cutoff_high, fs, order)
+        maxamp = np.argmax(y)
+        # -/+ audioFreq is necessary to actually hear something (1 second)
+        # when audio is played.
+        l = int(maxamp - 2 * audioFreq)
+        u = int(maxamp + 2 * audioFreq)
+        print(l / audioFreq, u / audioFreq)
+        sd.play(audio[l:u], audioFreq, blocking=True)
+        print('Clap playing second: ', np.ceil(maxamp / audioFreq))
+        # Capture keyboard events
+        while True:
+            print("Press [d] if clap was detected; [c] if sound was no clap; "
+                  "[a] to hear again; [e] to pass; press [n] if clap not "
+                  "found at all, second of clap will be set to one")
+            x = input()
+            if x == 'd':
+                second_of_clap = maxamp / audioFreq
+                print('Clap accepted; takes place at second ', second_of_clap)
+                foundStartClap = True
+                break
+            elif x == 'c':
+                print('Restart clap search')
+                audio[l:u] = 0
+                break
+            elif x == 'a':
+                sd.play(audio[l:u], audioFreq, blocking=True)
+            elif x == 'e':
+                return 0
+            elif x == 'n':
+                second_of_clap = 1
+                foundStartClap = True
+                print('Clap not found, video will start at ', second_of_clap)
+                break
+            else:
+                print('Key must be [d],[c], [a] or [e]')
+
+    return second_of_clap
+
+def clip_video_from_clap(video_path, output_path, start_time, end_time):
+
+    command_clip = (
+        f"ffmpeg -i {video_path} -ss {start_time} -to {end_time} -c copy {output_path}"
+    )
+    subprocess.call(command_clip, shell=True)
+
+
+if __name__ == "__main__":
+
+    """ 
+    This script processes video all in once. The default directory for videos
+    being processed is data/Videos. Frames are automatically extracted after
+    clap detection finished. The following steps are performed:
+        1) Create folder based on video name
+        2) Extract audio file
+        3) Clap detection (starting and ending clap) using low pass filter in 
+        video for synchronization of the three cameras.
+        Make sure sound is turned on. Clap detection results are written to 
+        offset.csv. Headphones necessary to confirm detection.
+        4) Extract video frames using ffmpeg
+
+        Args:
+            inputpath: Path to folder with videos.
+
+        Returns:
+            Folder containing audio file and video frames. Writes number of 
+            frame where clap takes place and second of clap to offset.csv. 
+    """
+
+    parser = argparse.ArgumentParser(description='Processes all Videos in '
+                                                 '--inputpath at once')
+    parser.add_argument('--inputpath',
+                        default='/home/valapil/Project/ForkCausal_Adithya/raw_video/pair3',
+                        help='Path to folder with videos to process')
+    parser.add_argument('--outputpath',
+                        default='/home/valapil/Project/ForkCausal_Adithya/clipped_vid',
+                        help='Path to folder with pair-folders')
+    parser.add_argument('--ExtractFrames', action='store_false',
+                        help='If true frames are extracted from video')
+    args = parser.parse_args()
+    print(args, args.ExtractFrames)
+
+    # Comment in if side view perspective is processed as well.
+    # videofiles = glob.glob(args.inputpath + '/*.MTS')
+    # videofiles.extend(glob.glob(args.inputpath + '/*.MOV'))
+    # videofiles = sorted(videofiles)
+
+    videofiles = sorted(glob.glob(args.inputpath + '/*.MTS'))
+
+    print('Videos to process: ', args.inputpath, videofiles)
+
+    f = open(os.path.join(args.outputpath, 'offset.csv'), 'a')
+
+    # Extract clap and audio
+    for video_dir in videofiles:
+
+        accept_Video = True
+        print('Processing video: ', video_dir)
+
+        video = os.path.basename(video_dir)
+        video_file = os.path.basename(video).split('.')[0]
+
+        video_file = video_file.replace('r', 'o')
+        print('PROCESS ', video_dir)
+
+        # Create directory for pair
+        pair = 'pair_' + str(video_file[0:3])
+        create_folder_pair = os.path.join(args.outputpath, pair)
+        command_folder = "mkdir -p " + create_folder_pair
+        subprocess.call(command_folder, shell=True)
+
+        # Create directory for images in outputpath
+        create_folder_img = os.path.join(args.outputpath, pair,
+                                         video_file, 'images')
+        command_folder_img = "mkdir -p " + create_folder_img
+        subprocess.call(command_folder_img, shell=True)
+
+        # Create directory for audio file in outputpath
+        create_folder_audio = os.path.join(args.outputpath, pair,
+                                           video_file, 'audio')
+        command_folder_audio = "mkdir -p " + create_folder_audio
+        subprocess.call(command_folder_audio, shell=True)
+
+        # Extract audio from video
+        aud = os.path.join(create_folder_audio,
+                           video_file + "_audio.wav")
+        command_audio = "ffmpeg -i " + video_dir + \
+                        " -ab 160k -ac 2 -ar 44100 -vn " + aud
+        subprocess.call(command_audio, shell=True)
+
+        # Detect strating and ending claps in video
+        for i in ['start', 'end']:
+            second_of_clap = find_start_frame_by_clap(path_to_audio_file
+                                                      =aud, where=i)
+            # Video camera from media department format
+            if os.path.basename(video).split('.')[1] == 'MTS':
+                frame_rate = 25
+                frame_of_clap = int(second_of_clap * frame_rate)
+            # ipad video format
+            elif os.path.basename(video).split('.')[1] == 'MOV':
+                # original frame rate is of .MOV is 29.984664, but side
+                # view should be synchronized with frontal perspective
+                frame_rate = 25
+                frame_of_clap = int(second_of_clap * frame_rate)
+
+            print(video_file, i, second_of_clap, frame_of_clap)
+
+            if i == 'start':
+                start_clip = second_of_clap
+            elif i == 'end':
+                if second_of_clap == start_clip:
+                    accept_Video = False
+                else:
+                    end_clip = second_of_clap
+
+        if accept_Video == True:
+            clip_output_path = os.path.join(create_folder_pair, video_file, f"{video}_clipped.mp4")
+            clip_video_from_clap(video_path=video_dir, output_path=clip_output_path, start_time=start_clip, end_time= end_clip)
+            # f.write(f'{video_file} {i} {second_of_clap} {frame_of_clap}\n')
+
+    f.close()
+
+    # # Extract video frames to pair/condition/images
+    # if args.ExtractFrames:
+    #     print('Start extracting frames ... ')
+    #     for video_dir in videofiles:
+    #         video = os.path.basename(video_dir)
+    #         video_file = os.path.basename(video).split('.')[0]
+    #         # you can exclude videos (e.g. sideview)
+    #         # processedVideos = ['I am already processed']
+    #         # if np.isin(video_file, processedVideos, invert=True):
+    #         pair = 'pair_' + str(video_file[0:3])
+    #         create_folder_pair = os.path.join(args.outputpath, pair)
+    #         create_folder = os.path.join(create_folder_pair, video_file,
+    #                                      'images')
+    #         vid = create_folder + '/' + video_file + '_%05d.png'
+    #
+    #         # Call from shell
+    #         command_frames = "ffmpeg -i " + video_dir + " -r 25 " + vid
+    #         subprocess.call(command_frames, shell=True)

02_process_offset.py

@@ -0,0 +1,68 @@
+# use the openface 2 to get the action units before this
+# go to- valapil@hemera4:~/Programs/OpenFace$
+# type- ./bin/FeatureExtraction -f "/home/valapil/Project/ForkCausal_Adithya/clipped_vid/pair_021/021S3o/021S3o.MTS_clipped.mp4"
+
+# also get the offset values
+
+import numpy as np
+import os
+
+
+AU = ["Person", "Case", "AU01", "AU02", "AU04", "AU05", "AU06", "AU07", "AU09","AU10","AU12","AU14","AU15","AU17","AU20","AU23",
+      "AU25","AU26", "AU45"]
+off = ["Case", "Start", "End", "Person"]
+
+# folder where output csv files of OpenFace are stored
+folder = '/home/valapil/Project/ForkCausal_Adithya/csv_files'
+csv_files = [file for file in os.listdir(folder) if file.endswith('.csv')]
+
+def loadOffset():
+    """
+    get offset values for video frames
+    :param basepath: path to where offset.csv is at
+    :return: dictionary with {identifier: offset}('002S1g': 7}
+    """
+
+    # opening the file where offset values are present
+    with open(os.path.join('/home/valapil/Project/ForkCausal_Adithya/offset.csv'), 'r') as f:
+        f.readline() # read header
+        offset = {}
+        for row in f:
+            id, StartEnd, second, frame = row.split(' ')[0:4]
+            if id not in offset:
+                offset[id] = {}
+            offset[id][StartEnd] = int(np.ceil(float(frame)))
+        return offset
+
+
+for file in csv_files:
+    base_name, extension = os.path.splitext(file)
+    new_name = base_name.split('.')[0]
+    f1 = new_name[:-2]
+    f2 = new_name[-1:]
+    file_path = os.path.join(folder, file)
+
+    offset = loadOffset()
+    if new_name[-5] != '0':
+        offset_file = np.hstack((new_name[:-3] + new_name[-1], offset[new_name]['start'], offset[new_name]['end'], new_name))
+    else:
+        offset_file = np.hstack((new_name[:-3] + new_name[-1], offset['0' + new_name]['start'], offset['0' + new_name]['end'], new_name))
+
+    off = np.row_stack((off, offset_file))
+
+# file where offset values are stored as npy format
+np.save(f'/home/valapil/Project/ForkCausal_Adithya/offset_values.npy', off)
+
+
+
+
+
+
+
+
+
+
+
+
+
+

03_process_data.py

@@ -0,0 +1,42 @@
+import numpy as np
+import os
+
+
+AU = ["Person", "Case", "Confidence", "AU01", "AU02", "AU04", "AU05", "AU06", "AU07", "AU09","AU10","AU12","AU14","AU15","AU17","AU20","AU23",
+  "AU25","AU26", "AU45"]
+
+# folder where output csv files of OpenFace are stored
+folder = '/home/valapil/Project/ForkCausal_Adithya/csv_files'
+csv_files = [file for file in os.listdir(folder) if file.endswith('.csv')]
+
+# load the offset values
+offset = np.load(f'/home/valapil/Project/ForkCausal_Adithya/offset_values.npy')
+length = {}
+same = np.unique(offset[1:, 0])
+for i in same:
+    same_rows = offset[offset[:, 0] == i]
+    length[i] = max(abs(int(same_rows[0][2]) - int(same_rows[0][1])), abs(int(same_rows[1][2]) - int(same_rows[1][1])))
+
+for file in csv_files:
+    base_name, extension = os.path.splitext(file)
+    new_name = base_name.split('.')[0]
+    f1 = new_name[:-2]
+    f2 = new_name[-1:]
+    file_path = os.path.join(folder, file)
+    len = length[new_name[:-3] + new_name[-1]]
+    start = np.where(offset[1:, 3]  == new_name)[0][0]
+    start = int(offset[start+1, 1])
+    data_all = np.genfromtxt(file_path, delimiter=',', skip_header=1)
+    data = data_all[start:len-200, 679:696]
+
+    selected_columns = np.hstack((np.full((data.shape[0], 2), [f1, f2]), data_all[start:len-200, [3]], data))
+
+    AU = np.row_stack((AU, selected_columns))
+
+    #  save each offset data here
+    np.save(f'/home/valapil/Project/ForkCausal_Adithya/data/{new_name}.npy', data)
+
+#  save all offset data with the particular case and confidence value here
+np.save(f'/home/valapil/Project/ForkCausal_Adithya/processed_data.npy', AU)
+
+

04_get_expression.py

@@ -0,0 +1,95 @@
+import numpy as np
+from numpy import save
+import os
+
+
+# AU = ["AU01", "AU02", "AU04", "AU05", "AU06", "AU07", "AU09","AU10","AU12","AU14","AU15","AU17","AU20","AU23",
+#       "AU25","AU26","AU45","AU01_c","AU02_c","AU04_c","AU05_c","AU06_c","AU07_c","AU09_c","AU10_c","AU12_c","AU14_c","AU15_c","AU17_c",
+#       "AU20_c","AU23_c","AU25_c","AU26_c","AU28_c","AU45_c"]
+
+reference_au = ["File", "AU01", "AU02", "AU04", "AU05", "AU06", "AU07", "AU09","AU10","AU12","AU14","AU15","AU17","AU20","AU23",
+      "AU25","AU26", "AU45"]
+
+
+# finding reference values
+ref_au = {}
+ref = []
+
+#  get all offset data
+data_all = np.load(f'/home/valapil/Project/ForkCausal_Adithya/processed_data.npy')
+# data = np.concatenate(data_all[:3], data_all[4:])
+data = np.delete(data_all, 2,1)
+activation = []
+
+
+unique = np.unique(data[1:, 0])
+
+for i in unique:
+    same_rows = data[data[:, 0] == i]
+    mean = np.mean(same_rows[:, 2:19].astype(np.float64), axis=0)
+    std = np.std(same_rows[:, 2:19].astype(np.float64), axis=0)
+    ref_au = mean + (0.5 * std)
+    # ref_au.append(mean + (0.5 * std))
+    ref = np.hstack((i, ref_au))
+    reference_au = np.row_stack((reference_au, ref))
+
+# folder where each offset data is stored
+folder = '/home/valapil/Project/ForkCausal_Adithya/data'
+csv_files = [file for file in os.listdir(folder) if file.endswith('.npy')]
+
+
+
+for file in csv_files:
+    EXP = ["Happiness upper", "Happiness lower", "Surprise upper", "Surprise lower", "Disgust lower", "Fear upper",
+            "Fear lower", "Sadness upper", "Sadness lower", "Anger upper", "Anger lower"]
+
+    # Get data
+    file_path = os.path.join(folder,file)
+    # data = np.genfromtxt(file_path, delimiter=',', skip_header=1)
+    data = np.load(file_path)
+
+    # find the particular person
+    base_name, extension = os.path.splitext(file)
+    new_name = base_name.split('.')[0]
+    f1 = new_name[:-2]
+    selected_columns = data
+
+    index = np.where(reference_au[1:, 0] == f1)[0]
+
+    # AU activation
+    change = np.subtract(selected_columns, reference_au[index+1, 1:].astype(np.float16))
+    activation = np.where(change >= 0, 1, 0)
+
+    expression_labels = np.zeros((activation.shape[0], 11), dtype=int)
+
+    # Define expressions based on specific AUs
+    expression_labels[:, 0] = np.any(activation[:, [4]], axis=1)  # AU06 is 4th au in order. [4] to get results with []
+    expression_labels[:, 1] = np.all(activation[:, [8, 14]], axis=1)  # AU12 and AU25
+    expression_labels[:, 2] = np.all(activation[:, [0, 1, 3]], axis=1)  # AU01, AU02, AU05
+    expression_labels[:, 3] = np.any(activation[:, [15]], axis=1)  # AU26
+    expression_labels[:, 4] = np.all(activation[:, [6, 7, 14]], axis=1)  # AU09, AU10, AU25
+    expression_labels[:, 5] = np.all(activation[:, [0, 1, 2, 3]], axis=1)  # AU01, AU02, AU04, AU05
+    expression_labels[:, 6] = np.all(activation[:, [12, 14]], axis=1)  # AU20, AU25
+    expression_labels[:, 7] = np.all(activation[:, [0, 2]], axis=1)  # AU01, AU04
+    expression_labels[:, 8] = np.all(activation[:, [10, 11]], axis=1)  # AU15, AU17
+    expression_labels[:, 9] = np.all(activation[:, [2, 3, 5]], axis=1)  # AU04, AU05, AU07
+    expression_labels[:, 10] = np.all(activation[:, [11, 13]], axis=1)  # AU17, AU23
+
+    EXP = np.row_stack((EXP, expression_labels))
+
+    # Save the result to a new NumPy file
+    # Each file has expressions for each frame based on activation of AU in that frame
+    np.save(f'/home/valapil/Project/ForkCausal_Adithya/test/{new_name}.npy', EXP)
+
+
+
+
+
+
+
+
+
+
+
+
+

05_wilcoxon.py

@@ -0,0 +1,128 @@
+# normalisation and results for wilcoxon test
+
+import numpy as np
+import pandas as pd
+from numpy import load
+from numpy import save
+import os
+from scipy.stats import wilcoxon
+
+# folder where expression activation per frame of each pair is stored
+expressions_folder = '/home/valapil/Project/ForkCausal_Adithya/test/'
+
+# EXP = ["CASE", "Happiness", "Surprise", "Disgust", "Fear", "Sadness", "Anger"]
+
+EXP = ["pair", "role", "case", "Happiness upper", "Happiness lower", "Surprise upper", "Surprise lower", "Disgust lower", "Fear upper",
+        "Fear lower", "Sadness upper", "Sadness lower", "Anger upper", "Anger lower"]
+
+
+np_files = [file for file in os.listdir(expressions_folder)]
+
+expr_sum = None
+
+frames = {}
+
+# get frames
+for file in np_files:
+    file_path = os.path.join(expressions_folder, file)
+    data = np.load(file_path)
+    base_name, extension = os.path.splitext(file)
+    new_name = base_name.split('.')[0]
+    frames[new_name] = data.shape[0]
+
+# normalisation 1
+for file in np_files:
+    file_path = os.path.join(expressions_folder, file)
+    data = np.load(file_path)
+    data = data[1:, :].astype(np.float64)
+
+    base_name, extension = os.path.splitext(file)
+    new_name = base_name.split('.')[0]
+    frame = min(frames[new_name[:-3] + 'S' + new_name[-2:]], frames[new_name[:-3] + 'E' + new_name[-2:]])
+    expr_sum = np.sum(data, axis=0)
+    norm_frames = np.divide(100 * expr_sum, frame)
+
+    total = np.hstack((new_name[:-3], new_name[-3], new_name[-1], norm_frames))
+    EXP = np.row_stack((EXP, total))
+
+# normalisation 2
+df = pd.DataFrame(EXP[1:, :])
+df.columns = ["pair", "role", "case", "Happiness upper", "Happiness lower", "Surprise upper", "Surprise lower",
+              "Disgust lower", "Fear upper", "Fear lower", "Sadness upper", "Sadness lower", "Anger upper", "Anger lower"]
+
+num_col = ["Happiness upper", "Happiness lower", "Surprise upper", "Surprise lower",
+              "Disgust lower", "Fear upper", "Fear lower", "Sadness upper", "Sadness lower", "Anger upper", "Anger lower"]
+df[num_col] = df[num_col].apply(pd.to_numeric, errors='coerce')
+
+# aumax1 = df.groupby(['pair', 'role']).sum()
+aumax2 = df.groupby(['pair', 'role'])[num_col].transform('sum')
+
+
+result = EXP[1:,:].copy()
+result[:, 3:] = df[num_col] / aumax2
+
+
+# need to group by case
+s1 = np.sum(EXP[1:, 3:].astype(np.float64), axis=0)
+s2 = np.multiply((np.divide(s1, sum(s1))), 100)
+
+output = ["Happiness upper", "Happiness lower", "Surprise upper", "Surprise lower", "Disgust lower", "Fear upper",
+        "Fear lower", "Sadness upper", "Sadness lower", "Anger upper", "Anger lower"]
+# output = np.row_stack((output , s2))
+
+print(output, s2)
+
+benjhoch = []
+testresult = {}
+intersect = []
+
+# Wilcoxon test
+# for expr in range(0, output.shape[1]):
+for expr, value in enumerate(output):
+    for cond in [['g', 'a'], ['a', 'o'], ['g', 'o']]:
+        index = np.where(result[:, 2] == cond[0])[0]
+        x = result[index, 3+expr].astype(np.float64)
+        index = np.where(result[:, 2] == cond[1])[0]
+        y = result[index, 3+expr].astype(np.float64)
+        res = wilcoxon(x=x, y=y, zero_method='pratt')
+        benjhoch += [[value, res[1], cond[0], cond[1]]]
+        testresult[value] = [cond, res]
+
+
+Q = 0.3
+benjhoch = np.array(benjhoch,dtype=object)
+benjhoch = np.column_stack((benjhoch[benjhoch[:, 1].argsort()],
+                            np.arange(1, len(benjhoch) + 1)))
+print('Benjamini Hochberg Table: \n', benjhoch)
+m = len(output) * 3
+
+# Find maximum rank where actual p-value is smaller
+# than hochberg-benjamini p-value
+maxi = 0
+for i, val in enumerate(benjhoch):
+    if float(val[1]) < (int(val[4]) / m) * Q:
+        maxi = i
+
+# Check if test results is significant when using a Benjamini Hochberg
+# p-value corretion
+print(maxi + 1, 'significantly changing expression(s): \n')
+sigexpr = []
+for val in benjhoch[:(maxi + 1)]:
+    au = val[0]
+    sigexpr += [au]
+    index = np.where(result[:, 2] == [val[2]])[0]
+    x = result[index, 3 + output.index(au)].astype(np.float64)
+    index = np.where(result[:, 2] == [val[3]])[0]
+    y = result[index, 3 + output.index(au)].astype(np.float64)
+    res = wilcoxon(x=x, y=y, zero_method='pratt')
+
+    print('Expression: ', val[0], 'significantly changes between '
+                                  'conditions', val[2], 'and', val[3])
+
+    print('Average expression activation \n',
+          val[2], ':', np.nanmean(x), val[3], ':', np.nanmean(y), '\n')
+
+
+
+
+# df.to_csv('/home/valapil/Project/ForkCausal_Adithya/analysis.csv', index=False)