Boundary problems using a numpy kernel and padding on a matrix within a C functions called from python

Question

I want to use an external c function a numpy kernel operation on an image. For this the goal is to go through the image image only within the padding so that the kernel does not go outside the image. In the following debugging code, I first generate a matrix filled with 1 in the python script. I then pass this matrix with its dimension to the c function. I iterate over rows and column and want to set the value to 9 if I am on the border of the matrix (in the padding) and to 8 if I am within the part of the matrix where I want to do the computation.

The eventual goal of the computation is to do a computation for every entry of the matrix taking all the neighbours (e.g. on a 3x3 matrix pattern) for that computation. For a 3x3 matrix, the search has to start before the stop after the first row/column and stop before the last row column.

The result I want to have is the following:

,0,1,2,3
0,9,9,9,9
1,9,8,8,9
2,9,8,8,9
3,9,9,9,9

Instead I get the following:

,0,1,2,3
0,9,9,9,9
1,8,8,8,8
2,8,8,8,8
3,8,1,1,1

Here is the .sh script debug_scriptsh.sh:

a=debug_scriptc && cc -fPIC -shared -o $a.so $a.c && python debug_scriptpy.py && open debug_table.csv

the .py script debug_scriptc.c:

import ctypes
import pandas as pd
import numpy as np
import sys
import pdb

DLL_NAME = "./debug_scriptc.{:s}".format("dll" if sys.platform[:3].lower() == "win" else "so")

def np_mat_type(rows, cols, element_type=float):
  return np.ctypeslib.ndpointer(dtype=element_type, shape=(rows, cols), flags="C_CONTIGUOUS")

rows0 = 4
cols0 = 4
kernel_size = 3

dll = ctypes.CDLL(DLL_NAME)
matrix_func = dll.matrixFunc4matin
matrix_func.argtypes = (
    ctypes.c_size_t,
    np_mat_type(rows0, cols0), 
    ctypes.c_size_t, 
    ctypes.c_size_t
    )
matrix_func.restype = np_mat_type(rows0, cols0)

dealloc_array = dll.deallocArray
dealloc_array.argtypes = (np_mat_type(rows0, cols0),)
dealloc_array.restype = None
u_rms_array = np.zeros((rows0, cols0))

velocity_image_ones = np.empty((rows0, cols0))
velocity_image_ones[:] = 1

u_rms_array = matrix_func(
    kernel_size, 
    velocity_image_ones,
    rows0, 
    cols0 
    )
debug_DF_u_rms_array_C = pd.DataFrame(u_rms_array.astype(int))
debug_DF_u_rms_array_C.to_csv("debug_table.csv")
dealloc_array(u_rms_array)

the .c script debug_scriptc.c:


#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#if defined(_WIN32)
#  define DLL00_EXPORT_API __declspec(dllexport)
#else
#  define DLL00_EXPORT_API
#endif

#define ELEMENT_TYPE double

#if defined(__cplusplus)
extern "C" {
#endif

DLL00_EXPORT_API ELEMENT_TYPE* matrixFunc4matin(
  size_t kernel_size,
  const ELEMENT_TYPE *pvelocity_image_ones, 
  size_t rows0, 
  size_t cols0
  );
DLL00_EXPORT_API void deallocArray(ELEMENT_TYPE *pMat);

#if defined(__cplusplus)
}
#endif

ELEMENT_TYPE* matrixFunc4matin(
  size_t kernel_size,
  const ELEMENT_TYPE *pvelocity_image_ones, 
  size_t rows0, 
  size_t cols0 
  )

{
  size_t matSize = sizeof(ELEMENT_TYPE) * cols0 * rows0;
  ELEMENT_TYPE *pRet = (ELEMENT_TYPE*)(malloc(matSize));

  memcpy(pRet, pvelocity_image_ones, matSize);  

  int start_kern_search = kernel_size - 2;
  int end_kern_search_rows = rows0 - (kernel_size-1)/2;
  int end_kern_search_cols = cols0 - (kernel_size-1)/2;

  for (int i_velmap = 0; i_velmap < rows0; i_velmap++){  
    for (int j_velmap = 0; j_velmap < cols0; j_velmap++){  
      if (i_velmap < start_kern_search | i_velmap > end_kern_search_rows | j_velmap < start_kern_search | j_velmap > end_kern_search_cols){
        pRet[i_velmap * cols0 + j_velmap] = 9;
        continue;
      }
      else
      {
        pRet[i_velmap * end_kern_search_cols + j_velmap] = 8;
      }
    }
  }
  return pRet;
}


