旅行好きなソフトエンジニアの備忘録

プログラミングや技術関連のメモを始めました

【異常検知】オートエンコーダーを用いた画像の異常検知

以下のサイトで画像の異常検知をやっていて面白そうなので自分でも試してみました。 qiita.com

--- 試した環境 ---
Windows10
Python 3.6
Keras 2.1.4
Tensorflow-gpu 1.5.0

使うデータセットは9クラスに分類されたキュウリの画像です。 github.com

以下の写真のように9クラスに分類されていて、最高品質と思われる2Lのみを学習させ、最低品質と思われるCを検出できるかどうか試します。 f:id:ni4muraano:20180708172007j:plain

まず、以下のようにautoencoder.pyにオートエンコーダーを定義します。画像に適用するので、Convolutionalオートエンコーダーを定義しています。

#!/usr/bin/env python
# -*- coding: utf-8 -*-
from keras.models import Model
from keras.layers import Input
from keras.layers.convolutional import Conv2D, MaxPooling2D, Conv2DTranspose

class AutoEncoder(object):
    def __init__(self, input_shape, first_layer_channels):
        self.CONV_FILTER_SIZE = 3
        self.CONV_STRIDE = 1
        self.DECONV_FILTER_SIZE = 3
        self.DECONV_STRIDE = 2

        # (32 x 32 x 3)
        inputs = Input(input_shape)

        # エンコーダーの作成
        # (16 x 16 x N)
        filter_count = first_layer_channels
        enc1 = self._add_encoding_layer(filter_count, inputs)

        # (8 x 8 x 2N)
        filter_count = first_layer_channels*2
        enc2 = self._add_encoding_layer(filter_count, enc1)

        # (4 x 4 x 4N)
        filter_count = first_layer_channels*4
        enc3 = self._add_encoding_layer(filter_count, enc2)

        # (8 x 8 x 2N)
        filter_count = first_layer_channels*2
        dec3 = self._add_decoding_layer(filter_count, enc3)

        # (16 x 16 x N)
        filter_count = first_layer_channels
        dec2 = self._add_decoding_layer(filter_count, dec3)

        # (32 x 32 x 3)
        filter_count = input_shape[2]
        dec1 = self._add_decoding_layer(filter_count, dec2)

        self.AutoEncoder = Model(input=inputs, output=dec1)
        #print(self.AutoEncoder.summary())

    def _add_encoding_layer(self, filter_count, sequence):
        new_sequence = Conv2D(filter_count, self.CONV_FILTER_SIZE, strides=self.CONV_STRIDE, padding='same', activation='relu')(sequence)
        new_sequence = MaxPooling2D()(new_sequence)
        return new_sequence

    def _add_decoding_layer(self, filter_count, sequence):
        new_sequence = Conv2DTranspose(filter_count, self.DECONV_FILTER_SIZE, strides=self.DECONV_STRIDE, padding='same',
                                       kernel_initializer='he_uniform', activation='relu')(sequence)
        return new_sequence

    def get_model(self):
        return self.AutoEncoder


次にtrain.pyを作成します。train.pyには画像読み込み処理、学習処理を書きます。

import numpy as np
from sklearn.model_selection import train_test_split
from autoencoder import AutoEncoder
from keras.optimizers import Adam, RMSprop, SGD

def load_image_and_label(pickled_files):
    # Each file contains 495 images
    IMAGE_COUNT_PER_FILE = 495
    # Image shape is 32x32x3
    ROW = 32
    COL = 32
    DIM = 3
    whole_images = np.empty((IMAGE_COUNT_PER_FILE*len(pickled_files), ROW, COL, DIM))
    whole_labels = np.empty(IMAGE_COUNT_PER_FILE*len(pickled_files))
    for i, pickled_file in enumerate(pickled_files):
        dict = _unpickle(pickled_file)
        images = dict['data'].reshape(IMAGE_COUNT_PER_FILE, DIM, ROW, COL).transpose(0, 2, 3, 1)
        whole_images[i*IMAGE_COUNT_PER_FILE:(i + 1)*IMAGE_COUNT_PER_FILE, :, :, :] = images
        labels = dict['labels']
        whole_labels[i*IMAGE_COUNT_PER_FILE:(i + 1)*IMAGE_COUNT_PER_FILE] = labels
    return (whole_images, whole_labels)

