Pythonを用いた粒子解析の基礎知識と応用例21選

import matplotlib.pyplot as plt
from skimage import io

# 画像の読み込み
image = io.imread('particle_image.jpg')

# 画像の表示
plt.imshow(image)
plt.axis('off')  # 軸を非表示に
plt.show()

from skimage import filters

# ガウシアンフィルタの適用
denoised = filters.gaussian(image, sigma=1)

# 結果の表示
plt.imshow(denoised)
plt.axis('off')
plt.show()

from skimage import filters

# 大津の二値化を適用
threshold = filters.threshold_otsu(denoised)
binary = denoised > threshold

# 結果の表示
plt.imshow(binary, cmap='gray')
plt.axis('off')
plt.show()

from skimage import morphology

# クロージング処理
closed = morphology.closing(binary, morphology.square(3))

# 結果の表示
plt.imshow(closed, cmap='gray')
plt.axis('off')
plt.show()

from skimage import measure

# ラベリング処理
labeled = measure.label(closed)

# 結果の表示
plt.imshow(labeled, cmap='nipy_spectral')
plt.axis('off')
plt.show()

import numpy as np
from skimage import measure

# ラベリング済みの画像を使用
props = measure.regionprops(labeled_image)

# 面積と周囲長の計算
areas = [prop.area for prop in props]
perimeters = [prop.perimeter for prop in props]

print(f"平均面積: {np.mean(areas):.2f}")
print(f"平均周囲長: {np.mean(perimeters):.2f}")

# アスペクト比と円形度の計算
aspect_ratios = [prop.major_axis_length / prop.minor_axis_length for prop in props]
circularities = [4 * np.pi * prop.area / (prop.perimeter ** 2) for prop in props]

print(f"平均アスペクト比: {np.mean(aspect_ratios):.2f}")
print(f"平均円形度: {np.mean(circularities):.2f}")

# 重心座標の計算
centroids = [prop.centroid for prop in props]

print("重心座標の例:")
for i, centroid in enumerate(centroids[:5]):
    print(f"粒子{i+1}: (y={centroid[0]:.2f}, x={centroid[1]:.2f})")

# 向きと主軸の計算
orientations = [prop.orientation for prop in props]
major_axis_lengths = [prop.major_axis_length for prop in props]

print(f"平均向き (ラジアン): {np.mean(orientations):.2f}")
print(f"平均主軸長: {np.mean(major_axis_lengths):.2f}")

from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

# 特徴量の準備
features = np.column_stack((areas, perimeters, aspect_ratios, circularities))

# 特徴量の標準化
scaler = StandardScaler()
features_scaled = scaler.fit_transform(features)

# K-meansクラスタリングの実行
kmeans = KMeans(n_clusters=3, random_state=42)
clusters = kmeans.fit_predict(features_scaled)

print("クラスタ別の粒子数:")
for i in range(3):
    print(f"クラスタ{i+1}: {np.sum(clusters == i)}個")

from scipy.optimize import linear_sum_assignment

def track_particles(centroids_t1, centroids_t2, max_distance=50):
    distances = np.sqrt(((centroids_t1[:, np.newaxis] - centroids_t2) ** 2).sum(axis=2))
    row_ind, col_ind = linear_sum_assignment(distances)
    valid_tracks = distances[row_ind, col_ind] < max_distance
    return row_ind[valid_tracks], col_ind[valid_tracks]

# 2つの時点での重心座標（例として）
centroids_t1 = np.array(centroids[:10])
centroids_t2 = centroids_t1 + np.random.randn(*centroids_t1.shape) * 5

tracked_indices = track_particles(centroids_t1, centroids_t2)

print(f"追跡された粒子数: {len(tracked_indices[0])}")

from skimage import io, measure

# 3D画像の読み込み（例として）
image_3d = io.imread('3d_particle_image.tif')

# 3Dラベリング
labels_3d = measure.label(image_3d)

# 3D特徴量の抽出
props_3d = measure.regionprops(labels_3d)

print(f"検出された3D粒子数: {len(props_3d)}")
print(f"最大体積の粒子: {max(prop.volume for prop in props_3d)}")

from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score

