{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom sklearn.decomposition import PCA\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.cluster import KMeans\nfrom keras.utils import Sequence\nfrom keras import Model\nfrom keras.optimizers import Adam\nfrom keras.losses import mean_absolute_error, mean_squared_error","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Dataset\n\nI use my [dataset of random lung 64x64 slices](https://www.kaggle.com/miklgr500/osic-random-slices-from-lung-regions), lung slice help autoencoder detect lung features and realise attention idea for autoencoder. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"paths = []\nfor root, subdirs, files in os.walk('../input/osic-random-slices-from-lung-regions/image_slice'):\n    for file in files:\n        if os.path.splitext(file)[1].lower() in ('.npz'):\n             paths.append(os.path.join(root, file))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"class ImageGenerator(Sequence):\n    def __init__(self, paths, batch_size=64, random_state=None):\n        self._paths = paths\n        self._batch_size = batch_size\n        self._random = np.random.RandomState(random_state)\n        \n    def __len__(self):\n        return len(self._paths) // self._batch_size\n    \n    def __getitem__(self, idx):\n        paths = self._random.choice(self._paths, size=self._batch_size)\n        imgs = np.expand_dims([np.load(p, allow_pickle=True)['arr_0'].tolist()['slice'] for p in paths], axis=-1) / 1_000\n        return imgs, imgs","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_gen = ImageGenerator(paths, random_state=43)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"batch = data_gen[0]\n\nfig, ax = plt.subplots(8, 8, figsize=(16, 16))\n\nfor i, _ax in enumerate(ax):\n    for j, __ax in enumerate(_ax):\n        __ax.imshow(batch[0][8*i + j, ..., 0])\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## AutoEncoder\n![](https://upload.wikimedia.org/wikipedia/commons/2/28/Autoencoder_structure.png)\nPipeline released in this kernel is very simple:\n\n* train encoder and used one for translate image in vector\n* visulize results vectors used PCA decomposition algorithm","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras.layers import (\n    Dense, Dropout, Activation, Flatten, Input, BatchNormalization, UpSampling2D, Add, Conv2D, AveragePooling2D, LeakyReLU\n)\n\nfrom tensorflow.keras import Model\nfrom tensorflow.keras.optimizers import Nadam\n\ndef get_encoder(shape=(512, 512, 1)):\n    def res_block(x, n_features):\n        _x = x\n        x = BatchNormalization()(x)\n        x = LeakyReLU()(x)\n    \n        x = Conv2D(n_features, kernel_size=(3, 3), strides=(1, 1), padding='same')(x)\n        x = Add()([_x, x])\n        return x\n    \n    inp = Input(shape=shape)\n    \n    # 64\n    x = Conv2D(64, kernel_size=(3, 3), strides=(1, 1), padding='same')(inp)\n    x = BatchNormalization()(x)\n    x = LeakyReLU()(x)\n    \n    x = Conv2D(64, kernel_size=(3, 3), strides=(1, 1), padding='same')(x)\n    x = BatchNormalization()(x)\n    x = LeakyReLU()(x)\n    \n    x = AveragePooling2D(pool_size=(2, 2), strides=(2, 2))(x)\n    \n    # 32\n    x = Conv2D(128, kernel_size=(3, 3), strides=(1, 1), padding='same')(x)\n    for _ in range(3):\n        x = res_block(x, 128)\n    x = AveragePooling2D(pool_size=(2, 2), strides=(2, 2))(x)\n    \n    # 16\n    x = Conv2D(256, kernel_size=(3, 3), strides=(1, 1), padding='same')(x)\n    for _ in range(3):\n        x = res_block(x, 256)\n    x = AveragePooling2D(pool_size=(2, 2), strides=(2, 2))(x)\n    \n    # 64\n    x = Conv2D(128, kernel_size=(3, 3), strides=(1, 1), padding='same')(x)\n    for _ in range(3):\n        x = res_block(x, 128)\n    x = AveragePooling2D(pool_size=(2, 2), strides=(2, 2))(x)\n    \n    # 8\n    x = Conv2D(32, kernel_size=(3, 3), strides=(1, 1), padding='same')(x)\n    for _ in range(2):\n        x = res_block(x, 32)    \n    \n    # 8\n    x = Conv2D(1, kernel_size=(1, 1), strides=(1, 1), padding='same')(x)\n    return Model(inp, x)\n\n\n\ndef get_decoder(shape=(8, 8, 1)):\n    inp = Input(shape=shape)\n\n    # 8\n    x = UpSampling2D((2, 2))(inp)\n    x = Conv2D(256, kernel_size=(3, 3), strides=(1, 1), padding='same')(x)\n    x = BatchNormalization()(x)\n    x = LeakyReLU()(x)\n    \n    # 16\n    x = UpSampling2D((2, 2))(x)\n    x = Conv2D(128, kernel_size=(3, 3), strides=(1, 1), padding='same')(x)\n    x = BatchNormalization()(x)\n    x = LeakyReLU()(x)\n    \n    # 32\n    x = UpSampling2D((2, 2))(x)\n    x = Conv2D(64, kernel_size=(3, 3), strides=(1, 1), padding='same')(x)\n    x = BatchNormalization()(x)\n    x = LeakyReLU()(x)\n\n    x = UpSampling2D((2, 2))(x)\n    x = Conv2D(32, kernel_size=(3, 3), strides=(1, 1), padding='same')(x)\n    x = BatchNormalization()(x)\n    x = LeakyReLU()(x)\n    \n    x = Conv2D(1, kernel_size=(1, 1), strides=(1, 1), padding='same')(x)\n    return Model(inp, x)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"encoder = get_encoder((64, 64, 1))\ndecoder = get_decoder((4, 4, 1))\n\ninp = Input((64, 64, 1))\ne = encoder(inp)\nd = decoder(e)\nmodel = Model(inp, d)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"encoder.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"decoder.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.compile(optimizer=Adam(lr=3*1e-3), loss='mae')\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_gen = ImageGenerator(paths, random_state=43, batch_size=512)\nmodel.fit_generator(data_gen, steps_per_epoch=len(data_gen), epochs=5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## PCA decomposition","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"x_encode = encoder.predict_generator(data_gen, verbose=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"e_samples, x_samples = [], []\nfor i in range(25):\n    x, _ = data_gen[i]\n    _x = encoder.predict(x)\n    e_samples.extend(_x.tolist())\n    x_samples.extend(x.tolist())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"e_stack, x_stack = np.array(e_samples), np.array(x_samples)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"x_encode = x_encode.reshape((-1, 16))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scaler = StandardScaler()\nx_scale = scaler.fit_transform(x_encode)\n\npca = PCA()\nx_pca = pca.fit_transform(x_scale)\nsample_pca = pca.transform(scaler.transform(e_stack.reshape((-1, 16))))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"kmeans = KMeans(n_clusters=10)\nx_c = kmeans.fit_predict(x_pca)\nsample_c = kmeans.predict(sample_pca.astype(x_pca.dtype))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,8))\nplt.plot(1 - pca.explained_variance_ratio_);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,8))\nplt.plot(x_pca[:, 0], x_pca[:, 1], '.', alpha=0.1);\nplt.plot(sample_pca[:, 0], sample_pca[:, 1], '.', alpha=0.1);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import random\ncolor = [\"#\"+''.join([random.choice('0123456789ABCDEF') for j in range(6)]) for i in range(10)]\n\nplt.figure(figsize=(10,8))\n\nfor c in range(10):\n    subset = x_pca[x_c==c]\n    plt.plot(subset[:, 0], subset[:, 1], '.', alpha=0.1, color=color[c]);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Cluster 0","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"x = x_stack[sample_c == 0]\n\nidx = np.random.randint(0, len(x), size=64)\nbatch = x[idx]\n\nfig, ax = plt.subplots(8, 8, figsize=(16, 16))\n\nfor i, _ax in enumerate(ax):\n    for j, __ax in enumerate(_ax):\n        __ax.imshow(batch[8*i + j, ..., 0])\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Cluster 1","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"x = x_stack[sample_c == 1]\n\nidx = np.random.randint(0, len(x), size=64)\nbatch = x[idx]\n\nfig, ax = plt.subplots(8, 8, figsize=(16, 16))\n\nfor i, _ax in enumerate(ax):\n    for j, __ax in enumerate(_ax):\n        __ax.imshow(batch[8*i + j, ..., 0])\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Cluster 2","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"x = x_stack[sample_c == 2]\n\nidx = np.random.randint(0, len(x), size=64)\nbatch = x[idx]\n\nfig, ax = plt.subplots(8, 8, figsize=(16, 16))\n\nfor i, _ax in enumerate(ax):\n    for j, __ax in enumerate(_ax):\n        __ax.imshow(batch[8*i + j, ..., 0])\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Cluster 3","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"x = x_stack[sample_c == 3]\n\nidx = np.random.randint(0, len(x), size=64)\nbatch = x[idx]\n\nfig, ax = plt.subplots(8, 8, figsize=(16, 16))\n\nfor i, _ax in enumerate(ax):\n    for j, __ax in enumerate(_ax):\n        __ax.imshow(batch[8*i + j, ..., 0])\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Cluster 4","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"x = x_stack[sample_c == 4]\n\nidx = np.random.randint(0, len(x), size=64)\nbatch = x[idx]\n\nfig, ax = plt.subplots(8, 8, figsize=(16, 16))\n\nfor i, _ax in enumerate(ax):\n    for j, __ax in enumerate(_ax):\n        __ax.imshow(batch[8*i + j, ..., 0])\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Cluster 5","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"x = x_stack[sample_c == 5]\n\nidx = np.random.randint(0, len(x), size=64)\nbatch = x[idx]\n\nfig, ax = plt.subplots(8, 8, figsize=(16, 16))\n\nfor i, _ax in enumerate(ax):\n    for j, __ax in enumerate(_ax):\n        __ax.imshow(batch[8*i + j, ..., 0])\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Cluster 6","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"x = x_stack[sample_c == 6]\n\nidx = np.random.randint(0, len(x), size=64)\nbatch = x[idx]\n\nfig, ax = plt.subplots(8, 8, figsize=(16, 16))\n\nfor i, _ax in enumerate(ax):\n    for j, __ax in enumerate(_ax):\n        __ax.imshow(batch[8*i + j, ..., 0])\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Cluster 7","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"x = x_stack[sample_c == 7]\n\nidx = np.random.randint(0, len(x), size=64)\nbatch = x[idx]\n\nfig, ax = plt.subplots(8, 8, figsize=(16, 16))\n\nfor i, _ax in enumerate(ax):\n    for j, __ax in enumerate(_ax):\n        __ax.imshow(batch[8*i + j, ..., 0])\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Cluster 8","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"x = x_stack[sample_c == 8]\n\nidx = np.random.randint(0, len(x), size=64)\nbatch = x[idx]\n\nfig, ax = plt.subplots(8, 8, figsize=(16, 16))\n\nfor i, _ax in enumerate(ax):\n    for j, __ax in enumerate(_ax):\n        __ax.imshow(batch[8*i + j, ..., 0])\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Cluster 9","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"x = x_stack[sample_c == 9]\n\nidx = np.random.randint(0, len(x), size=64)\nbatch = x[idx]\n\nfig, ax = plt.subplots(8, 8, figsize=(16, 16))\n\nfor i, _ax in enumerate(ax):\n    for j, __ax in enumerate(_ax):\n        __ax.imshow(batch[8*i + j, ..., 0])\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Conclusion \nOn this notebook shows that using using lung slice attempt neural networks extracting features with lung property. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Reference:\n* [Image2Vec: AutoEncoder](https://www.kaggle.com/miklgr500/image2vec-autoencoder)\n* [[OSIC] Random slices from lung regions](https://www.kaggle.com/miklgr500/osic-random-slices-from-lung-regions)","execution_count":null}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}