def _unpickle(pickled_file):
    import pickle

    with open(pickled_file, 'rb') as file:
        # You'll have an error without "encoding='latin1'"
        dict = pickle.load(file, encoding='latin1')
    return dict

# Function to load cucumber-9 dataset and split it into training and test data
def load():
    (X1, y1) = load_image_and_label(['Train\\data_batch_1',
                                     'Train\\data_batch_2',
                                     'Train\\data_batch_3',
                                     'Train\\data_batch_4',
                                     'Train\\data_batch_5'])
    (X2, y2) = load_image_and_label(['Test\\test_batch'])
    X = np.concatenate((X1, X2), axis=0)
    y = np.concatenate((y1, y2), axis=0)
    # 2L as normal
    normal_index = np.where((y == 0))
    # C as anomaly
    anomaly_index = np.where(y == 8)
    X_normal = X[normal_index]
    X_anomaly = X[anomaly_index]
    y_normal = y[normal_index]
    y_anomaly = y[anomaly_index]
    # split normal images into train and test data
    X_train, X_test, y_train, y_test = train_test_split(X_normal, y_normal, test_size=0.2, stratify=y_normal, random_state=0)
    X_test = np.concatenate((X_test, X_anomaly), axis=0)
    y_test = np.concatenate((y_test, y_anomaly), axis=0)
    y_test = y_test == 8
    return X_train/255, X_test/255, y_test

if __name__ == '__main__':
    batch_size = 32
    epochs = 500
    first_layer_channels = 32

    X_train, X_test, y_test = load()
    input_shape = X_train[0].shape
    auto_encoder = AutoEncoder(input_shape, first_layer_channels)
    model = auto_encoder.get_model()
    model.compile(optimizer=Adam(lr=1e-3, amsgrad=True), loss='mse')
    model.fit(X_train, X_train, batch_size=batch_size, epochs=epochs)
    model.save_weights('weights.hdf5')


最後に予測処理部をtest.pyとして書きます。

from train_cucumber1 import load
from autoencoder import AutoEncoder
import numpy as np

if __name__ == '__main__':
    batch_size = 64
    first_layer_channels = 32

    # load data
    X_train, X_test, y_test = load()
    # load model
    input_shape = X_train[0].shape
    auto_encoder = AutoEncoder(input_shape, first_layer_channels)
    model = auto_encoder.get_model()
    model.load_weights('weights.hdf5')
    # Apply prediction to training data to determine threshold
    y_train_predict = model.predict(X_train, batch_size=batch_size)
    # Check difference between inputs and outputs
    diff_train = np.sum(np.abs(X_train - y_train_predict), axis=(1,2,3))
    diff_train_mean = np.mean(diff_train)
    diff_train_std = np.std(diff_train)
    #print(str(diff_train_mean))
    #print(str(diff_train_std))
    # Determine threshold
    thresh = diff_train_mean + 2*diff_train_std

    # Apply prediction to test data
    y_test_predict = model.predict(X_test, batch_size=batch_size)
    diff_test = np.sum(np.abs(X_test - y_test_predict), axis=(1,2,3))
    # First 66 data is 2L quality so these values are expected to have small values
    diff_test_normal = diff_test[:66]
    # Remaining data is C quality so these values are expected to have large values
    diff_test_anomaly = diff_test[66:]
    true_positive = np.sum(diff_test_anomaly >= thresh)
    false_positive = np.sum(diff_test_normal >= thresh)
    false_negative = np.sum(diff_test_anomaly < thresh)
    precision = true_positive/(true_positive + false_positive)
    recall = true_positive/(true_positive + false_negative)
    f_measure = 2*recall*precision/(recall + precision)
    print(f_measure)


正常画像(2L)に対する復元誤差(青色)と異常画像(C)に対する復元誤差(橙色)の分布を以下に描画します。F値は0.96と良好な値ではあるものの、分布を見ると結構被ってる印象を受けるため、実戦で使うには一工夫必要そうです。 f:id:ni4muraano:20180708190744p:plain

【Python】 x-meansを使ってクラスタリングする例