# 特徴量とラベルの準備（例として）
X = np.column_stack((areas, perimeters, aspect_ratios, circularities))
y = np.random.choice(['A', 'B', 'C'], size=len(X))  # ダミーラベル

# データの分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# ランダムフォレストモデルの構築と学習
clf = RandomForestClassifier(n_estimators=100, random_state=42)
clf.fit(X_train, y_train)

# テストデータでの予測と精度評価
y_pred = clf.predict(X_test)
accuracy = accuracy_score(y_test, y_pred)

print(f"分類精度: {accuracy:.2f}")

import matplotlib.pyplot as plt
import seaborn as sns

# ヒストグラムの作成
plt.figure(figsize=(10, 6))
sns.histplot(areas, bins=30, kde=True)
plt.title('粒子サイズ分布')
plt.xlabel('面積')
plt.ylabel('頻度')
plt.show()

# 散布図の作成
plt.figure(figsize=(10, 10))
plt.scatter([c[1] for c in centroids], [c[0] for c in centroids], alpha=0.5)
plt.title('粒子の空間分布')
plt.xlabel('X座標')
plt.ylabel('Y座標')
plt.gca().invert_yaxis()  # Y軸を反転（画像座標系に合わせる）
plt.show()

# 箱ひげ図の作成
plt.figure(figsize=(12, 6))
data = [areas, perimeters, aspect_ratios, circularities]
labels = ['面積', '周囲長', 'アスペクト比', '円形度']
plt.boxplot(data, labels=labels)
plt.title('粒子特徴量の分布')
plt.ylabel('値')
plt.show()

import numpy as np

# 星座図の作成
def radar_chart(data, labels, title):
    angles = np.linspace(0, 2*np.pi, len(labels), endpoint=False)
    data = np.concatenate((data, [data[0]]))  # 閉じた図形にするため
    angles = np.concatenate((angles, [angles[0]]))  # 同上

    fig, ax = plt.subplots(figsize=(6, 6), subplot_kw=dict(projection='polar'))
    ax.plot(angles, data)
    ax.fill(angles, data, alpha=0.25)
    ax.set_xticks(angles[:-1])
    ax.set_xticklabels(labels)
    ax.set_title(title)
    plt.show()

# 例として最初の粒子のデータを使用
sample_data = [props[0].area, props[0].perimeter, props[0].major_axis_length / props[0].minor_axis_length, 
               4 * np.pi * props[0].area / (props[0].perimeter ** 2)]
sample_labels = ['面積', '周囲長', 'アスペクト比', '円形度']

radar_chart(sample_data, sample_labels, '粒子形状の特徴')

from skimage import exposure

# コントラスト強調
def enhance_contrast(image):
    p2, p98 = np.percentile(image, (2, 98))
    return exposure.rescale_intensity(image, in_range=(p2, p98))

# コントラスト強調を適用
enhanced_image = enhance_contrast(image)

# 結果の表示
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
ax1.imshow(image, cmap='gray')
ax1.set_title('原画像')
ax2.imshow(enhanced_image, cmap='gray')
ax2.set_title('コントラスト強調後')
plt.show()

from scipy import ndimage as ndi
from skimage.segmentation import watershed
from skimage.feature import peak_local_max

def separate_overlapping_particles(binary_image):
    # 距離変換
    distance = ndi.distance_transform_edt(binary_image)

    # ローカルマキシマムの検出
    coords = peak_local_max(distance, footprint=np.ones((3, 3)), labels=binary_image)
    mask = np.zeros(distance.shape, dtype=bool)
    mask[tuple(coords.T)] = True
    markers, _ = ndi.label(mask)

    # ウォーターシェッド法による分離
    labels = watershed(-distance, markers, mask=binary_image)

    return labels

# オーバーラップ粒子の分離を適用
separated_particles = separate_overlapping_particles(binary)

# 結果の表示
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
ax1.imshow(binary, cmap='gray')
ax1.set_title('二値化画像')
ax2.imshow(separated_particles, cmap='nipy_spectral')
ax2.set_title('分離後')
plt.show()

