<commit_before># -*- coding: utf-8 -*-
"""
Created on Thu Apr  2 16:33:34 2015

@author: ajaver
"""


import numpy as np
import tables
from math import sqrt
import cv2
from skimage.filters import threshold_otsu
import os
import sys

from sklearn.utils.linear_assignment_ import linear_assignment  # hungarian algorithm
from scipy.spatial.distance import cdist

from MWTracker.helperFunctions.timeCounterStr import timeCounterStr
from MWTracker.compressVideos.extractMetaData import readAndSaveTimestamp

table_filters = tables.Filters(complevel=5, complib='zlib', shuffle=True)


class plate_worms(tables.IsDescription):
    # class for the pytables
    worm_index_blob = tables.Int32Col(pos=0)
    worm_index_joined = tables.Int32Col(pos=1)

    frame_number = tables.Int32Col(pos=2)
    coord_x = tables.Float32Col(pos=3)
    coord_y = tables.Float32Col(pos=4)
    area = tables.Float32Col(pos=5)
    perimeter = tables.Float32Col(pos=6)
    box_length = tables.Float32Col(pos=7)
    box_width = tables.Float32Col(pos=8)
    quirkiness = tables.Float32Col(pos=9)
    compactness = tables.Float32Col(pos=10)
    box_orientation = tables.Float32Col(pos=11)
    solidity = tables.Float32Col(pos=12)
    intensity_mean = tables.Float32Col(pos=13)
    intensity_std = tables.Float32Col(pos=14)

    threshold = tables.Int32Col(pos=15)
    bounding_box_xmin = tables.Int32Col(pos=16)
    bounding_box_xmax = tables.Int32Col(pos=17)
    bounding_box_ymin = tables.Int32Col(pos=18)
    bounding_box_ymax = tables.Int32Col(pos=19)

    # deprecated, probably it would be good to remove it in the future
    segworm_id = tables.Int32Col(pos=20)

    hu0 = tables.Float32Col(pos=21)
    hu1 = tables.Float32Col(pos=22)
    hu2 = tables.Float32Col(pos=23)
    hu3 = tables.Float32Col(pos=24)
    hu4 = tables.Float32Col(pos=25)
    hu5 = tables.Float32Col(pos=26)
    hu6 = tables.Float32Col(pos=27)


def _getWormThreshold(pix_valid):
    
    # calculate otsu_threshold as lower limit. Otsu understimates the threshold.
    try:
        otsu_thresh = threshold_otsu(pix_valid)
    except:
        return np.nan

    # calculate the histogram
    pix_hist = np.bincount(pix_valid)

    # the higher limit is the most frequent value in the distribution
    # (background)
    largest_peak = np.argmax(pix_hist)
    if otsu_thresh < largest_peak and otsu_thresh + 2 < len(pix_hist) - 1:
        # smooth the histogram to find a better threshold
        pix_hist = np.convolve(pix_hist, np.ones(3), 'same')
        cumhist = np.cumsum(pix_hist)

        xx = np.arange(otsu_thresh, cumhist.size)
        try:
            # the threshold is calculated as the pixel level where there would be
            # larger increase in the object area.
            hist_ratio = pix_hist[xx] / cumhist[xx]
            thresh = np.where(
                (hist_ratio[3:] - hist_ratio[:-3]) > 0)[0][0] + otsu_thresh
        except IndexError:
            thresh = np.argmin(
                pix_hist[
                    otsu_thresh:largest_peak]) + otsu_thresh
    else:
        # if otsu is larger than the maximum peak keep otsu threshold
        thresh = otsu_thresh

    return thresh


def _getWormContours(ROI_image, threshold, strel_size=(5, 5), is_light_background=True):
    # get the border of the ROI mask, this will be used to filter for valid
    # worms
    ROI_valid = (ROI_image != 0).astype(np.uint8)
    _, ROI_border_ind, _ = cv2.findContours(
        ROI_valid, cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)

    if len(ROI_border_ind) <= 1:
        valid_ind = 0
    else:
        # consider the case where there is more than one contour in the blob
        # i.e. there is a neighboring ROI in the square, just keep the largest area
        ROI_area = [cv2.contourArea(x) for x in ROI_border_ind]
        valid_ind = np.argmax(ROI_area)
        ROI_valid = np.zeros_like(ROI_valid)
        ROI_valid = cv2.drawContours(
            ROI_valid, ROI_border_ind, valid_ind, 1, -1)
        ROI_image = ROI_image * ROI_valid

    # get binary image, 
    ROI_mask = ROI_image < threshold if is_light_background else ROI_image > threshold
    ROI_mask = (ROI_mask & (ROI_image != 0)).astype(np.uint8)
    
    # clean it using morphological closing - make this optional by setting strel_size to 0
    if np.all(strel_size):
        strel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, strel_size)
        ROI_mask = cv2.morphologyEx(ROI_mask, cv2.MORPH_CLOSE, strel)

    # get worms, assuming each contour in the ROI is a worm
    [_, ROI_worms, hierarchy] = cv2.findContours(
        ROI_mask.copy(), cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)

    return ROI_worms, hierarchy


