Commit 0ddfc91d authored by Josef Brandt's avatar Josef Brandt

Particle Image Features, Curvature DFT Features

parent 7ca2dd08
......@@ -61,34 +61,53 @@ def get_n_points_closest_to_point(points: np.ndarray, n: int, refPoint: np.ndarr
return list(sortedIndices[:n])
def get_particle_featurematrix(particleContainer: ParticleContainer) -> np.ndarray:
def get_particle_featurematrix(particleContainer: ParticleContainer, fullimg: np.ndarray = None) -> np.ndarray:
"""
:return: np.ndarray, numRows: Features, numCols: Particles
"""
vectors: list = []
for particle in particleContainer.particles:
vectors.append(get_characteristic_vector(particle))
firstVecLength: int = 0
for index, particle in enumerate(particleContainer.particles):
particleImg: np.ndarray = None
if fullimg is not None:
particleImg = pc.getParticleImageFromFullimage(particle.contour, fullimg)
vector: np.ndarray = get_characteristic_vector(particle, particleImg)
if index == 0:
firstVecLength = vector.shape[0]
else:
assert vector.shape[0] == firstVecLength, 'particle feature vectors have non-uniform lengths...'
vectors.append(vector)
vectors: np.ndarray = np.array(vectors)
assert vectors.shape[0] == len(particleContainer.particles)
return vectors
def get_characteristic_vector(particle: pm.Particle) -> np.ndarray:
def get_characteristic_vector(particle: pm.Particle, particleImg: np.ndarray = None) -> np.ndarray:
vector: list = []
# vector += list(get_log_hu_moments(particle.contour))
vector += list(get_log_hu_moments(particle.contour))
vector.append(float(get_color_index(particle.color)))
vector.append(get_solidity(particle.contour))
vector.append(get_aspect_ratio(particle.contour))
vector.append(get_extent(particle.contour))
vector.append(cv2.contourArea(particle.contour))
# vector.append(get_shape_index(particle.shape))
# vector.append(cv2.arcLength(particle.contour, True))
vector.append(get_shape_index(particle.shape))
vector.append(cv2.arcLength(particle.contour, True))
for ftHarmonic in get_curvature_ft(particle.contour):
vector.append(ftHarmonic)
if particleImg is not None:
image_features: np.ndarray = get_image_feature_vec(particleImg)
for val in image_features:
vector.append(val)
# vector: np.ndarray = np.hstack((log_hu, color))
# if len(vector) != 11:
# print('error')
# assert len(vector) == 7 + 4, f'wrong feature vector: {vector} with shape: {vector.shape}'
for entry in vector:
assert type(entry) in [float, int, np.int, np.float64], f'{entry} is {type(entry)}'
assert not np.isnan(entry)
return np.array(vector)
......@@ -141,6 +160,18 @@ def get_log_hu_moments(contour: np.ndarray) -> np.ndarray:
return resultMoments[:, 0]
def get_image_feature_vec(particleImg: np.ndarray) -> np.ndarray:
meanStdev: np.ndarray = get_mean_and_stdev(particleImg)
return meanStdev.reshape(meanStdev.size)
def get_mean_and_stdev(img: np.ndarray) -> np.ndarray:
meanStd: tuple = cv2.meanStdDev(img)
colorMean: np.ndarray = np.array([i[0] for i in meanStd[0]])
colorStd: np.ndarray = np.array([i[0] for i in meanStd[1]])
return np.vstack((colorMean, colorStd))
def get_color_hash(color: str, desiredLength: int = 4) -> np.ndarray:
colorArray: list = [int(i) for i in str(abs(hash(color)))[:desiredLength]]
return np.transpose(np.array(colorArray))
......@@ -153,10 +184,61 @@ def get_color_index(color: str) -> int:
return colors.index(color)
# def get_shape_index(shape: str) -> int:
# shapes: list = ['spherule', 'fibre', 'flake', 'irregular']
# assert shape in shapes
# return shapes.index(shape)
def get_shape_index(shape: str) -> int:
shapes: list = ['spherule', 'fibre', 'flake', 'irregular', None]
assert shape in shapes
return shapes.index(shape)
def get_curvature_ft(contour: np.ndarray, angularSegment: float = 20, numHarmonics: int = 8) -> np.ndarray:
"""
Calculates curvature of the contour and applies discrete fourier transform (dft) to it.
Adapted from Xu et al, "Comparison of shape features for the classification of wear particles",
doi:10.1016/S0952-1976(97)00017-1
:param angularSegment: In degree, section of 360 degree to consider for each angle
:param contour: Particle Contour
:param numHarmonics: Number of FT Harmonics to return
:return: the computed
"""
m: int = int(round(contour.shape[0] / (360 / angularSegment))) # smoothing factor, is this reasonable? That's, how it's described in the paper...
numPoints: int = contour.shape[0]
curvature: np.ndarray = np.zeros(numPoints)
for i in range(numPoints):
x, y = contour[i, 0, 0], contour[i, 0, 1]
index_before = i-m-1
index_after = i+m+1
if index_after > numPoints-1:
index_after -= numPoints
x_before, y_before = contour[index_before, 0, 0], contour[index_before, 0, 1]
x_after, y_after = contour[index_after, 0, 0], contour[index_after, 0, 1]
tan_before: float = 0
if x - x_before != 0:
tan_before = np.rad2deg(np.arctan((y-y_before) / (x-x_before)))
tan_after: float = 0
if x_after - x != 0:
tan_after: float = np.rad2deg(np.arctan((y_after-y) / (x_after-x)))
curvature[i] = tan_after - tan_before
dft: np.ndarray = np.fft.fft(curvature)
freqs: np.ndarray = np.fft.fftfreq(numPoints)
mask: np.ndarray = freqs > 0
dft_true: np.ndarray = 2 * np.abs(dft/numPoints)
dft_true = dft_true[mask]
numFTPoints = dft_true.shape[0]
# let's only take the odd ones (they seem to carry more information).
# Also, make sure that always numHarmonics numbers come out. Fill up with zeros, if contour was very short
indicesToTake: np.ndarray = np.arange(numHarmonics)*2 + 1
final_dfts = np.zeros(numHarmonics)
for i, indexToTake in enumerate(indicesToTake):
if indexToTake < numFTPoints:
final_dfts[i] = dft_true[indexToTake]
return final_dfts
class TrainedSubsampling(SubsamplingMethod):
......@@ -168,10 +250,6 @@ class TrainedSubsampling(SubsamplingMethod):
self.clfPath: str = path
self.fraction = desiredFraction
# @property
# def fraction(self) -> float:
# return self.desiredFraction/2
def equals(self, otherMethod) -> bool:
isEqual: bool = False
if type(otherMethod) == TrainedSubsampling and otherMethod.fraction == self.fraction:
......@@ -199,9 +277,6 @@ class TrainedSubsampling(SubsamplingMethod):
return selectedParticles
# def _make_subparticles_match_fraction(self, subParticles: list) -> list:
# return subParticles
def _load_classifier(self) -> None:
assert os.path.exists(self.clfPath)
fname: str = self.clfPath
......
......@@ -12,9 +12,12 @@ from sklearn.feature_selection import chi2
import pickle
import time
import sys
import os
sys.path.append("C://Users//xbrjos//Desktop//Python")
from gepard import dataset
from gepard.analysis.particleContainer import ParticleContainer
from gepard.analysis.particleCharacterization import getParticleImageFromFullimage
from gepard.helperfunctions import cv2imread_fix
from input_output import get_pkls_from_directory
from chemometricMethods import get_log_hu_moments, get_color_index, get_pca, get_characteristic_vector
......@@ -58,10 +61,10 @@ def test_classification_models(dataset: tuple) -> None:
if __name__ == '__main__':
recreateNew: bool = False
recreateNew: bool = True
if recreateNew:
pklsInFolders: dict = get_pkls_from_directory(r'C:\Users\xbrjos\Desktop\temp MP\NewDatasets')
fullimgpath: str = r'C:\Users\xbrjos\Desktop\temp MP\Fullimages'
X: list = []
y: list = []
......@@ -71,11 +74,20 @@ if __name__ == '__main__':
if counter < 100:
dset: dataset.DataSet = dataset.loadData(pklPath)
print('loaded', dset.name)
imgPath: str = os.path.join(fullimgpath, dset.name + '.tif')
fullimg = cv2imread_fix(imgPath)
print('loaded fullimg', imgPath)
partContainer: ParticleContainer = dset.particleContainer
for particle in partContainer.particles:
features: np.ndarray = get_characteristic_vector(particle)
# features: list = [abs(i) for i in get_log_hu_moments(particle.contour)]
# features.append(get_color_index(particle.color))
firstVecLength: int = 0
for index, particle in enumerate(partContainer.particles):
partImg: np.ndarray = getParticleImageFromFullimage(particle.contour, fullimg)
features: np.ndarray = get_characteristic_vector(particle, partImg)
if index == 0:
firstVecLength = features.shape[0]
else:
assert features.shape[0] == firstVecLength
X.append(features)
y.append(int(is_MP_particle(particle)))
......@@ -100,6 +112,14 @@ if __name__ == '__main__':
X_equalized: np.ndarray = np.array(X_MP + X_nonMP)
y_equalized: np.ndarray = np.array(y_MP + y_nonMP)
sum = X_MP + X_nonMP
for i in range(len(sum)):
if i== 0:
continue
try:
np.min(np.array(sum[:i]))
except:
print('error')
dset: tuple = (X_equalized, y_equalized)
......@@ -111,10 +131,10 @@ if __name__ == '__main__':
X, y = dset
with open(r'C:\Users\xbrjos\Desktop\Python\Subsampling\chemometrics\RandomForestClassifier, score 0.72.pkl', "rb") as fp:
clf: RandomForestClassifier = pickle.load(fp)
# with open(r'C:\Users\xbrjos\Desktop\Python\Subsampling\chemometrics\RandomForestClassifier, score 0.72.pkl', "rb") as fp:
# clf: RandomForestClassifier = pickle.load(fp)
y_predicted = clf.predict(X)
# y_predicted = clf.predict(X)
# np.savetxt('Data.txt', X)
# np.savetxt('Assignments.txt', y)
......@@ -127,4 +147,4 @@ if __name__ == '__main__':
# print(X_equalized.shape)
# test_classification_models((X, y))
test_classification_models((X, y))
import sys
sys.path.append("C://Users//xbrjos//Desktop//Python")
from gepard.helperfunctions import cv2imread_fix
from gepard.dataset import loadData
from gepard.analysis.particleContainer import ParticleContainer
import cv2
import numpy as np
from scipy import spatial
import os
import matplotlib.pyplot as plt
if __name__ == '__main__':
# imgpath: str = r'C:\Users\xbrjos\Desktop\temp MP\Fullimages'
# imgname: str = '181120_MCI_2_ds1+2_all_ kleiner500_10_1.tif'
# imgname: str = '190619_5_PTPH_sld_190321_ds1_50_1_neu.tif'
#
# img: np.ndarray = cv2imread_fix(os.path.join(imgpath, imgname))
# gray: np.ndarray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
# _, binimg = cv2.threshold(gray, 20, 1, cv2.THRESH_BINARY_INV)
# distmap: np.ndarray = cv2.distanceTransform(binimg, cv2.DIST_L1, 3)
# plt.imshow(distmap, cmap='gray')
paths: list = [r'C:\Users\xbrjos\Desktop\temp MP\NewDatasets\water\181120_MCI_2_ds1+2_all_ kleiner500_10_1.pkl',
r'C:\Users\xbrjos\Desktop\temp MP\NewDatasets\wastewater, slush\190619_5_PTPH_sld_190321_ds1_50_1_neu.pkl']
# paths.append(r'C:\Users\xbrjos\Desktop\temp MP\NewDatasets\water\190306_MCII_1_2_50.pkl')
# paths.append(r'C:\Users\xbrjos\Desktop\temp MP\NewDatasets\water\190222_MCII_1_1_50_1.pkl')
distances: list = []
for path in paths:
dset = loadData(path)
particleContainer: ParticleContainer = dset.particleContainer
particleCenters: list = []
for particle in particleContainer.particles:
particleCenters.append([np.mean(particle.contour[:, 0, 0]), np.mean(particle.contour[:, 0, 1])])
closest_particle_distances: np.ndarray = np.zeros(len(particleCenters))
particleCenters: np.ndarray = np.array(particleCenters)
print('particle centers done')
distMat: np.ndarray = spatial.distance_matrix(particleCenters, particleCenters)
print('distmat computed')
for i in range(distMat.shape[0]):
if i == 0:
closest_particle_distances[i] = np.min(distMat[i, 1:])
elif i == distMat.shape[0]-1:
closest_particle_distances[i] = np.min(distMat[i, :-1])
else:
closest_particle_distances[i] = np.min([np.min(distMat[i, :i]), np.min(distMat[i, i+1:])])
distances.append(closest_particle_distances)
plt.boxplot(distances)
......@@ -21,8 +21,33 @@ from evaluation import SubsamplingResult
class TestParticleFeatures(unittest.TestCase):
def test_image_features(self):
img: np.ndarray = np.zeros((100, 100, 3), dtype=np.uint8)
img[:, :, 0] = 255 # we just have a plain red image
imgFeatureVector: np.ndarray = cmeth.get_image_feature_vec(img)
self.assertEqual(imgFeatureVector.shape[0], 6)
def test_get_mean_color_and_stdev(self):
img: np.ndarray = np.zeros((100, 100, 3), dtype=np.uint8)
img[:, :, 0] = 255 # we just have a plain red image
meanStd = cmeth.get_mean_and_stdev(img)
self.assertTrue(np.array_equal(meanStd[0], np.array([255, 0, 0])))
self.assertTrue(np.array_equal(meanStd[1], np.array([0, 0, 0])))
img[:, :, 1] = 255
meanStd = cmeth.get_mean_and_stdev(img)
self.assertTrue(np.array_equal(meanStd[0], np.array([255, 255, 0])))
self.assertTrue(np.array_equal(meanStd[1], np.array([0, 0, 0])))
img[:50, :50, 0] = 128
meanStd = cmeth.get_mean_and_stdev(img)
meanRed: float = np.mean([128, 255, 255, 255])
stdRed: float = np.std([128, 255, 255, 255])
self.assertTrue(np.array_equal(meanStd[0], np.array([meanRed, 255, 0])))
self.assertTrue(np.array_equal(meanStd[1], np.array([stdRed, 0, 0])))
def test_get_contour_moments(self):
imgs = []
imgs: list = []
imgA: np.ndarray = np.zeros((200, 200), dtype=np.uint8)
cv2.putText(imgA, 'A', (25, 175), fontFace=cv2.FONT_HERSHEY_SIMPLEX,
fontScale=7, color=1, thickness=5)
......@@ -51,6 +76,25 @@ class TestParticleFeatures(unittest.TestCase):
diff: np.ndarray = moments[i, :] - np.mean(moments[i, :])
self.assertFalse(np.any(diff > 0.1))
def test_get_curvature_ft(self):
def get_cnt_from_img(binImg: np.ndarray) -> np.ndarray:
contours, hierarchy = cv2.findContours(binImg, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_NONE)
return contours[0]
harmonics: list = []
for shape in [cv2.MORPH_CROSS, cv2.MORPH_ELLIPSE, cv2.MORPH_RECT]:
for size in [50, 500]:
padding = round(size/10)
img: np.ndarray = np.zeros((size + 2*padding, size + 2*padding))
img[padding:size+padding, padding:size+padding] = cv2.getStructuringElement(shape, (size, size))
cnt: np.ndarray = get_cnt_from_img(img.astype(np.uint8))
for numHarmonics in [2, 5, 15]:
dft: np.ndarray = cmeth.get_curvature_ft(cnt, numHarmonics=numHarmonics)
harmonics.append(dft)
self.assertEqual(dft.shape[0], numHarmonics)
def test_get_color_hash(self):
for color in ['red', 'green', 'violet', 'blue', 'Blue', 'non-determinable', None]:
for numDigits in [4, 6, 8]:
......@@ -70,6 +114,8 @@ class TestParticleFeatures(unittest.TestCase):
particleContainer: ParticleContainer = get_default_ParticleContainer()
features: np.ndarray = cmeth.get_particle_featurematrix(particleContainer)
self.assertEqual(features.shape[0], len(particleContainer.particles))
for entry in features[0, :]:
self.assertTrue(type(entry) in [float, np.float64])
class TestTrainedSubsampling(unittest.TestCase):
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment