SolveAcoustic.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "SolveAcoustic.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "machine_shape": "hm"
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rOVxI6WnFkS7"
      },
      "source": [
        "# Create dataset using artificial skull models"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XadC4xIjFw_Z"
      },
      "source": [
        "This notebook runs wavefield simulations containing artificial skull models, and captures images of the borders of the wavefield to create a dataset."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "4k2atFH31hPm",
        "outputId": "ae6cb730-d347-4d72-90ab-700a863d8424"
      },
      "source": [
        "!pip install git+https://github.com/devitocodes/devito.git --quiet\n",
        "!pip install sympy==1.7.0 --quiet\n",
        "!pip install numpy==1.17.0 --quiet"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "\u001b[K     |████████████████████████████████| 81kB 2.2MB/s \n",
            "\u001b[K     |████████████████████████████████| 102kB 6.3MB/s \n",
            "\u001b[K     |████████████████████████████████| 143kB 18.8MB/s \n",
            "\u001b[K     |████████████████████████████████| 51kB 6.6MB/s \n",
            "\u001b[K     |████████████████████████████████| 194kB 36.1MB/s \n",
            "\u001b[K     |████████████████████████████████| 51kB 8.0MB/s \n",
            "\u001b[K     |████████████████████████████████| 71kB 9.2MB/s \n",
            "\u001b[K     |████████████████████████████████| 71kB 8.3MB/s \n",
            "\u001b[?25h  Building wheel for devito (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Building wheel for py-cpuinfo (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Building wheel for cgen (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Building wheel for codepy (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Building wheel for pyrevolve (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Building wheel for pytools (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Building wheel for contexttimer (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "\u001b[31mERROR: pytest-cov 2.12.1 has requirement coverage>=5.2.1, but you'll have coverage 3.7.1 which is incompatible.\u001b[0m\n",
            "\u001b[31mERROR: pytest-cov 2.12.1 has requirement pytest>=4.6, but you'll have pytest 3.6.4 which is incompatible.\u001b[0m\n",
            "\u001b[K     |████████████████████████████████| 5.9MB 2.9MB/s \n",
            "\u001b[K     |████████████████████████████████| 20.3MB 1.4MB/s \n",
            "\u001b[31mERROR: tensorflow 2.5.0 has requirement numpy~=1.19.2, but you'll have numpy 1.17.0 which is incompatible.\u001b[0m\n",
            "\u001b[31mERROR: kapre 0.3.5 has requirement numpy>=1.18.5, but you'll have numpy 1.17.0 which is incompatible.\u001b[0m\n",
            "\u001b[31mERROR: datascience 0.10.6 has requirement folium==0.2.1, but you'll have folium 0.8.3 which is incompatible.\u001b[0m\n",
            "\u001b[31mERROR: albumentations 0.1.12 has requirement imgaug<0.2.7,>=0.2.5, but you'll have imgaug 0.2.9 which is incompatible.\u001b[0m\n",
            "\u001b[?25h"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "sDYMkPwEQgWT",
        "outputId": "cbf316a1-a9d8-44b5-9d28-091bc2d6578d"
      },
      "source": [
        "from google.colab import drive\n",
        "drive.mount('/content/drive', force_remount=True)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Mounted at /content/drive\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "lo4U6gon4gt4",
        "outputId": "a603cc53-ca6d-4d61-cf9f-8744c486e066"
      },
      "source": [
        "%cd /content/drive/My\\ Drive/ML_Ultrasound_Project/UROPdata"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/content/drive/My Drive/ML_Ultrasound_Project/UROPdata\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "iI9flbIi_8pQ",
        "outputId": "85b1e035-af56-458e-d925-4c02d4f76b34"
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import copy\n",
        "\n",
        "import devito\n",
        "from sympy import Symbol\n",
        "from devito import Operator\n",
        "from devito import Eq, solve\n",
        "from devito import TimeFunction\n",
        "\n",
        "from medical.plotting import plot_shotrecord, plot_velocity\n",
        "from medical.source import ToneBurstSource, Receiver\n",
        "from medical.source import TimeAxis\n",
        "from medical.model import Model\n",
        "\n",
        "from tqdm.notebook import tqdm\n",
        "\n",
        "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
        "\n",
        "# devito.parameters.configuration['dse'] = 'advanced'\n",
        "# devito.parameters.configuration['dle'] = 'speculative'\n",
        "devito.parameters.configuration['openmp'] = True\n",
        "devito.parameters.configuration['mpi'] = False\n",
        "\n",
        "def create_eliptic_array(n_sources, boundary_offset, dx, shape):\n",
        "    a = (shape[1] / 2)\n",
        "    b = (shape[0] / 2)\n",
        "\n",
        "    pos = np.empty((2, n_sources), dtype=np.float32)\n",
        "    angle = np.zeros((n_sources))\n",
        "    for idx in range(n_sources):\n",
        "        angle[idx] = ((idx) * 360.0 / n_sources) * 2 * np.pi / 360.0  # transducer angle\n",
        "        pos[1, idx] = (a * np.cos(angle[idx]) + shape[1] / 2)  # spherical equation of ellipse\n",
        "        pos[0, idx] = (b * np.sin(angle[idx]) + shape[0] / 2)  # spherical equation of ellipse\n",
        "\n",
        "    pos[0, :] += boundary_offset[0]\n",
        "    pos[1, :] += boundary_offset[1]\n",
        "    pos *= dx\n",
        "    # np.array((np.divide(grid[0] ** 2, major_ax ** 2) + np.divide(grid[1] ** 2, minor_ax ** 2)) < 1, dtype=np.float32)\n",
        "    return pos"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Trying to access deprecated config `openmp`. Using `language` instead\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "okAmaDkw1rnh"
      },
      "source": [
        "**Functions for creating samples**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "z9uOLMXV94Gi"
      },
      "source": [
        "def get_images(basedata, t, grid_shape):\n",
        "\n",
        "  # Here we get 9 images of the wavefield from the past, 6 timesteps apart,\n",
        "  # i.e. t, t-6, t-12, ..., t-48.\n",
        "  \n",
        "  if t < 6:\n",
        "    image0 = np.zeros(grid_shape)\n",
        "    image1 = np.zeros(grid_shape)\n",
        "    image2 = np.zeros(grid_shape)\n",
        "    image3 = np.zeros(grid_shape)\n",
        "    image4 = np.zeros(grid_shape)\n",
        "    image5 = np.zeros(grid_shape)\n",
        "    image6 = np.zeros(grid_shape)\n",
        "    image7 = np.zeros(grid_shape)\n",
        "    image8 = basedata[t]\n",
        "    \n",
        "  elif t < 12:\n",
        "    image0 = np.zeros(grid_shape)\n",
        "    image1 = np.zeros(grid_shape)\n",
        "    image2 = np.zeros(grid_shape)\n",
        "    image3 = np.zeros(grid_shape)\n",
        "    image4 = np.zeros(grid_shape)\n",
        "    image5 = np.zeros(grid_shape)\n",
        "    image6 = np.zeros(grid_shape)\n",
        "    image7 = basedata[t-6]\n",
        "    image8 = basedata[t]    \n",
        "\n",
        "  elif t < 18:\n",
        "    image0 = np.zeros(grid_shape)\n",
        "    image1 = np.zeros(grid_shape)\n",
        "    image2 = np.zeros(grid_shape)\n",
        "    image3 = np.zeros(grid_shape)\n",
        "    image4 = np.zeros(grid_shape)\n",
        "    image5 = np.zeros(grid_shape)\n",
        "    image6 = basedata[t-12]\n",
        "    image7 = basedata[t-6]\n",
        "    image8 = basedata[t]\n",
        "\n",
        "  elif t < 24:\n",
        "    image0 = np.zeros(grid_shape)\n",
        "    image1 = np.zeros(grid_shape)\n",
        "    image2 = np.zeros(grid_shape)\n",
        "    image3 = np.zeros(grid_shape)\n",
        "    image4 = np.zeros(grid_shape)\n",
        "    image5 = basedata[t-18]\n",
        "    image6 = basedata[t-12]\n",
        "    image7 = basedata[t-6]\n",
        "    image8 = basedata[t]\n",
        "\n",
        "  elif t < 30:\n",
        "    image0 = np.zeros(grid_shape)\n",
        "    image1 = np.zeros(grid_shape)\n",
        "    image2 = np.zeros(grid_shape)\n",
        "    image3 = np.zeros(grid_shape)\n",
        "    image4 = basedata[t-24]\n",
        "    image5 = basedata[t-18]\n",
        "    image6 = basedata[t-12]\n",
        "    image7 = basedata[t-6]\n",
        "    image8 = basedata[t]\n",
        "\n",
        "  elif t < 36:\n",
        "    image0 = np.zeros(grid_shape)\n",
        "    image1 = np.zeros(grid_shape)\n",
        "    image2 = np.zeros(grid_shape)\n",
        "    image3 = basedata[t-30]\n",
        "    image4 = basedata[t-24]\n",
        "    image5 = basedata[t-18]\n",
        "    image6 = basedata[t-12]\n",
        "    image7 = basedata[t-6]\n",
        "    image8 = basedata[t]\n",
        "\n",
        "  elif t < 42:\n",
        "    image0 = np.zeros(grid_shape)\n",
        "    image1 = np.zeros(grid_shape)\n",
        "    image2 = basedata[t-36]\n",
        "    image3 = basedata[t-30]\n",
        "    image4 = basedata[t-24]\n",
        "    image5 = basedata[t-18]\n",
        "    image6 = basedata[t-12]\n",
        "    image7 = basedata[t-6]\n",
        "    image8 = basedata[t]\n",
        "\n",
        "  elif t < 48:\n",
        "    image0 = np.zeros(grid_shape)\n",
        "    image1 = basedata[t-42]\n",
        "    image2 = basedata[t-36]\n",
        "    image3 = basedata[t-30]\n",
        "    image4 = basedata[t-24]\n",
        "    image5 = basedata[t-18]\n",
        "    image6 = basedata[t-12]\n",
        "    image7 = basedata[t-6]\n",
        "    image8 = basedata[t]\n",
        "\n",
        "  else:\n",
        "    image0 = basedata[t-48]\n",
        "    image1 = basedata[t-42]\n",
        "    image2 = basedata[t-36]\n",
        "    image3 = basedata[t-30]\n",
        "    image4 = basedata[t-24]\n",
        "    image5 = basedata[t-18]\n",
        "    image6 = basedata[t-12]\n",
        "    image7 = basedata[t-6]\n",
        "    image8 = basedata[t]\n",
        "\n",
        "  return image0, image1, image2, image3, image4, image5, image6, image7, image8"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_1U4iIC-BmDa"
      },
      "source": [
        "def create_samples(image0, image1, image2, image3, image4, image5, image6, image7, image8, i, slices, rot_num):\n",
        "\n",
        "  # This function creates network input images and targets from wavefield borders\n",
        "\n",
        "  # Create square 23x23 network input images\n",
        "  roi0 = np.array(np.rot90(image0[slices],rot_num))\n",
        "  square_timestep0 = roi0[i:i+23, 11:34]\n",
        "  roi1 = np.array(np.rot90(image1[slices],rot_num))\n",
        "  square_timestep1 = roi1[i:i+23, 11:34]    \n",
        "  roi2 = np.array(np.rot90(image2[slices],rot_num))\n",
        "  square_timestep2 = roi2[i:i+23, 11:34]\n",
        "  roi3 = np.array(np.rot90(image3[slices],rot_num))\n",
        "  square_timestep3 = roi3[i:i+23, 11:34]\n",
        "  roi4 = np.array(np.rot90(image4[slices],rot_num))\n",
        "  square_timestep4 = roi4[i:i+23, 11:34]\n",
        "  roi5 = np.array(np.rot90(image5[slices],rot_num))\n",
        "  square_timestep5 = roi5[i:i+23, 11:34]\n",
        "  roi6 = np.array(np.rot90(image6[slices],rot_num))\n",
        "  square_timestep6 = roi6[i:i+23, 11:34]\n",
        "  roi7 = np.array(np.rot90(image7[slices],rot_num))\n",
        "  square_timestep7 = roi7[i:i+23, 11:34]\n",
        "  roi8 = np.array(np.rot90(image8[slices],rot_num))\n",
        "  square_timestep8 = roi8[i:i+23, 11:34]\n",
        "\n",
        "  #link the 9 images together\n",
        "  newdata = np.vstack((square_timestep0, square_timestep1))\n",
        "  newdata = np.vstack((newdata, square_timestep2))\n",
        "  newdata = np.vstack((newdata, square_timestep3))\n",
        "  newdata = np.vstack((newdata, square_timestep4))\n",
        "  newdata = np.vstack((newdata, square_timestep5))\n",
        "  newdata = np.vstack((newdata, square_timestep6))\n",
        "  newdata = np.vstack((newdata, square_timestep7))\n",
        "  newdata = np.vstack((newdata, square_timestep8))\n",
        "\n",
        "  newdata = np.reshape(newdata, (9,23,23))\n",
        "\n",
        "  # Create corresponding labels\n",
        "  label = roi8[i+11:i+12, 0:11]\n",
        "  label_T = np.transpose(label)\n",
        "\n",
        "  return newdata, label_T\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XsU633y510nn"
      },
      "source": [
        "skulls = np.load('/content/drive/My Drive/ML_Ultrasound_Project/UROPdata/data/fakeSkulls.npy')\n",
        "\n",
        "y_length = 480+66+4  #480+11\n",
        "x_length = 350+66+4\n",
        "\n",
        "num = 2500              # Simulation last number of time steps - CHANGE THIS FOR THE AMOUNT OF TIMESTEPS\n",
        "\n",
        "shape = (1001, 1001) #(1001, 1001)            # Number of grid point (nx, nz) - CHANGE THIS FOR THE MEDIUM SIZE\n",
        "size_diff = int((shape[0]-1001)/2)\n",
        "\n",
        "top = (slice(225-1-11+size_diff, -226-11+11-size_diff), slice(290+size_diff, 290+34+size_diff))\n",
        "right = (slice(-226-34-size_diff, -226-size_diff), slice(290+19-11+size_diff, -291-11+11-size_diff))\n",
        "bottom = (slice(225-1-11+size_diff, -226-11+11-size_diff), slice(-291-34-size_diff, -291-size_diff))\n",
        "left = (slice(225+size_diff, 225+34+size_diff), slice(290+19-11+size_diff, -291-11+11-size_diff))\n",
        "\n",
        "area = (slice(225+size_diff, -226-size_diff), slice(290+size_diff, -291-size_diff))\n",
        "\n",
        "skull_start = 0\n",
        "for n, m in tqdm(enumerate(range(skull_start, len(skulls)))):\n",
        "  n+=skull_start\n",
        "  print(n)\n",
        "  model = skulls[m]\n",
        "  src_index = np.random.randint(512)  # Choose random point source location\n",
        "\n",
        "  # Define a physical size\n",
        "  spacing = (0.5e-3, 0.5e-3)      # Grid spacing in m.\n",
        "  origin = (0., 0.)               # What is the location of the top left corner. This is necessary to define\n",
        "                                # the absolute location of the source and receivers\n",
        "\n",
        "  # Define a velocity profile. The velocity is in m/s\n",
        "  diff_x = shape[0] - model.shape[0]\n",
        "  diff_y = shape[1] - model.shape[1]\n",
        "  offset_x = int(np.floor(diff_x/2))\n",
        "  offset_y = int(np.floor(diff_y/2))\n",
        "\n",
        "  model = np.pad(model,\n",
        "              ((int(np.floor(diff_x/2)), int(np.ceil(diff_x/2))), (int(np.floor(diff_y/2)), int(np.ceil(diff_y/2)))),\n",
        "              'edge')\n",
        "\n",
        "  plt.imshow(np.transpose(model), origin='lower')\n",
        "\n",
        "  # With the velocity and model size defined, we can create the seismic model that\n",
        "  # encapsulates this properties. We also define the size of the absorbing layer as 10 grid points\n",
        "  model = Model(vp=model, origin=origin, shape=shape, spacing=spacing, space_order=4, nbpml=0)\n",
        "\n",
        "  t0 = 0.                 # Simulation starts a t=0\n",
        "  dt = 0.08e-6            # Time step from model grid spacing\n",
        "\n",
        "  time_range = TimeAxis(start=t0, num=num, step=dt)\n",
        "\n",
        "  f0 = 0.40e6  # Source peak frequency is 400kHz\n",
        "  src = ToneBurstSource(name='src', grid=model.grid, f0=f0, npoint=1, time_range=time_range)\n",
        "\n",
        "  # First, position source centrally in all dimensions\n",
        "  src_positions = create_eliptic_array(512, (offset_x, offset_y), spacing[0], (480, 350))\n",
        "  src.coordinates.data[0, :] = src_positions[:, src_index]\n",
        "\n",
        "  # plt.scatter(src_positions[0, :]/spacing[0], src_positions[1, :]/spacing[1])\n",
        "  # plt.show()\n",
        "\n",
        "  # # We can plot the time signature to see the wavelet\n",
        "  # src.show()\n",
        "\n",
        "  # We can now show the source and receivers within our domain:\n",
        "  # Red dot: Source location\n",
        "  plot_velocity(model, source=src.coordinates.data)\n",
        "\n",
        "  # In order to represent the wavefield u and the square slowness we need symbolic objects\n",
        "  # corresponding to time-space-varying field (u, TimeFunction) and\n",
        "  # space-varying field (m, Function)\n",
        "\n",
        "  # Define the wavefield with the size of the model and the time dimension\n",
        "  u = TimeFunction(name=\"u\", grid=model.grid, time_order=2, space_order=4, save=num)\n",
        "\n",
        "  # Create a temporary symbol for H to avoid expensive sympy solve\n",
        "  H = Symbol('H')\n",
        "\n",
        "  # Define PDE\n",
        "  eq = model.m * u.dt2 - H\n",
        "\n",
        "  # Solve the symbolic equation for the field to be updated\n",
        "  eq_time = solve(eq, u.forward)\n",
        "\n",
        "  # Get the spacial FD\n",
        "  biharmonic = u.biharmonic(1/model.m)\n",
        "  laplacian = u.laplace + dt**2/12 * biharmonic\n",
        "\n",
        "  # This discrete PDE can be solved in a time-marching way updating u(t+dt) from the previous time step\n",
        "  # Devito as a shortcut for u(t+dt) which is u.forward. We can then rewrite the PDE as\n",
        "  # a time marching updating equation known as a stencil using customized SymPy functions\n",
        "  stencil = Eq(u.forward, eq_time.subs({H: laplacian}))\n",
        "\n",
        "  # Finally we define the source injection function to generate the corresponding code\n",
        "  src_term = src.inject(field=u.forward, expr=src * dt**2 / model.m)\n",
        "\n",
        "  # Execute the operator\n",
        "  u.data.fill(0.)\n",
        "  op = Operator([stencil] + src_term, subs=model.spacing_map)\n",
        "  op(dt=dt)  # AFTER THIS LINE u.data HAS THE WAVEFIELD AT ALL TIMESTEPS WITH SHAPE (num, nx, ny)\n",
        "\n",
        "  skull_images = []\n",
        "  skull_labels = []\n",
        "\n",
        "  for t in tqdm(range(0, num, 20)):\n",
        "    if np.any(u.data[t][top] != 0):\n",
        "      image0, image1, image2, image3, image4, image5, image6, image7, image8 = get_images(u.data, t, shape)\n",
        "      for i in range(540):\n",
        "        ts_images, ts_label = create_samples(image0, image1, image2, image3, image4, image5, image6, image7, image8, i, top, 0)\n",
        "        skull_images.append(ts_images)\n",
        "        skull_labels.append(ts_label)\n",
        "        if (i+1)%30 == 0:\n",
        "          if np.all(np.abs(np.array(skull_images[-30:])) < 0.5e-9):\n",
        "            skull_images = skull_images[:-30]\n",
        "            skull_labels = skull_labels[:-30]\n",
        "          else:\n",
        "            print('top')\n",
        "            plt.figure()\n",
        "            plt.imshow(skull_images[-12][0,:,:])\n",
        "            plt.show()\n",
        "            plt.figure()\n",
        "            plt.imshow(skull_images[-12][8,:,:])\n",
        "            plt.show()\n",
        "          #   plt.figure()\n",
        "          #   plt.imshow(np.hstack((np.squeeze(skull_labels[-23:]), skull_images[-12][8,:,:])))\n",
        "          #   plt.colorbar()\n",
        "          #   plt.show()\n",
        "\n",
        "    if np.any(u.data[t][right] != 0):\n",
        "      image0, image1, image2, image3, image4, image5, image6, image7, image8 = get_images(u.data, t, shape)\n",
        "      for i in range(390):\n",
        "        ts_images, ts_label = create_samples(image0, image1, image2, image3, image4, image5, image6, image7, image8, i, right, 3)\n",
        "        skull_images.append(ts_images)\n",
        "        skull_labels.append(ts_label)\n",
        "        if (i+1)%30 == 0:\n",
        "          if np.all(np.abs(np.array(skull_images[-30:])) < 0.5e-9):\n",
        "            skull_images = skull_images[:-30]\n",
        "            skull_labels = skull_labels[:-30]\n",
        "          # else:\n",
        "          #   print('right')\n",
        "          #   plt.figure()\n",
        "          #   plt.imshow(np.hstack((np.squeeze(skull_labels[-23:]), skull_images[-12][8,:,:])))\n",
        "          #   plt.colorbar()\n",
        "          #   plt.show()\n",
        "\n",
        "    if np.any(u.data[t][bottom] != 0):\n",
        "      image0, image1, image2, image3, image4, image5, image6, image7, image8 = get_images(u.data, t, shape)\n",
        "      for i in range(540):\n",
        "        ts_images, ts_label = create_samples(image0, image1, image2, image3, image4, image5, image6, image7, image8, i, bottom, 2)\n",
        "        skull_images.append(ts_images)\n",
        "        skull_labels.append(ts_label)\n",
        "        if (i+1)%30 == 0:\n",
        "          if np.all(np.abs(np.array(skull_images[-30:])) < 0.5e-9):\n",
        "            skull_images = skull_images[:-30]\n",
        "            skull_labels = skull_labels[:-30]\n",
        "          # else:\n",
        "          #   print('bottom')\n",
        "          #   plt.figure()\n",
        "          #   plt.imshow(np.hstack((np.squeeze(skull_labels[-23:]), skull_images[-12][8,:,:])))\n",
        "          #   plt.colorbar()\n",
        "          #   plt.show()\n",
        "\n",
        "    if np.any(u.data[t][left] != 0):\n",
        "      image0, image1, image2, image3, image4, image5, image6, image7, image8 = get_images(u.data, t, shape)\n",
        "      for i in range(390):\n",
        "        ts_images, ts_label = create_samples(image0, image1, image2, image3, image4, image5, image6, image7, image8, i, left, 1)\n",
        "        skull_images.append(ts_images)\n",
        "        skull_labels.append(ts_label)\n",
        "        if (i+1)%30 == 0:\n",
        "          if np.all(np.abs(np.array(skull_images[-30:])) < 0.5e-9):\n",
        "            skull_images = skull_images[:-30]\n",
        "            skull_labels = skull_labels[:-30]\n",
        "          # else:\n",
        "          #   print('left')\n",
        "          #   plt.figure()\n",
        "          #   plt.imshow(np.hstack((np.squeeze(skull_labels[-23:]), skull_images[-12][8,:,:])))\n",
        "          #   plt.colorbar()\n",
        "          #   plt.show()\n",
        "\n",
        "  skull_images = np.array(skull_images)\n",
        "  skull_labels = np.array(skull_labels)\n",
        "\n",
        "  np.savez_compressed('/content/drive/My Drive/ML_Ultrasound_Project/UROPdata/skull_images_{}'.format(n), a=skull_images)\n",
        "  np.savez_compressed('/content/drive/My Drive/ML_Ultrasound_Project/UROPdata/skull_labels_{}'.format(n), a=skull_labels)\n",
        "\n",
        "  print('{} samples'.format(len(skull_labels)))"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}