6_pycuda_conv.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "6_pycuda_conv.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyMLXhYq89c015k4zkiMAUZW",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/geoffwoollard/gpu-speedups-mbptechtalk2020/blob/master/6_pycuda_conv.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WEEBPaNpyWZX",
        "colab_type": "text"
      },
      "source": [
        "# MBP Tech Talk 2020 :: PyCUDA Image Convolution\n",
        "From a [pedagogical GitHub repo](https://github.com/jtc42/pycuda-convolution) last updated in 2017"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Jn92jD91zc-2",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!pip install pycuda"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8u8xXE2D0Giv",
        "colab_type": "text"
      },
      "source": [
        "# Image"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vTnfAH7F1Usc",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import skimage.data\n",
        "from skimage.color import rgb2gray\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline\n",
        "\n",
        "full_image = rgb2gray(skimage.data.coffee()).astype(np.float32) #* 255\n",
        "plt.figure()\n",
        "plt.imshow(full_image, cmap='gray')\n",
        "plt.title(\"Full size image:\")\n",
        "image = full_image[150:350, 200:400].copy() # We don't want a view but an array and therefore use copy()\n",
        "plt.figure()\n",
        "plt.imshow(image, cmap='gray')\n",
        "plt.title(\"Part of the image we use:\")\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pzyrwgpGwd4g",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eLNTlJW4yKRV",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# kernel.cu in the repo\n",
        "kernel_cu = '''\n",
        "__global__ void conv(const float *A, const float *B, int aw, int ah, int bw, int bh, float b_sum, float *C){\n",
        "\n",
        "    /*Get row and column to operate on from thread coordinates*/\n",
        "    int tx = threadIdx.x;\n",
        "    int ty = threadIdx.y;\n",
        "    \n",
        "    int bx = blockIdx.x;\n",
        "    int by = blockIdx.y;\n",
        "    \n",
        "    int row = by*blockDim.y + ty;\n",
        "    int col = bx*blockDim.x + tx;\n",
        "    \n",
        "    /*Calculate \"padding\" radius of convolution kernel (distance around central pixel)*/\n",
        "    int pw = (bw-1)/2;\n",
        "    int ph = (bh-1)/2;\n",
        "\n",
        "    /*If within the range of C (ie A - padding)*/\n",
        "    if( row < (ah-2*ph) && col < (aw-2*pw) ) {\n",
        "        \n",
        "        /*Set initial pixel value*/\n",
        "        float val = 0; // change to float so that does not round down normalized image to all zeros\n",
        "        \n",
        "         /*For each vertical position on the kernel matrix, relative to the central pixel*/\n",
        "        for(int i=-ph; i<=ph; i=i+1){\n",
        "            /*Calculate zero-indexed row ID on kernel matrix*/\n",
        "            int b_row = i+ph; \n",
        "\n",
        "            /*For each horizontal position on the kernel matrix, relative to the central pixel*/\n",
        "            for(int j=-pw; j<=pw; j=j+1){\n",
        "                /*Calculate zero-indexed column ID on kernel matrix*/\n",
        "                int b_col = j+pw;\n",
        "\n",
        "                /*Add product of kernel value and corresponding image value to running total*/\n",
        "                val += A[ (row+ph +i)*aw + (col+pw +j) ] * B[ b_row*bw + b_col ];\n",
        "            }\n",
        "        }\n",
        "        \n",
        "        /*Copy appropriately normalised resulting pixel value to position on C matrix*/\n",
        "        C[row*(aw-2*pw) + col] = val/b_sum;\n",
        "    }\n",
        "}\n",
        "'''"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3X2HwUnizmUC",
        "colab_type": "text"
      },
      "source": [
        "Setup the kernel"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yC2I2A3fzOd8",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "\n",
        "import pycuda.driver as cuda\n",
        "import pycuda.autoinit\n",
        "from pycuda.compiler import SourceModule\n",
        "\n",
        "# DEVICE SETUP\n",
        "BLOCK_SIZE = 32  # Max 32. 32**2 = 1024, max for GTX1060\n",
        "    \n",
        "# Compile kernel\n",
        "#mod = SourceModule(open(\"kernel.cu\", \"r\").read())\n",
        "mod = SourceModule(kernel_cu)\n",
        "\n",
        "# Get functions\n",
        "conv = mod.get_function(\"conv\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ItzbgjEUzfn2",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# this python function is a wrapper around the conv kernel function\n",
        "def convolve(a, b):\n",
        "    global BLOCK_SIZE\n",
        "    global conv\n",
        "    \n",
        "    a, b = [np.array(i).astype(np.float32) for i in [a, b]]\n",
        "    \n",
        "    # Matrix A \n",
        "    aw = np.int32(a.shape[1])  # Widthof in matrix\n",
        "    ah = np.int32(a.shape[0])  # Height of in matrix\n",
        "    \n",
        "    # Matrix B (kernel)\n",
        "    bw = np.int32(b.shape[1])  # Widthof in matrix\n",
        "    if bw % 2 == 0:\n",
        "        print(\"Kernel width is not an odd number! Strange things will happen...\")\n",
        "    bh = np.int32(b.shape[0])  # Height of in matrix\n",
        "    if bh % 2 == 0:\n",
        "        print(\"Kernel height is not an odd number! Strange things will happen...\")\n",
        "    b_sum = np.absolute(b).sum() #np.int32(np.absolute(b).sum())\n",
        "    \n",
        "    # Matrix C, subtract 2*padding, *2 because it's taken off all sides\n",
        "    c = np.empty([ah-(bh-1), aw-(bw-1)])\n",
        "    c = c.astype(np.float32)\n",
        "    \n",
        "    # Allocate memory on device\n",
        "    a_gpu = cuda.mem_alloc(a.nbytes)\n",
        "    b_gpu = cuda.mem_alloc(b.nbytes)\n",
        "    c_gpu = cuda.mem_alloc(c.nbytes)\n",
        "    \n",
        "    # Copy matrix to memory\n",
        "    cuda.memcpy_htod(a_gpu, a)\n",
        "    cuda.memcpy_htod(b_gpu, b)\n",
        "\n",
        "    # Set grid size from A matrix\n",
        "    grid = (int(aw/BLOCK_SIZE+(0 if aw % BLOCK_SIZE is 0 else 1)), \n",
        "            int(ah/BLOCK_SIZE+(0 if ah % BLOCK_SIZE is 0 else 1)), \n",
        "                          1)\n",
        "    \n",
        "    # Call gpu function\n",
        "    conv(a_gpu, b_gpu, aw, ah, bw, bh, b_sum, c_gpu, block=(BLOCK_SIZE, BLOCK_SIZE, 1), grid=grid)\n",
        "    \n",
        "    # Copy back the result\n",
        "    cuda.memcpy_dtoh(c, c_gpu)\n",
        "    \n",
        "    # Free memory. May not be useful? Ask about this.\n",
        "    a_gpu.free()\n",
        "    b_gpu.free()\n",
        "    c_gpu.free()\n",
        "    \n",
        "    # Return the result\n",
        "    return c"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RYW32GDXzkij",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "mask = np.ones((3,3))\n",
        "c = convolve(image, mask)\n",
        "\n",
        "plt.figure()\n",
        "plt.imshow(c, cmap='gray')\n",
        "plt.title(\"Convolved image\")\n",
        "plt.show()\n",
        "c"
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}