particleCharacterization.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
GEPARD - Gepard-Enabled PARticle Detection
Copyright (C) 2018  Lars Bittrich and Josef Brandt, Leibniz-Institut für 
Polymerforschung Dresden e. V. <bittrich-lars@ipfdd.de>    

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program, see COPYING.  
If not, see <https://www.gnu.org/licenses/>.
"""

import numpy as np
import cv2
from copy import deepcopy

from .particleClassification.colorClassification import ColorClassifier
from .particleClassification.shapeClassification import ShapeClassifier
from .ftirAperture import findRotatedMaximalRectangle
from ..segmentation import closeHolesOfSubImage
from ..errors import InvalidParticleError
from ..helperfunctions import cv2imread_fix


class FTIRAperture(object):
    """
    Configuration for an FTIR aperture. CenterCoords, width and height in Pixels
    """
    centerX: float = None
    centerY: float = None
    centerZ: float = None
    width: float = None
    height: float = None
    angle: float = None
    rectCoords: list = None


class ParticleStats(object):
    longSize: float = None
    shortSize: float = None
    height: float = None
    area: float = None
    shape: str = None
    color: str = None
    aperture: FTIRAperture = None


def particleIsValid(particle):
    if particle.longSize == 0 or particle.shortSize == 0:
        return False
    
    if cv2.contourArea(particle.contour) == 0:
        return False
    return True


def updateStatsOfParticlesIfNotManuallyEdited(particleContainer):
    dataset = particleContainer.datasetParent
    fullimage = loadFullimageFromDataset(dataset)
    zimg = loadZValImageFromDataset(dataset)
    
    for particle in particleContainer.particles:
        if not hasattr(particle, "wasManuallyEdited"):
            particle.wasManuallyEdited = False
        
        if not particle.wasManuallyEdited:
            newStats = getParticleStatsWithPixelScale(particle.contour, dataset, fullimage, zimg)
            particle.__dict__.update(newStats.__dict__)


def getFTIRAperture(partImg: np.ndarray) -> FTIRAperture:
    rectPoints, angle, width, heigth = findRotatedMaximalRectangle(partImg, nbre_angle=4, limit_image_size=300)
    xvalues: list = [i[0] for i in rectPoints]
    yvalues: list = [i[1] for i in rectPoints]

    aperture = FTIRAperture()
    aperture.height = heigth
    aperture.width = width
    aperture.angle = angle
    aperture.centerX = int(round(min(xvalues) + (max(xvalues) - min(xvalues)) / 2))
    aperture.centerY = int(round(min(yvalues) + (max(yvalues) - min(yvalues)) / 2))
    aperture.rectCoords = rectPoints

    return aperture


def getParticleStatsWithPixelScale(contour, dataset, fullimage=None, zimg=None):
    if fullimage is None:
        fullimage = loadFullimageFromDataset(dataset)
    if zimg is None:
        zimg = loadZValImageFromDataset(dataset)
        
    cnt = deepcopy(contour)
    pixelscale = dataset.getPixelScale()
    
    newStats = ParticleStats()
    newStats.longSize, newStats.shortSize, newStats.area = getContourStats(cnt)
    newStats.longSize *= pixelscale
    newStats.shortSize *= pixelscale
    newStats.area *= (pixelscale**2)
    
    if 0 in [newStats.longSize, newStats.shortSize, newStats.area]:
        raise InvalidParticleError
    
    newStats.height = getParticleHeight(cnt, zimg, dataset)    
    newStats.shape = getParticleShape(cnt, newStats.height)
    if newStats.shape == 'fibre':
        newStats.longSize, newStats.shortSize = getFibreDimension(cnt)
        newStats.longSize *= pixelscale
        newStats.shortSize *= pixelscale

    partImg, extrema = getParticleImageFromFullimage(cnt, fullimage)
    newStats.color = getParticleColor(partImg)

    padding: int = int(2)
    imgforAperture = np.zeros((partImg.shape[0]+2*padding, partImg.shape[1]+2*padding))
    imgforAperture[padding:partImg.shape[0]+padding, padding:partImg.shape[1]+padding] = np.mean(partImg, axis=2)
    imgforAperture[imgforAperture > 0] = 1  # convert to binary image
    imgforAperture = np.uint8(1 - imgforAperture)  # background has to be 1, particle has to be 0
    aperture: FTIRAperture = getFTIRAperture(imgforAperture)
    aperture.centerX = aperture.centerX + extrema[0] - padding
    aperture.centerY = aperture.centerY + extrema[2] - padding
    aperture.centerZ = newStats.height
    newStats.aperture = aperture

    return newStats


def getFibreDimension(contour):
    longSize = cv2.arcLength(contour, True)/2
    img = contoursToImg([contour])[0]
    dist = cv2.distanceTransform(img, cv2.DIST_L2, 3)
    maxThickness = np.max(dist)*2
    return longSize, maxThickness


def getParticleColor(imgRGB, colorClassifier=None):
    img = cv2.cvtColor(imgRGB, cv2.COLOR_RGB2HSV_FULL)
    meanHSV = cv2.mean(img)
    if colorClassifier is None:
        colorClassifier = ColorClassifier()
    color = colorClassifier.classifyColor(meanHSV)
    return color


def getParticleShape(contour, particleHeight, shapeClassifier=None):
    if shapeClassifier is None:
        shapeClassifier = ShapeClassifier()
    shape = shapeClassifier.classifyShape(contour, particleHeight)
    return shape


def getParticleHeight(contour, fullZimg, dataset):
    zimg, _extrema = getParticleImageFromFullimage(contour, fullZimg)
    if zimg.shape[0] == 0 or zimg.shape[1] == 0:
        raise InvalidParticleError
        
    zimg = cv2.medianBlur(zimg, 5)
    avg_ZValue = np.mean(zimg[zimg > 0])
    if np.isnan(avg_ZValue):  #i.e., only zeros in zimg
        avg_ZValue = 0
    z0, z1 = dataset.zpositions.min(), dataset.zpositions.max()
    height = avg_ZValue/255.*(z1-z0) + z0
    return height


def getContourStats(cnt):
    ##characterize particle
        if cnt.shape[0] >= 5:       ##at least 5 points required for ellipse fitting...
            ellipse = cv2.fitEllipse(cnt)
            short, long = ellipse[1]
        else:
            rect = cv2.minAreaRect(cnt)
            long, short = rect[1]
            
        if short>long:
            long, short = short, long
        
        area = cv2.contourArea(cnt)
    
        return long, short, area


def mergeContours(contours):
    img, xmin, ymin, padding = contoursToImg(contours)
    return imgToCnt(img, xmin, ymin, padding)


def getParticleImageFromFullimage(contour, fullimage):
    contourCopy = deepcopy(contour)
    xmin, xmax, ymin, ymax = getContourExtrema(contourCopy)

    img = fullimage[ymin:ymax, xmin:xmax]
    img = img.copy()
    mask = np.zeros(img.shape[:2])
    
    for i in range(len(contourCopy)):
        contourCopy[i][0][0] -= xmin
        contourCopy[i][0][1] -= ymin
        
    cv2.drawContours(mask, [contourCopy], -1, (255, 255, 255), -1)
    cv2.drawContours(mask, [contourCopy], -1, (255, 255, 255), 1)
    
    img[mask == 0] = 0
    img = np.array(img, dtype = np.uint8)
    return img, (xmin, xmax, ymin, ymax)


def contoursToImg(contours, padding=0):
    contourCopy = deepcopy(contours)
    xmin, xmax, ymin, ymax = getContourExtrema(contourCopy)
    
    padding = padding   #pixel in each direction
    rangex = int(np.round((xmax-xmin)+2*padding))
    rangey = int(np.round((ymax-ymin)+2*padding))
    if rangex == 0 or rangey == 0:
        raise InvalidParticleError

    img = np.zeros((rangey, rangex))
    for curCnt in contourCopy:
        for i in range(len(curCnt)):
            curCnt[i][0][0] -= xmin-padding
            curCnt[i][0][1] -= ymin-padding
        cv2.drawContours(img, [curCnt], -1, 255, -1)
        cv2.drawContours(img, [curCnt], -1, 255, 1)
    img = np.uint8(cv2.morphologyEx(img, cv2.MORPH_CLOSE, np.ones((3, 3))))
    return img, xmin, ymin, padding


def imgToCnt(img, xmin, ymin, padding=0):
    def getContour(img, contourMode):
        if cv2.__version__ > '3.5':
            contours, hierarchy = cv2.findContours(img, cv2.RETR_EXTERNAL, contourMode)
        else:
            temp, contours, hierarchy = cv2.findContours(img, cv2.RETR_EXTERNAL, contourMode)
        
        if len(contours) == 0:   #i.e., no contour found
            raise InvalidParticleError
        elif len(contours) == 1:   #i.e., exactly one contour found
            contour = contours[0]
        else:       #i.e., multiple contours found
            contour = getLargestContour(contours)
            
        return contour
    
    def getLargestContour(contours):
        areas = []
        for contour in contours:
            areas.append(cv2.contourArea(contour))
        maxIndex = areas.index(max(areas))
        print(f'{len(contours)} contours found, getting the largest one. Areas are: {areas}, taking contour at index {maxIndex}')
        return contours[maxIndex]

    img = closeHolesOfSubImage(img)
    contour = getContour(img, contourMode=cv2.CHAIN_APPROX_NONE)
    
    for i in range(len(contour)):
        contour [i][0][0] += xmin-padding
        contour [i][0][1] += ymin-padding
        
    return contour


def getContourExtrema(contours):
    try:
        cnt = np.vstack(tuple(contours))
        xmin, xmax = cnt[:,0,:][:, 0].min(), cnt[:,0,:][:, 0].max()
        ymin, ymax = cnt[:,0,:][:, 1].min(), cnt[:,0,:][:, 1].max()        
    except IndexError: #i.e., not a list of contours was passed, but an individual contour. Hence, the above indexing does not work
        xmin, xmax = cnt[:, 0].min(), cnt[:, 0].max()
        ymin, ymax = cnt[:, 1].min(), cnt[:, 1].max()
        
    xmin, xmax = int(round(xmin)), int(round(xmax))
    ymin, ymax = int(round(ymin)), int(round(ymax))
    
    return xmin, xmax, ymin, ymax


def getParticleCenterPoint(contour):
    img, xmin, ymin, padding = contoursToImg(contour)
    dist = cv2.distanceTransform(img, cv2.DIST_L2, 3)
    ind = np.argmax(dist)
    y = ind//dist.shape[1]-1
    x = ind%dist.shape[1]-1
    x += xmin
    y += ymin
    return x, y


def loadFullimageFromDataset(dataset):
    return cv2imread_fix(dataset.getImageName())


def loadZValImageFromDataset(dataset):
    return cv2imread_fix(dataset.getZvalImageName(), cv2.IMREAD_GRAYSCALE)