particleCharacterization.py 7.58 KB
Newer Older
JosefBrandt's avatar
JosefBrandt committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
#!/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
25
from copy import deepcopy
JosefBrandt's avatar
JosefBrandt committed
26

27 28
from .particleClassification.colorClassification import ColorClassifier
from .particleClassification.shapeClassification import ShapeClassifier
29 30
from ..segmentation import closeHolesOfSubImage
from ..errors import InvalidParticleError
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

class ParticleStats(object):
    longSize = None
    shortSize = None
    height = None
    area = None
    shape = None
    color = None

def particleIsValid(particle):
    if particle.longSize == 0 or particle.shortSize == 0:
        return False
    
    if cv2.contourArea(particle.contour) == 0:
        return False
    return True
    
48 49
def getParticleStatsWithPixelScale(contour, fullimage, dataset):
    cnt = deepcopy(contour)
50
    pixelscale = dataset.getPixelScale()
JosefBrandt's avatar
JosefBrandt committed
51
    
52
    newStats = ParticleStats()
53 54 55 56
    newStats.longSize, newStats.shortSize, newStats.area = getContourStats(cnt)
    newStats.longSize *= pixelscale
    newStats.shortSize *= pixelscale
    newStats.area *= (pixelscale**2)
57 58 59 60
    
    if 0 in [newStats.longSize, newStats.shortSize, newStats.area]:
        raise InvalidParticleError
    
61 62
    newStats.height = getParticleHeight(cnt, dataset)    
    newStats.shape = getParticleShape(cnt, newStats.height)
JosefBrandt's avatar
JosefBrandt committed
63 64 65 66
    if newStats.shape == 'fibre':
        newStats.longSize, newStats.shortSize = getFibreDimension(cnt)
        newStats.longSize *= pixelscale
        newStats.shortSize *= pixelscale
67

68 69 70
    partImg = getParticleImageFromFullimage(cnt, fullimage)
    newStats.color = getParticleColor(partImg)
    return newStats
71

JosefBrandt's avatar
JosefBrandt committed
72 73 74 75 76 77 78
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
    
79 80 81 82 83 84 85 86
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

87 88 89
def getParticleShape(contour, particleHeight, shapeClassifier=None):
    if shapeClassifier is None:
        shapeClassifier = ShapeClassifier()
JosefBrandt's avatar
JosefBrandt committed
90
    shape = shapeClassifier.classifyShape(contour, particleHeight)
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
    return shape

def getParticleHeight(contour, dataset):
    zimg = getParticleImageFromFullimage(contour, dataset.getZvalImg())
    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

122 123 124 125 126 127 128 129 130
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]
131
    img = img.copy()
132 133 134 135 136 137 138 139 140 141 142 143 144
    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
    
145
def contoursToImg(contours, padding=0):
146 147 148 149 150 151
    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))
152 153
    if rangex == 0 or rangey == 0:
        raise InvalidParticleError
154 155 156 157 158 159 160 161 162 163 164

    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

165
def imgToCnt(img, xmin, ymin, padding=0):
166
    def getContour(img, contourMode):
167
        if cv2.__version__ > '3.5':
168
            contours, hierarchy = cv2.findContours(img, cv2.RETR_EXTERNAL, contourMode)
169
        else:
170
            temp, contours, hierarchy = cv2.findContours(img, cv2.RETR_EXTERNAL, contourMode)
171
        
172 173 174 175 176 177 178
        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)
            
179 180
        return contour
    
181 182 183 184 185 186 187 188 189
    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)
190
    contour = getContour(img, contourMode=cv2.CHAIN_APPROX_NONE)
191
    
192 193 194
    for i in range(len(contour)):
        contour [i][0][0] += xmin-padding
        contour [i][0][1] += ymin-padding
195
        
196
    return contour
197 198 199 200 201 202 203 204

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()
205 206 207 208 209
        ymin, ymax = cnt[:, 1].min(), cnt[:, 1].max()
        
    xmin, xmax = int(round(xmin)), int(round(xmax))
    ymin, ymax = int(round(ymin)), int(round(ymax))
    
210 211
    return xmin, xmax, ymin, ymax

212 213 214 215 216 217 218 219 220
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