Python-сообщество

from cgriddata import griddata
...
...
# x = ndarray, dim=2, type=float32
# y = ndarray, dim=2, type=float32
# v = ndarray, dim=2, type=float32
# griddata(x, y, v, minimum_value, scale_factor) -> dict({'min_x':float, 'min_y':float, 'max_x':float, 'max_y':float, 'values':ndarray_float32}) OR None
data = griddata(x, y, v, 0, 0.001)
if data:  # Здесь БЫВАЕТ ошибка, а бывает всё нормально
    ...

#include <Python.h>
#include <numpy/arrayobject.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>
typedef struct
{
    int x;
    int y;
} TPoint;
typedef struct
{
    float x;
    float y;
    float v;
} TPointF;
float min3(float a, float b, float c)
{
    if (a < b) return (a<c?a:c);
    else return (b<c?b:c);
}
float max3(float a, float b, float c)
{
    if (a > b) return (a>c?a:c);
    else return (b>c?b:c);
}
float min4(float a, float b, float c, float d)
{
    if (a < b)
    {
        if (a < c) return (a<d?a:d);
        else return (c<d?c:d);
    }
    else
    {
        if (b < c) return (b<d?b:d);
        else return (c<d?c:d);
    }
}
float max4(float a, float b, float c, float d)
{
    if (a > b)
    {
        if (a > c) return (a>d?a:d);
        else return (c>d?c:d);
    }
    else
    {
        if (b > c) return (b>d?b:d);
        else return (c>d?c:d);
    }
}
bool in_triangle(TPoint p, TPointF a, TPointF b, TPointF c)
{
    float ab = (a.x-p.x)*(b.y-a.y) - (b.x-a.x)*(a.y-p.y);
    float bc = (b.x-p.x)*(c.y-b.y) - (c.x-b.x)*(b.y-p.y);
    float ca = (c.x-p.x)*(a.y-c.y) - (a.x-c.x)*(c.y-p.y);
    if ((ab<=0 && bc<=0 && ca<=0) || (ab>=0 && bc>=0 && ca>=0)) return true;
    else return false;
}
float interpolation(TPoint p, TPointF a, TPointF b, TPointF c)
{
    float x_delta = p.x - a.x;
    float y_delta = p.y - a.y;
    float w_a = 1 / sqrtf(x_delta*x_delta + y_delta*y_delta);
    x_delta = p.x - b.x;
    y_delta = p.y - b.y;
    float w_b = 1 / sqrtf(x_delta*x_delta + y_delta*y_delta);
    x_delta = p.x - c.x;
    y_delta = p.y - c.y;
    float w_c = 1 / sqrtf(x_delta*x_delta + y_delta*y_delta);
    return (a.v*w_a + b.v*w_b + c.v*w_c) / (w_a+w_b+w_c);
}
static PyObject*
griddata(PyObject* self, PyObject* args)
{
    PyArrayObject* x_array = NULL;
    PyArrayObject* y_array = NULL;
    PyArrayObject* values_array = NULL;
    float minimum_value;
    float scale_factor;
    if (PyArg_ParseTuple(args, "O!O!O!ff",
        &PyArray_Type, &x_array,
        &PyArray_Type, &y_array,
        &PyArray_Type, &values_array,
        &minimum_value, &scale_factor) == 0) return NULL;
    /* Input validation */
    if (PyArray_TYPE(x_array) != NPY_FLOAT ||
        PyArray_TYPE(y_array) != NPY_FLOAT ||
        PyArray_TYPE(values_array) != NPY_FLOAT) return NULL;
    if (PyArray_NDIM(x_array) != 2 || PyArray_NDIM(y_array) != 2 || PyArray_NDIM(values_array) != 2) return NULL;
    int dims[2];
    dims[0] = x_array->dimensions[0];
    dims[1] = x_array->dimensions[1];
    if (! (dims[0] > 0 && dims[1] > 0)) return NULL;
    if (y_array->dimensions[0] != dims[0] || y_array->dimensions[1] != dims[1] ||
        values_array->dimensions[0] != dims[0] || values_array->dimensions[1] != dims[1]) return NULL;
    /* Find the minimum and maximum coordinates */
    float minimum_x;
    float minimum_y;
    float maximum_x;
    float maximum_y;
    {
        int i = 0;
        bool finded = false;
        for (; i < dims[0]; ++i)
        {
            for (int j = 0; j < dims[1]; ++j)
            {
                if (*(float*)PyArray_GETPTR2(values_array, i, j) < minimum_value) continue;
                float x = *(float*)PyArray_GETPTR2(x_array, i, j);
                float y = *(float*)PyArray_GETPTR2(y_array, i, j);
                minimum_x = maximum_x = x;
                minimum_y = maximum_y = y;
                finded = true;
                break;
            }
            if (finded) break;
        }
        if (! finded) Py_RETURN_NONE;
        finded = false;
        for (; i < dims[0]; ++i)
        {
            for (int j = 0; j < dims[1]; ++j)
            {
                if (*(float*)PyArray_GETPTR2(values_array, i, j) < minimum_value) continue;
                float x = *(float*)PyArray_GETPTR2(x_array, i, j);
                float y = *(float*)PyArray_GETPTR2(y_array, i, j);
                if (x < minimum_x) minimum_x = x;
                if (y < minimum_y) minimum_y = y;
                if (x > maximum_x) maximum_x = x;
                if (y > maximum_y) maximum_y = y;
                finded = true;
            }
        }
        if (! finded) Py_RETURN_NONE;
    }
    /* Transformation of coordinates */
    float* new_x_array = malloc(sizeof(float)*dims[0]*dims[1]);
    if (new_x_array == NULL) return NULL;
    float* new_y_array = malloc(sizeof(float)*dims[0]*dims[1]);
    if (new_y_array == NULL)
    {
        free(new_x_array);
        return NULL;
    }
    for (int i=0; i < dims[0]; ++i)
    {
        for (int j=0; j < dims[1]; ++j)
        {
            if (*(float*)PyArray_GETPTR2(values_array, i, j) < minimum_value) continue;
            float x = *(float*)PyArray_GETPTR2(x_array, i, j);
            float y = *(float*)PyArray_GETPTR2(y_array, i, j);
            int index = dims[1]*i + j;
            new_x_array[index] = (x - minimum_x) * scale_factor;
            new_y_array[index] = (y - minimum_y) * scale_factor;
        }
    }
    /* Creating the output array */
    int width = (int) ceilf((maximum_x - minimum_x) * scale_factor);
    int height = (int) ceilf((maximum_y - minimum_y) * scale_factor);
    if (width <= 0 || height <= 0) return NULL;
    int new_dims[2];
    new_dims[0] = height;
    new_dims[1] = width;
    PyArrayObject* output_array = PyArray_FromDims(2, new_dims, NPY_FLOAT);
    for (int i=0; i<height; ++i) for (int j=0; j<=width; ++j) *(float*)PyArray_GETPTR2(output_array, i,j) = minimum_value-1;
    /* Interpolation */
    bool is_empty = true;
    for (int i=0; i+1 < dims[0]; ++i)
    {
        for (int j=0; j+1 < dims[1]; ++j)
        {
            int index00 = dims[1]*i + j;
            int index01 = dims[1]*i + j+1;
            int index10 = dims[1]*(i+1) + j;
            int index11 = dims[1]*(i+1) + j+1;
            float v00 = *(float*)PyArray_GETPTR2(values_array, i,j);
            float v01 = *(float*)PyArray_GETPTR2(values_array, i,j+1);
            float v10 = *(float*)PyArray_GETPTR2(values_array, i+1,j);
            float v11 = *(float*)PyArray_GETPTR2(values_array, i+1,j+1);
            TPointF p00 = {.x = new_x_array[index00], .y = new_y_array[index00], .v = v00};
            TPointF p01 = {.x = new_x_array[index01], .y = new_y_array[index01], .v = v01};
            TPointF p10 = {.x = new_x_array[index10], .y = new_y_array[index10], .v = v10};
            TPointF p11 = {.x = new_x_array[index11], .y = new_y_array[index11], .v = v11};
            if (v00 < minimum_value)
            {
                if (v01 < minimum_value || v10 < minimum_value || v11 < minimum_value) continue;
                int begin_x = (int) ceilf(min3(p01.x, p10.x, p11.x));
                int end_x = (int) floorf(max3(p01.x, p10.x, p11.x));
                int begin_y = (int) ceilf(min3(p01.y, p10.y, p11.y));
                int end_y = (int) floorf(max3(p01.y, p10.y, p11.y));
                if (begin_x < 0) begin_x = 0;
                if (end_x >= width) end_x = width-1;
                if (begin_y < 0) begin_y = 0;
                if (end_y >= height) end_y = height-1;
                for (int y = begin_y; y <= end_y; ++y)
                {
                    for (int x = begin_x; x <= end_x; ++x)
                    {
                        TPoint p = {.x = x, .y = y};
                        if (in_triangle(p, p01, p10, p11))
                        {
                            *(float*)PyArray_GETPTR2(output_array, y, x) = interpolation(p, p01, p10, p11);
                            is_empty = false;
                        }
                    }
                }
            }
            else if (v01 < minimum_value)
            {
                if (v00 < minimum_value || v10 < minimum_value || v11 < minimum_value) continue;
                int begin_x = (int) ceilf(min3(p00.x, p10.x, p11.x));
                int end_x = (int) floorf(max3(p00.x, p10.x, p11.x));
                int begin_y = (int) ceilf(min3(p00.y, p10.y, p11.y));
                int end_y = (int) floorf(max3(p00.y, p10.y, p11.y));
                if (begin_x < 0) begin_x = 0;
                if (end_x >= width) end_x = width-1;
                if (begin_y < 0) begin_y = 0;
                if (end_y >= height) end_y = height-1;
                for (int y = begin_y; y <= end_y; ++y)
                {
                    for (int x = begin_x; x <= end_x; ++x)
                    {
                        TPoint p = {.x = x, .y = y};
                        if (in_triangle(p, p00, p10, p11))
                        {
                            *(float*)PyArray_GETPTR2(output_array, y, x) = interpolation(p, p00, p10, p11);
                            is_empty = false;
                        }
                    }
                }
            }
            else if (v10 < minimum_value)
            {
                if (v00 < minimum_value || v01 < minimum_value || v11 < minimum_value) continue;
                int begin_x = (int) ceilf(min3(p00.x, p01.x, p11.x));
                int end_x = (int) floorf(max3(p00.x, p01.x, p11.x));
                int begin_y = (int) ceilf(min3(p00.y, p01.y, p11.y));
                int end_y = (int) floorf(max3(p00.y, p01.y, p11.y));
                if (begin_x < 0) begin_x = 0;
                if (end_x >= width) end_x = width-1;
                if (begin_y < 0) begin_y = 0;
                if (end_y >= height) end_y = height-1;
                for (int y = begin_y; y <= end_y; ++y)
                {
                    for (int x = begin_x; x <= end_x; ++x)
                    {
                        TPoint p = {.x = x, .y = y};
                        if (in_triangle(p, p00, p01, p11))
                        {
                            *(float*)PyArray_GETPTR2(output_array, y, x) = interpolation(p, p00, p01, p11);
                            is_empty = false;
                        }
                    }
                }
            }
            else if (v11 < minimum_value)
            {
                if (v00 < minimum_value || v01 < minimum_value || v10 < minimum_value) continue;
                int begin_x = (int) ceilf(min3(p00.x, p01.x, p10.x));
                int end_x = (int) floorf(max3(p00.x, p01.x, p10.x));
                int begin_y = (int) ceilf(min3(p00.x, p01.y, p10.y));
                int end_y = (int) floorf(max3(p00.x, p01.y, p10.y));
                if (begin_x < 0) begin_x = 0;
                if (end_x >= width) end_x = width-1;
                if (begin_y < 0) begin_y = 0;
                if (end_y >= height) end_y = height-1;
                for (int y = begin_y; y <= end_y; ++y)
                {
                    for (int x = begin_x; x <= end_x; ++x)
                    {
                        TPoint p = {.x = x, .y = y};
                        if (in_triangle(p, p00, p01, p10))
                        {
                            *(float*)PyArray_GETPTR2(output_array, y, x) = interpolation(p, p00, p01, p10);
                            is_empty = false;
                        }
                    }
                }
            }
            else
            {
                int begin_x = (int) ceilf(min4(p00.x, p01.x, p10.x, p11.x));
                int end_x = (int) floorf(max4(p00.x, p01.x, p10.x, p11.x));
                int begin_y = (int) ceilf(min4(p00.y, p01.y, p10.y, p11.y));
                int end_y = (int) floorf(max4(p00.y, p01.y, p10.y, p11.y));
                if (begin_x < 0) begin_x = 0;
                if (end_x >= width) end_x = width-1;
                if (begin_y < 0) begin_y = 0;
                if (end_y >= height) end_y = height-1;
                for (int y = begin_y; y <= end_y; ++y)
                {
                    for (int x = begin_x; x <= end_x; ++x)
                    {
                        TPoint p = {.x = x, .y = y};
                        if (in_triangle(p, p00, p01, p10))
                        {
                            *(float*)PyArray_GETPTR2(output_array, y, x) = interpolation(p, p00, p01, p10);
                            is_empty = false;
                        }
                        else if(in_triangle(p, p10, p01, p11))
                        {
                            *(float*)PyArray_GETPTR2(output_array, y, x) = interpolation(p, p10, p01, p11);
                            is_empty = false;
                        }
                    }
                }
            }
        }
    }
    free(new_x_array);
    free(new_y_array);
    if (is_empty) Py_RETURN_NONE;
    return Py_BuildValue("{sfsfsfsfsO}",
                         "min_x", minimum_x,
                         "min_y", minimum_y,
                         "max_x", maximum_x,
                         "max_y", maximum_y,
                         "values", PyArray_Return(output_array));
}
static PyMethodDef
cgriddata_methods[] = {
    {"griddata", griddata, METH_VARARGS, "linear interpolation"},
    {NULL, NULL, 0, NULL}
};
PyMODINIT_FUNC
initcgriddata(void)
{
    Py_InitModule("cgriddata", cgriddata_methods);
    import_array();
}

if (PyArg_ParseTuple(args, "O!O!O!ff",
        &PyArray_Type, &x_array,
        &PyArray_Type, &y_array,
        &PyArray_Type, &values_array,
        &minimum_value, &scale_factor) == 0) return NULL;

PyArrayObject* output_array = PyArray_FromDims(2, new_dims, NPY_FLOAT);

return Py_BuildValue("{sfsfsfsfsO}",
                         "min_x", minimum_x,
                         "min_y", minimum_y,
                         "max_x", maximum_x,
                         "max_y", maximum_y,
                         "values", PyArray_Return(output_array));

PyArrayObject* output_array = PyArray_FromDims(2, new_dims, NPY_FLOAT);

return Py_BuildValue("{sfsfsfsfsO}",
                         "min_x", minimum_x,
                         "min_y", minimum_y,
                         "max_x", maximum_x,
                         "max_y", maximum_y,
                         "values", PyArray_Return(output_array));

PyObject* result = Py_BuildValue("{sfsfsfsfsO}",
                                 "min_x", minimum_x,
                                 "min_y", minimum_y,
                                 "max_x", maximum_x,
                                 "max_y", maximum_y,
                                 "values", PyArray_Return(output_array));
Py_DECREF(output_array);
return result;