particleCharacterization.py 11.9 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
26
from typing import TYPE_CHECKING
JosefBrandt's avatar
JosefBrandt committed
27

28 29
from .particleClassification.colorClassification import ColorClassifier
from .particleClassification.shapeClassification import ShapeClassifier
Josef Brandt's avatar
Josef Brandt committed
30
from .ftirAperture import findRotatedMaximalRectangle
31 32
from ..segmentation import closeHolesOfSubImage
from ..errors import InvalidParticleError
33
from ..helperfunctions import cv2imread_fix
34

35 36 37 38
if TYPE_CHECKING:
    from ..analysis.particleAndMeasurement import Particle
    from ..scenePyramid import ScenePyramid

Josef Brandt's avatar
Josef Brandt committed
39 40 41 42 43

class FTIRAperture(object):
    """
    Configuration for an FTIR aperture. CenterCoords, width and height in Pixels
    """
44 45 46 47 48 49
    centerX: float = None  # in px
    centerY: float = None  # in px
    centerZ: float = None  # in px
    width: float = None  # in px
    height: float = None  # in px
    angle: float = None  # in degree
Josef Brandt's avatar
Josef Brandt committed
50 51 52
    rectCoords: list = None


53
class ParticleStats(object):
Josef Brandt's avatar
Josef Brandt committed
54 55 56 57 58 59 60 61
    longSize: float = None
    shortSize: float = None
    height: float = None
    area: float = None
    shape: str = None
    color: str = None
    aperture: FTIRAperture = None

62 63 64 65 66 67 68 69

def particleIsValid(particle):
    if particle.longSize == 0 or particle.shortSize == 0:
        return False
    
    if cv2.contourArea(particle.contour) == 0:
        return False
    return True
Josef Brandt's avatar
Josef Brandt committed
70 71


72 73 74 75 76 77 78 79 80 81 82 83 84
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__)

Josef Brandt's avatar
Josef Brandt committed
85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101

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


Josef Brandt's avatar
Josef Brandt committed
102
def getParticleStatsWithPixelScale(contour, dataset, fullimage=None, scenePyramid=None, zimg=None, ftir: bool = False):
103 104 105
    if zimg is None:
        zimg = loadZValImageFromDataset(dataset)
        
106
    cnt = deepcopy(contour)
107
    pixelscale = dataset.getPixelScale()
JosefBrandt's avatar
JosefBrandt committed
108
    
109
    newStats = ParticleStats()
110 111 112 113
    newStats.longSize, newStats.shortSize, newStats.area = getContourStats(cnt)
    newStats.longSize *= pixelscale
    newStats.shortSize *= pixelscale
    newStats.area *= (pixelscale**2)
114 115 116 117
    
    if 0 in [newStats.longSize, newStats.shortSize, newStats.area]:
        raise InvalidParticleError
    
118
    newStats.height = getParticleHeight(cnt, zimg, dataset)    
119
    newStats.shape = getParticleShape(cnt, newStats.height)
JosefBrandt's avatar
JosefBrandt committed
120 121 122 123
    if newStats.shape == 'fibre':
        newStats.longSize, newStats.shortSize = getFibreDimension(cnt)
        newStats.longSize *= pixelscale
        newStats.shortSize *= pixelscale
124

Josef Brandt's avatar
Josef Brandt committed
125
    partImg = None
126 127
    assert scenePyramid is not None or fullimage is not None

Josef Brandt's avatar
Josef Brandt committed
128 129 130 131
    if scenePyramid is None and fullimage is not None:
        partImg, extrema = getParticleImageFromFullimage(cnt, fullimage)
    elif scenePyramid is not None and fullimage is None:
        partImg, extrema = getParticleImageFromScenePyramid(cnt, scenePyramid)
132

Josef Brandt's avatar
Josef Brandt committed
133
    assert partImg is not None, "error in getting particle image"
134
    newStats.color = getParticleColor(partImg)
Josef Brandt's avatar
Josef Brandt committed
135

Josef Brandt's avatar
Josef Brandt committed
136
    if ftir:
137
        newStats.aperture = getApertureObjectForParticleImage(partImg, extrema, newStats.height)
Josef Brandt's avatar
Josef Brandt committed
138

139
    return newStats
140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158


def calculateFTIRAperturesForParticle(particle: 'Particle', pyramid: 'ScenePyramid') -> 'FTIRAperture':
    partImg, extrema = getParticleImageFromScenePyramid(particle.contour, pyramid)
    aperture: 'FTIRAperture' = getApertureObjectForParticleImage(partImg, extrema, particle.height)
    return aperture