def process_large_image(image_path, chunk_size=1000):
    # 画像を小さなチャンクに分割して処理
    with rasterio.open(image_path) as src:
        height = src.height
        width = src.width

        results = []

        for y in range(0, height, chunk_size):
            for x in range(0, width, chunk_size):
                window = rasterio.windows.Window(x, y, 
                                                 min(chunk_size, width - x),
                                                 min(chunk_size, height - y))
                chunk = src.read(1, window=window)

                # チャンクごとの処理（例：平均値の計算）
                chunk_result = np.mean(chunk)
                results.append(chunk_result)

        return np.mean(results)

# 大規模画像の処理
average_value = process_large_image('large_image.tif')
print(f"画像全体の平均値: {average_value}")

import cv2
import numpy as np
from skimage import measure

def count_cells(image_path):
    # 画像の読み込みとグレースケール変換
    image = cv2.imread(image_path)
    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

    # 二値化
    _, binary = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)

    # ノイズ除去
    kernel = np.ones((3,3), np.uint8)
    opening = cv2.morphologyEx(binary, cv2.MORPH_OPEN, kernel, iterations=2)

    # 細胞の検出
    labels = measure.label(opening)
    cell_count = labels.max()

    return cell_count

# 細胞計数の実行
image_path = 'cell_image.jpg'
cell_count = count_cells(image_path)
print(f"検出された細胞数: {cell_count}")

import numpy as np
import matplotlib.pyplot as plt
from skimage import measure, io

def analyze_grain_size(image_path):
    # 画像の読み込みと二値化
    image = io.imread(image_path, as_gray=True)
    binary = image > threshold_otsu(image)

    # ラベリングと領域プロパティの取得
    labels = measure.label(binary)
    props = measure.regionprops(labels)

    # 粒子サイズの取得
    sizes = [prop.area for prop in props]

    # 粒度分布のプロット
    plt.hist(sizes, bins=20)
    plt.title('粒度分布')
    plt.xlabel('粒子サイズ')
    plt.ylabel('頻度')
    plt.show()

    return np.mean(sizes), np.std(sizes)

# 粒度分布解析の実行
image_path = 'material_image.jpg'
mean_size, std_size = analyze_grain_size(image_path)
print(f"平均粒子サイズ: {mean_size:.2f}")
print(f"粒子サイズの標準偏差: {std_size:.2f}")

import cv2
import numpy as np
from skimage import measure, color

def detect_microplastics(image_path):
    # 画像の読み込み
    image = cv2.imread(image_path)
    hsv = cv2.cvtColor(image, cv2.COLOR_BGR2HSV)

    # 色相に基づくマイクロプラスチックの検出
    lower_bound = np.array([0, 50, 50])
    upper_bound = np.array([10, 255, 255])
    mask = cv2.inRange(hsv, lower_bound, upper_bound)

    # ノイズ除去
    kernel = np.ones((5,5), np.uint8)
    opening = cv2.morphologyEx(mask, cv2.MORPH_OPEN, kernel, iterations=2)

    # マイクロプラスチックの計数と面積計算
    labels = measure.label(opening)
    props = measure.regionprops(labels)

    count = len(props)
    total_area = sum(prop.area for prop in props)

    return count, total_area

# マイクロプラスチック検出の実行
image_path = 'water_sample.jpg'
count, area = detect_microplastics(image_path)
print(f"検出されたマイクロプラスチック数: {count}")
print(f"マイクロプラスチックの総面積: {area} ピクセル")

import cv2
import numpy as np
from skimage import measure

def quality_control(image_path, defect_threshold):
    # 画像の読み込みとグレースケール変換
    image = cv2.imread(image_path)
    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

    # エッジ検出
    edges = cv2.Canny(gray, 50, 150)

    # 欠陥の検出
    contours, _ = cv2.findContours(edges, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)

    defects = []
    for contour in contours:
        area = cv2.contourArea(contour)
        if area > defect_threshold:
            defects.append(contour)

    # 結果の可視化
    result = image.copy()
    cv2.drawContours(result, defects, -1, (0, 0, 255), 2)

    cv2.imshow('Defect Detection', result)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

    return len(defects)

# 品質管理の実行
image_path = 'product_image.jpg'
defect_count = quality_control(image_path, defect_threshold=100)
print(f"検出された欠陥数: {defect_count}")