def _getWormFeatures(
        worm_cnt,
        ROI_image,
        ROI_bbox,
        current_frame,
        thresh,
        min_area):
    SEGWORM_ID_DEFAULT = -1  # default value for the column segworm_id

    area = float(cv2.contourArea(worm_cnt))
    if area < min_area:
        return None  # area too small to be a worm

    worm_bbox = cv2.boundingRect(worm_cnt)
    # find use the best rotated bounding box, the fitEllipse function produces bad results quite often
    # this method is better to obtain an estimate of the worm length than
    # eccentricity
    (CMx, CMy), (L, W), angle = cv2.minAreaRect(worm_cnt)
    
    if L == 0 or W == 0:
        return None #something went wrong abort
    
    if W > L:
        L, W = W, L  # switch if width is larger than length
    quirkiness = sqrt(1 - W**2 / L**2)

    hull = cv2.convexHull(worm_cnt)  # for the solidity
    solidity = area / cv2.contourArea(hull)
    perimeter = float(cv2.arcLength(worm_cnt, True))
    compactness = area / (4 * np.pi * perimeter**2)

    # calculate the mean intensity of the worm
    worm_mask = np.zeros(ROI_image.shape, dtype=np.uint8)
    cv2.drawContours(worm_mask, [worm_cnt], 0, 255, -1)
    intensity_mean, intensity_std = cv2.meanStdDev(ROI_image, mask=worm_mask)

    # calculate hu moments, they are scale and rotation invariant
    hu_moments = cv2.HuMoments(cv2.moments(worm_cnt))

    # save everything into the the proper output format
    mask_feat = (
        current_frame,
        CMx +
        ROI_bbox[0],
        CMy +
        ROI_bbox[1],
        area,
        perimeter,
        L,
        W,
        quirkiness,
        compactness,
        angle,
        solidity,
        intensity_mean[
            0,
            0],
        intensity_std[<commit_msg>Corrected feature calculation for compactness to be the isoperimetric quotient correctly, namely 4 pi A / L^2 (previously the 4 pi was in the denominator).<commit_after># -*- coding: utf-8 -*-
"""
Created on Thu Apr  2 16:33:34 2015

@author: ajaver
"""


import numpy as np
import tables
from math import sqrt
import cv2
from skimage.filters import threshold_otsu
import os
import sys

from sklearn.utils.linear_assignment_ import linear_assignment  # hungarian algorithm
from scipy.spatial.distance import cdist

from MWTracker.helperFunctions.timeCounterStr import timeCounterStr
from MWTracker.compressVideos.extractMetaData import readAndSaveTimestamp

table_filters = tables.Filters(complevel=5, complib='zlib', shuffle=True)


class plate_worms(tables.IsDescription):
    # class for the pytables
    worm_index_blob = tables.Int32Col(pos=0)
    worm_index_joined = tables.Int32Col(pos=1)

    frame_number = tables.Int32Col(pos=2)
    coord_x = tables.Float32Col(pos=3)
    coord_y = tables.Float32Col(pos=4)
    area = tables.Float32Col(pos=5)
    perimeter = tables.Float32Col(pos=6)
    box_length = tables.Float32Col(pos=7)
    box_width = tables.Float32Col(pos=8)
    quirkiness = tables.Float32Col(pos=9)
    compactness = tables.Float32Col(pos=10)
    box_orientation = tables.Float32Col(pos=11)
    solidity = tables.Float32Col(pos=12)
    intensity_mean = tables.Float32Col(pos=13)
    intensity_std = tables.Float32Col(pos=14)

    threshold = tables.Int32Col(pos=15)
    bounding_box_xmin = tables.Int32Col(pos=16)
    bounding_box_xmax = tables.Int32Col(pos=17)
    bounding_box_ymin = tables.Int32Col(pos=18)
    bounding_box_ymax = tables.Int32Col(pos=19)

    # deprecated, probably it would be good to remove it in the future
    segworm_id = tables.Int32Col(pos=20)

    hu0 = tables.Float32Col(pos=21)
    hu1 = tables.Float32Col(pos=22)
    hu2 = tables.Float32Col(pos=23)
    hu3 = tables.Float32Col(pos=24)
    hu4 = tables.Float32Col(pos=25)
    hu5 = tables.Float32Col(pos=26)
    hu6 = tables.Float32Col(pos=27)


def _getWormThreshold(pix_valid):
    
    # calculate otsu_threshold as lower limit. Otsu understimates the threshold.
    try:
        otsu_thresh = threshold_otsu(pix_valid)
    except:
        return np.nan

    # calculate the histogram
    pix_hist = np.bincount(pix_valid)

    # the higher limit is the most frequent value in the distribution
    # (background)
    largest_peak = np.argmax(pix_hist)
    if otsu_thresh < largest_peak and otsu_thresh + 2 < len(pix_hist) - 1:
        # smooth the histogram to find a better threshold
        pix_hist = np.convolve(pix_hist, np.ones(3), 'same')
        cumhist = np.cumsum(pix_hist)

        xx = np.arange(otsu_thresh, cumhist.size)
        try:
            # the threshold is calculated as the pixel level where there would be
            # larger increase in the object area.
            hist_ratio = pix_hist[xx] / cumhist[xx]
            thresh = np.where(
                (hist_ratio[3:] - hist_ratio[:-3]) > 0)[0][0] + otsu_thresh
        except IndexError:
            thresh = np.argmin(
                pix_hist[
                    otsu_thresh:largest_peak]) + otsu_thresh
    else:
        # if otsu is larger than the maximum peak keep otsu threshold
        thresh = otsu_thresh

    return thresh


def _getWormContours(ROI_image, threshold, strel_size=(5, 5), is_light_background=True):
    # get the border of the ROI mask, this will be used to filter for valid
    # worms
    ROI_valid = (ROI_image != 0).astype(np.uint8)
    _, ROI_border_ind, _ = cv2.findContours(
        ROI_valid, cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)

    if len(ROI_border_ind) <= 1:
        valid_ind = 0
    else:
        # consider the case where there is more than one contour in the blob
        # i.e. there is a neighboring ROI in the square, just keep the largest area
        ROI_area = [cv2.contourArea(x) for x in ROI_border_ind]
        valid_ind = np.argmax(ROI_area)
        ROI_valid = np.zeros_like(ROI_valid)
        ROI_valid = cv2.drawContours(
            ROI_valid, ROI_border_ind, valid_ind, 1, -1)
        ROI_image = ROI_image * ROI_valid

    # get binary image, 
    ROI_mask = ROI_image < threshold if is_light_background else ROI_image > threshold
    ROI_mask = (ROI_mask & (ROI_image != 0)).astype(np.uint8)
    
    # clean it using morphological closing - make this optional by setting strel_size to 0
    if np.all(strel_size):
        strel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, strel_size)
        ROI_mask = cv2.morphologyEx(ROI_mask, cv2.MORPH_CLOSE, strel)

    # get worms, assuming each contour in the ROI is a worm
    [_, ROI_worms, hierarchy] = cv2.findContours(
        ROI_mask.copy(), cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)

    return ROI_worms, hierarchy


def _getWormFeatures(
        worm_cnt,
        ROI_image,
        ROI_bbox,
        current_frame,
        thresh,
        min_area):
    SEGWORM_ID_DEFAULT = -1  # default value for the column segworm_id

    area = float(cv2.contourArea(worm_cnt))
    if area < min_area:
        return None  # area too small to be a worm

    worm_bbox = cv2.boundingRect(worm_cnt)
    # find use the best rotated bounding box, the fitEllipse function produces bad results quite often
    # this method is better to obtain an estimate of the worm length than
    # eccentricity
    (CMx, CMy), (L, W), angle = cv2.minAreaRect(worm_cnt)
    
    if L == 0 or W == 0:
        return None #something went wrong abort
    
    if W > L:
        L, W = W, L  # switch if width is larger than length
    quirkiness = sqrt(1 - W**2 / L**2)

    hull = cv2.convexHull(worm_cnt)  # for the solidity
    solidity = area / cv2.contourArea(hull)
    perimeter = float(cv2.arcLength(worm_cnt, True))
    compactness = 4 * np.pi * area / (perimeter**2)

    # calculate the mean intensity of the worm
    worm_mask = np.zeros(ROI_image.shape, dtype=np.uint8)
    cv2.drawContours(worm_mask, [worm_cnt], 0, 255, -1)
    intensity_mean, intensity_std = cv2.meanStdDev(ROI_image, mask=worm_mask)

    # calculate hu moments, they are scale and rotation invariant
    hu_moments = cv2.HuMoments(cv2.moments(worm_cnt))

    # save everything into the the proper output format
    mask_feat = (
        current_frame,
        CMx +
        ROI_bbox[0],
        CMy +
        ROI_bbox[1],
        area,
        perimeter,
        L,
        W,
        quirkiness,
        compactness,
        angle,
        solidity,
        intensity_mean[
            0,
            0],
        intensity_std[
{'fraction_non_alphanumeric': 0.08614059375480695, 'fraction_numerical': 0.028226426703584065, 'mean_word_length': 2.859602255862274, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3943158', 'n_tokens_mistral': 4713, 'n_tokens_neox': 4361, 'n_words': 1120}