k-meansのクラスタ数を自動で決めてくれるx-meansという手法があることを以下で知ったので試してみようとしました。 qiita.com ただリンク先のコードは自分のデータセットに対しては「行列式が0で計算できない」といったようなエラーが出ました。それで代わりにx-meansを実装しているpyclusteringというライブラリで試したので、pyclusteringのx-meansの使い方をメモとして残します(pyclusteringはpip install pyclusteringでインストールできます)。

# -*- coding: utf-8 -*-
import numpy as np
from pyclustering.cluster.xmeans import xmeans, kmeans_plusplus_initializer
from pyclustering.utils import draw_clusters

if __name__ == "__main__":
    # データの準備
    x = np.array([np.random.normal(loc, 0.1, 20) for loc in np.repeat([1,2], 2)]).flatten() #ランダムな80個の数を生成
    y = np.array([np.random.normal(loc, 0.1, 20) for loc in np.tile([1,2], 2)]).flatten() #ランダムな80個の数を生成
    points = np.concatenate((x[np.newaxis,:], y[np.newaxis,:]), axis=0).T.tolist()

    # クラスタ数2から探索させてみる
    initial_centers = kmeans_plusplus_initializer(points, 2).initialize()
    # クラスタリングの実行
    instances = xmeans(points, initial_centers, ccore=True)
    instances.process()
    # クラスタはget_clustersで取得できる
    clusters = instances.get_clusters()
    # クラスタの平均は以下で計算できる
    centers = []
    points = np.array(points)
    for cluster in clusters:
        c = points[cluster, :]
        centers.append(np.mean(c, axis=0))
    # クラスタリング結果を表示する
    draw_clusters(points, clusters)

f:id:ni4muraano:20180617101517p:plain

【C#】byte*をbyte[]に変換する

C++のDLLからJPEGエンコードされた画像データ(型はunsinged charのポインタ)が送られてくるので、それをC#側でbyteポインタで受け取りbyte[]に変換後、OpenCvSharpのMatに変換したかった話です。byteポインタをbyte配列に変換するにはMarshal.Copyを利用すると簡単とのことでした。

// byte*のサイズはlenと知っているものとする
byte[] arr = new byte[len];
// ptrがbyte*
Marshal.Copy((IntPtr)ptr, arr, 0, len);


stackoverflow.com

【C#】エンコードされたbyte配列をMatに変換する

C++のDLLからJPEGエンコードされた画像データ(型はunsinged charのポインタ)が送られてくるので、それをC#側でbyteポインタで受け取りbyte[]に変換後、OpenCvSharpのMatに変換したかった話です。調べてみるとOpenCvSharpにはImDecodeというメソッドがあり、それを使えば良いだけでした。

// bytesの型はbyte[]
Mat mat = Cv2.ImDecode(bytes, ImreadModes.Color); 

github.com

【統計】母欠点数に関する検定

ある製品を検査すると従来N個の欠点が見つかっていた。この対策が行われ、効果を検証するため製品をM個ピックアップしたところ、L個の欠点が見つかった。さて、欠点は減ったと言えるかどうか、というお題があって書籍「入門統計解析法」で調べると226ページからがまさにその内容だったので判定部のコードをメモします。

static void Main(string[] args)
{
    // 従来製品毎に3個の欠点が見つかっていた
    int conventionalDefectNumberPerProduct = 3;
    // 工程工夫後の製品を10個ピックアップしたら18個欠点が見つかった
    int currentDefectNumber = 18;
    int productNumber = 10;
    // 欠点は減った?
    Console.WriteLine(DoesDefectNumberDecrease(conventionalDefectNumberPerProduct, currentDefectNumber, productNumber));
}

static bool DoesDefectNumberDecrease(int conventionalDefectNumberPerProduct, int currentDefectNumber, int productNumber)
{
    // 有意水準0.05の時の臨界値
    const double criticalValue = -1.645;

    double rambda = (currentDefectNumber + 0.5) / productNumber;
    double u0 = (rambda - conventionalDefectNumberPerProduct) / Math.Sqrt((double)conventionalDefectNumberPerProduct / productNumber);
    return u0 <= criticalValue;
}

入門 統計解析法

入門 統計解析法

【C#】TimeSpanに割り算を適用する方法

複数のTimeSpanの平均値を求めたかったため、TimeSpanに割り算を適用する方法を調べたのですが、以下のやり方で出来ることが分かりました。