void deallocArray(ELEMENT_TYPE *pMat)
{
    if (pMat)
        free(pMat);
}

Thank you for the comment @hyde. I changed the title and the start of the question. Is it clearer now? — ecjb, Oct 23 '22 at 14:22
Expected results can't be obtained using current data. `rows0 - kernel_size = 2` (so that you have 1 before and 1 after the kernel size). Same for *cols0*. What would be the expected result for greater differences (e.g. 5)? — CristiFati, Oct 23 '22 at 17:31
Thanks for your comment @CristiFati. I changed the line by `rows0 - (kernel_size-1)/2`. I also simplified the numbers filled in and the precision so that the code is more readable. You cannot obtain the same result as I do with the code??? — ecjb, Oct 23 '22 at 17:58
Kernel size seems irrelevant here. What does it represent? Is the inner matrix (with the useful data) ***kernel\_size X kernel\_size*** dimensioned? The outer matrix is ***rows X cols*** large. Could you clarify on that? What I meant in my prev comment is that the problem is incorrectly formed (assuming what I said above is correct). — CristiFati, Oct 23 '22 at 17:59
@CristiFati I added the following explanation to the question "The eventual goal of the computation is to do a computation for every entry of the matrix taking all the neighbours (e.g. on a 3x3 matrix pattern) for that computation. For a 3x3 matrix, the search has to start before the stop after the first row/column and stop before the last row column." is it clearer now? — ecjb, Oct 23 '22 at 18:06
I don't know how to initiate chat. Unfortunately I didn't get it. (I still don;t know what the *kernel\_size* is). Are we always talking about a square matrix? What would the required output be on *2x2*, *3x3*, *6x6* (*5x4* if square is not required) matrices? — CristiFati, Oct 23 '22 at 18:12
@CristiFati. Yes exactly! we are always talking about a square matrix (with a side (in that case the kernel size) uneven, e.g. 3x3, 5x5, 7x7 which corresonds to the "kernel_size" of the code to 3, 5 or 7, respectively). Does that get clearer? is there somebody still not clear? — ecjb, Oct 23 '22 at 18:42

CristiFati · Accepted Answer · 2022-10-23T19:43:57.890

Continuation of [SO]: How to write a function in C which takes as input two matrices as numpy array.

The problem is index calculation.
When framing (centering) a (small) rectangle (inner), into a larger one (outer), the edge formula (on each axis (dimension)) is:

edge = (outer_dim - inner_dim) / 2

dll00.c:

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#if defined(_WIN32)
#  define DLL00_EXPORT_API __declspec(dllexport)
#else
#  define DLL00_EXPORT_API
#endif

#define ELEMENT_TYPE double

#if defined(__cplusplus)
extern "C" {
#endif

DLL00_EXPORT_API ELEMENT_TYPE* matrixFormat(const ELEMENT_TYPE *pMat, size_t rows, size_t cols, size_t kernel);
DLL00_EXPORT_API void deallocArray(ELEMENT_TYPE *pMat);

#if defined(__cplusplus)
}
#endif

#define MIN_EDGE 1
#define PAD_VAL 9
#define KERNEL_VAL 8


static size_t edge(size_t dim, size_t kernel, size_t minEdge)
{
    if (kernel >= dim)
        return minEdge;
    size_t ret = (dim - kernel) / 2;
    return (ret > minEdge) ? ret : minEdge;
}


ELEMENT_TYPE* matrixFormat(
    const ELEMENT_TYPE *pMat,
    size_t rows,
    size_t cols,
    size_t kernel)
{
    size_t matSize = sizeof(ELEMENT_TYPE) * cols * rows;
    ELEMENT_TYPE *pRet = (ELEMENT_TYPE*)(malloc(matSize));

    printf("\nC: - matrix (%zu X %zu (%zu)):\n", rows, cols, kernel);
    size_t minRow = edge(rows, kernel, MIN_EDGE),
           minCol = edge(cols, kernel, MIN_EDGE);
    size_t maxRow = rows - minRow - 1,
           maxCol = cols - minCol - 1;
    for (size_t i = 0; i < rows; ++i) {
        for (size_t j = 0; j < cols; ++j) {
            pRet[i * cols + j] = ((i >= minRow) && (i <= maxRow) && (j >= minCol) && (j <= maxCol)) ? KERNEL_VAL : PAD_VAL;
        }
    }
    return pRet;
}