def getApertureObjectForParticleImage(particleImage: np.ndarray, extrema: tuple,
                                      particleHeight: float, padding: int = 2) -> 'FTIRAperture':
    imgforAperture = np.zeros((particleImage.shape[0] + 2 * padding, particleImage.shape[1] + 2 * padding))
    imgforAperture[padding:particleImage.shape[0] + padding, padding:particleImage.shape[1] + padding] = np.mean(particleImage, 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 = particleHeight
    return aperture
159

Josef Brandt's avatar
Josef Brandt committed
160

JosefBrandt's avatar
JosefBrandt committed
161 162 163 164 165 166
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
Josef Brandt's avatar
Josef Brandt committed
167 168


169 170 171 172 173 174 175 176
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

Josef Brandt's avatar
Josef Brandt committed
177

178 179 180
def getParticleShape(contour, particleHeight, shapeClassifier=None):
    if shapeClassifier is None:
        shapeClassifier = ShapeClassifier()
JosefBrandt's avatar
JosefBrandt committed
181
    shape = shapeClassifier.classifyShape(contour, particleHeight)
182 183
    return shape

Josef Brandt's avatar
Josef Brandt committed
184

185
def getParticleHeight(contour, fullZimg, dataset):
Josef Brandt's avatar
Josef Brandt committed
186
    zimg, _extrema = getParticleImageFromFullimage(contour, fullZimg)
187 188 189 190 191
    if zimg.shape[0] == 0 or zimg.shape[1] == 0:
        raise InvalidParticleError
        
    zimg = cv2.medianBlur(zimg, 5)
    avg_ZValue = np.mean(zimg[zimg > 0])
Josef Brandt's avatar
Josef Brandt committed
192
    if np.isnan(avg_ZValue):  # i.e., only zeros in zimg
193 194 195 196 197
        avg_ZValue = 0
    z0, z1 = dataset.zpositions.min(), dataset.zpositions.max()
    height = avg_ZValue/255.*(z1-z0) + z0
    return height

Josef Brandt's avatar
Josef Brandt committed
198

199
def getContourStats(cnt):
Josef Brandt's avatar
Josef Brandt committed
200 201 202 203 204 205 206 207 208 209 210 211 212 213
    # 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
214

Josef Brandt's avatar
Josef Brandt committed
215

216 217 218 219
def mergeContours(contours):
    img, xmin, ymin, padding = contoursToImg(contours)
    return imgToCnt(img, xmin, ymin, padding)

Josef Brandt's avatar
Josef Brandt committed
220

221 222 223 224 225
def getParticleImageFromFullimage(contour, fullimage):
    contourCopy = deepcopy(contour)
    xmin, xmax, ymin, ymax = getContourExtrema(contourCopy)

    img = fullimage[ymin:ymax, xmin:xmax]
226
    img = img.copy()
227 228 229 230 231 232 233 234 235 236 237
    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)
Josef Brandt's avatar
Josef Brandt committed
238 239 240
    return img, (xmin, xmax, ymin, ymax)


Josef Brandt's avatar
Josef Brandt committed
241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
def getParticleImageFromScenePyramid(contour, scenePyramid):
    contourCopy = deepcopy(contour)
    xmin, xmax, ymin, ymax = getContourExtrema(contourCopy)

    img = scenePyramid.getImagePart(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)


261
def contoursToImg(contours, padding=0):
262 263
    contourCopy = deepcopy(contours)
    xmin, xmax, ymin, ymax = getContourExtrema(contourCopy)
Josef Brandt's avatar
Josef Brandt committed
264

265 266
    rangex = int(np.round((xmax-xmin)+2*padding))
    rangey = int(np.round((ymax-ymin)+2*padding))
267 268
    if rangex == 0 or rangey == 0:
        raise InvalidParticleError
269 270 271 272 273 274 275 276 277 278 279

    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

Josef Brandt's avatar
Josef Brandt committed
280

281
def imgToCnt(img, xmin, ymin, padding=0):
Josef Brandt's avatar
Josef Brandt committed
282

283
    def getContour(img, contourMode):
Josef Brandt's avatar
Josef Brandt committed
284 285 286 287 288 289 290 291 292 293

        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}, '
                  f'taking contour at index {maxIndex}')
            return contours[maxIndex]

294
        if cv2.__version__ > '3.5':
295
            contours, hierarchy = cv2.findContours(img, cv2.RETR_EXTERNAL, contourMode)
296
        else:
297
            temp, contours, hierarchy = cv2.findContours(img, cv2.RETR_EXTERNAL, contourMode)
298
        
Josef Brandt's avatar
Josef Brandt committed
299
        if len(contours) == 0:  # i.e., no contour found
300
            raise InvalidParticleError
Josef Brandt's avatar
Josef Brandt committed
301
        elif len(contours) == 1:  # i.e., exactly one contour found
302
            contour = contours[0]
Josef Brandt's avatar
Josef Brandt committed
303
        else:  # i.e., multiple contours found
304 305
            contour = getLargestContour(contours)
            
306
        return contour
307 308

    img = closeHolesOfSubImage(img)
309
    contour = getContour(img, contourMode=cv2.CHAIN_APPROX_NONE)
Josef Brandt's avatar
Josef Brandt committed
310

311
    for i in range(len(contour)):
Josef Brandt's avatar
Josef Brandt committed
312 313 314
        contour[i][0][0] += xmin-padding
        contour[i][0][1] += ymin-padding

315
    return contour
316

Josef Brandt's avatar
Josef Brandt committed
317

318 319 320 321 322 323 324
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()
325 326 327 328 329
        ymin, ymax = cnt[:, 1].min(), cnt[:, 1].max()
        
    xmin, xmax = int(round(xmin)), int(round(xmax))
    ymin, ymax = int(round(ymin)), int(round(ymax))
    
330 331
    return xmin, xmax, ymin, ymax

Josef Brandt's avatar
Josef Brandt committed
332

333 334 335 336 337 338 339 340 341
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
342

Josef Brandt's avatar
Josef Brandt committed
343

344
def loadFullimageFromDataset(dataset):
345
    return cv2imread_fix(dataset.getImageName())
346

Josef Brandt's avatar
Josef Brandt committed
347

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