// 60秒
var t1 = new TimeSpan(0, 0, 60);
// Ticksに割り算を適用してTimeSpanのコンストラクタに入れる
var t2 = new TimeSpan(t1.Ticks/10);
// 6と表示される
Console.WriteLine(t2.Seconds);

情報源は以下。 stackoverflow.com

【C#】動的にChartを追加する

円グラフを描きたいけれども、アプリを動作させてからでないと何個の円グラフを描くか決めれない状況のため、動的にChartを追加する方法を調べました。グラフを描くために使ったライブラリはLiveCharts.Wpfになります。

まずはxamlです。以下で定義したStackPanel1に動的にChartを追加します。

<Window x:Class="LiveChartsExample.MainWindow"
        xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
        xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"
        xmlns:d="http://schemas.microsoft.com/expression/blend/2008"
        xmlns:mc="http://schemas.openxmlformats.org/markup-compatibility/2006"
        xmlns:local="clr-namespace:LiveChartsExample"
        xmlns:lvc="clr-namespace:LiveCharts.Wpf;assembly=LiveCharts.Wpf"
        mc:Ignorable="d"
        Title="MainWindow" Height="350" Width="525">
    <StackPanel Name="StackPanel1" Orientation="Horizontal">
    </StackPanel>
</Window>


次にC#側のソースです。

using System;
using System.Windows;
using System.Windows.Controls;
using System.Windows.Media;

using LiveCharts;
using LiveCharts.Wpf;

namespace LiveChartsExample
{
    /// <summary>
    /// MainWindow.xaml の相互作用ロジック
    /// </summary>
    public partial class MainWindow : Window
    {
        public MainWindow()
        {
            InitializeComponent();

            var random = new Random();
            // 例として円グラフを4つ描いてみる
            for (int i = 0; i < 4; ++i)
            {
                // グリッドの定義
                var grid = new Grid();
                grid.Width = 150;
                var row1 = new RowDefinition();
                row1.Height = new GridLength(30);
                var row2 = new RowDefinition();
                row2.Height = new GridLength(150);
                grid.RowDefinitions.Add(row1);
                grid.RowDefinitions.Add(row2);

                // 表題用TextBlock
                var title = new TextBlock();
                title.Text = "Chart" + i.ToString();
                title.FontSize = 20;
                title.Foreground = new SolidColorBrush(Colors.Black);
                title.HorizontalAlignment = HorizontalAlignment.Center;

                // 円グラフの定義
                var chart = new PieChart();
                chart.Name = "Chart" + i.ToString();
                chart.StartingRotationAngle = 0;
                chart.Width = 150;
                chart.Height = 150;
                chart.MouseDown += Chart_MouseDown;
                chart.LegendLocation = LegendLocation.Bottom;

                var series1 = new PieSeries()
                {
                    Title = "A",
                    Values = new ChartValues<int> { random.Next(100) },
                    DataLabels = true,
                    LabelPoint = point => string.Format("{0} ({1:P})", point.Y, point.Participation)
                };
                var series2 = new PieSeries()
                {
                    Title = "B",
                    Values = new ChartValues<int> { random.Next(100) },
                    DataLabels = true,
                    LabelPoint = point => string.Format("{0} ({1:P})", point.Y, point.Participation)
                };
                chart.Series.Add(series1);
                chart.Series.Add(series2);

                // グリッドの0行目に表題、1行目にグラフを設置
                Grid.SetRow(title, 0);
                Grid.SetRow(chart, 1);
                grid.Children.Add(title);
                grid.Children.Add(chart);

                // スタックパネルにグリッドを追加
                StackPanel1.Children.Add(grid);
            }
        }

        private void Chart_MouseDown(object sender, System.Windows.Input.MouseButtonEventArgs e)
        {
            PieChart chart = sender as PieChart;

            // 以下チャートがクリックされた時の動作を書く
        }

    }
}

上記を実行すると以下のように円グラフが4つ表示されます。 f:id:ni4muraano:20180513121255p:plain
それにしてもLiveChartsはグラフに表題を付けるのにこんな書き方をしないとダメなのでしょうか。PieChartクラスにTitleというプロパティが無く、サイトの例を見ても分からなかったので今回はGridにTextBlockとChartを置くことでTextBlockがChartのタイトルに見えるようにしています。