void deallocArray(ELEMENT_TYPE *pMat)
{
    if (pMat)
        free(pMat);
}

code00.py:

#!/usr/bin/env python

import ctypes as cts
import sys

import numpy as np


DLL_NAME = "./dll00.{:s}".format("dll" if sys.platform[:3].lower() == "win" else "so")

ELEMENT_TYPE = cts.c_double


def np_mat_type(rows, cols, element_type=ELEMENT_TYPE):
    return np.ctypeslib.ndpointer(dtype=element_type, shape=(rows, cols), flags="C_CONTIGUOUS")


def main(*argv):
    np.set_printoptions(precision=3)

    dll = cts.CDLL(DLL_NAME)
    matrix_format = dll.matrixFormat
    dealloc_array = dll.deallocArray
    dealloc_array.restype = None

    values = (
        (4, 4, 3),
        (6, 6, 3),
        #(3, 3, 3),
    )

    for rows, cols, kernel in values:
        print("\nRows: {:d}, Cols: {:d}, Kernel: {:d}".format(rows, cols, kernel))
        mat = np.ones((rows, cols), dtype=ELEMENT_TYPE)

        matrix_format.argtypes = (np_mat_type(rows, cols), cts.c_size_t, cts.c_size_t, cts.c_size_t)
        matrix_format.restype = np_mat_type(rows, cols)
        dealloc_array.argtypes = (np_mat_type(rows, cols),)

        print("mat:")
        print(mat)
        mat_res = matrix_format(mat, rows, cols, kernel)
        print("result:")
        print(mat_res)

        dealloc_array(mat_res)


if __name__ == "__main__":
    print("Python {:s} {:03d}bit on {:s}\n".format(" ".join(elem.strip() for elem in sys.version.split("\n")),
                                                   64 if sys.maxsize > 0x100000000 else 32, sys.platform))
    rc = main(*sys.argv[1:])
    print("\nDone.\n")
    sys.exit(rc)

Output:

(qaic-env) [cfati@cfati-5510-0:/mnt/e/Work/Dev/StackOverflow/q074171783]> ~/sopr.sh
### Set shorter prompt to better fit when pasted in StackOverflow (or other) pages ###

[064bit prompt]> ls
code00.py  dll00.c
[064bit prompt]> gcc -fPIC -shared -o dll00.so dll00.c
[064bit prompt]> ls
code00.py  dll00.c  dll00.so
[064bit prompt]> python ./code00.py
Python 3.8.10 (default, Jun 22 2022, 20:18:18) [GCC 9.4.0] 064bit on linux


Rows: 4, Cols: 4, Kernel: 3
mat:
[[1. 1. 1. 1.]
 [1. 1. 1. 1.]
 [1. 1. 1. 1.]
 [1. 1. 1. 1.]]

C: - matrix (4 X 4 (3)):
result:
[[9. 9. 9. 9.]
 [9. 8. 8. 9.]
 [9. 8. 8. 9.]
 [9. 9. 9. 9.]]

Rows: 6, Cols: 6, Kernel: 3
mat:
[[1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1.]]

C: - matrix (6 X 6 (3)):
result:
[[9. 9. 9. 9. 9. 9.]
 [9. 8. 8. 8. 8. 9.]
 [9. 8. 8. 8. 8. 9.]
 [9. 8. 8. 8. 8. 9.]
 [9. 8. 8. 8. 8. 9.]
 [9. 9. 9. 9. 9. 9.]]

Done.

I did not understand anything yet. I have to process it more. But until then... MANY THANKS!!!!! :) @CristiFati — ecjb, Oct 23 '22 at 19:53
You're welcome! Everything is in the *edge* function. It's a bit complex to avoid any inconveniences that could be introduced by unsigned types. The principle is the same as when trying to position a window (h X w) to be centered on the screen (hs X ws). Top/left coordinates are calculated based on that formula. — CristiFati, Oct 23 '22 at 19:58

Boundary problems using a numpy kernel and padding on a matrix within a C functions called from python

1 Answers1

Linked