/* -*- coding: utf-8;mode:c++;c-file-style:"stroustrup" -*- */

/*

Licensed under the terms of the BSD 3-Clause

(see plotpy/__init__.py for details)

*/

#include <Python.h>

#undef NO_IMPORT_ARRAY

#define PY_ARRAY_UNIQUE_SYMBOL PyScalerArray

#include <numpy/arrayobject.h>

#ifdef _MSC_VER

#include <float.h>

#define FE_TOWARDZERO _RC_CHOP

#define fegetround() (_controlfp(0, 0) & _MCW_RC)

int fesetround(int r)

{

if ((r & _MCW_RC) == r)

{ /* check the supplied value is one of

those allowed */

int result = _controlfp(r, _MCW_RC) & _MCW_RC;

return !(r == result);

}

return !0;

}

#else

#include <fenv.h>

#endif

#include <math.h>

#if defined(_MSC_VER) || defined(__MINGW32__)

#define isnan(x) _isnan(x)

#endif

#include <stdio.h>

#include <algorithm>

#include <iostream>

#include <vector>

#include "points.hpp"

#include "arrays.hpp"

#include "scaler.hpp"

using std::max;

using std::min;

using std::swap;

using std::vector;

enum

{

INTERP_NEAREST = 0,

INTERP_LINEAR = 1,

INTERP_AA = 2

};

typedef union

{

npy_uint32 v;

npy_uint8 c[4];

} rgba_t;

typedef XYTransform<Array1D<double>> XYScale;

template <class Transform>

struct params

{

typedef Transform transform_type;

params(PyArrayObject *_src,

PyArrayObject *_dst, PyObject *_dst_data,

PyObject *_lut, PyObject *_interp,

Transform &_trans) : p_dst(_dst), p_dst_data(_dst_data),

p_src(_src),

p_lut(_lut), p_interpolation(_interp),

trans(_trans)

{

}

// Source, dest, coordinate transformation

PyArrayObject *p_dst;

PyObject *p_dst_data;

PyArrayObject *p_src;

PyObject *p_lut; // Pixel value transformation tuple

PyObject *p_interpolation;

Transform &trans;

int dx1, dx2, dy1, dy2;

};

template <class T, class TR>

struct NearestInterpolation

{

T operator()(const Array2D<T> &src, const TR &tr, const typename TR::point &p)

{

int nx = p.ix();

int ny = p.iy();

return src.value(nx, ny);

}

};

template <class T, class TR>

struct LinearInterpolation

{

T operator()(const Array2D<T> &src, const TR &tr, const typename TR::point &p)

{

int nx = p.ix();

int ny = p.iy();

double v = src.value(nx, ny);

double a = 0;

// The following couple of lines were commented out to avoid disabling

// the linear interpolation on image edges. Demonstrating the effect of

// this change is quite easy: just try to show a very small image

// (e.g. 10x10) with pyplot.imshow for example.

// if (nx==0||nx==src.nj-1) return (T)v;

// if (ny==0||ny==src.ni-1) return (T)v;

if (nx < src.nj - 1)

{

a = p.x() - nx;

v = (1 - a) * v + a * src.value(nx + 1, ny);

}

if (ny >= src.ni - 1)

return (T)v;

double v2 = src.value(nx, ny + 1);

double b = p.y() - ny;

if (nx < src.nj - 1)

{

v2 = (1 - a) * v2 + a * src.value(nx + 1, ny + 1);

}

return (T)(v * (1 - b) + b * v2);

}

};

template <class TR>

struct LinearInterpolation<npy_uint32, TR>

{

npy_uint32 operator()(const Array2D<npy_uint32> &src, const TR &tr, const typename TR::point &p)

{

int k;

int nx = p.ix();

int ny = p.iy();

rgba_t p1, p2, p3, p4, r;

p1.v = src.value(nx, ny);

float v[4], v2[4];

double a = 0;

if (nx < src.nj - 1)

{

p2.v = src.value(nx + 1, ny);

a = p.x() - nx;

for (k = 0; k < 4; ++k)

{

v[k] = (1 - a) * p1.c[k] + a * p2.c[k];

}

}

else

{

for (k = 0; k < 4; ++k)

{

v[k] = p1.c[k];

}

}

if (ny >= src.ni - 1)

{

for (k = 0; k < 4; ++k)

{

r.c[k] = (npy_uint8)(v[k]);

}

return r.v;

}

p3.v = src.value(nx, ny + 1);

double b = p.y() - ny;

if (nx < src.nj - 1)

{

p4.v = src.value(nx + 1, ny + 1);

for (k = 0; k < 4; ++k)

{

v2[k] = (1 - a) * p3.c[k] + a * p4.c[k];

}

}

else

{

for (k = 0; k < 4; ++k)

{

v2[k] = p3.c[k];

}

}

for (k = 0; k < 4; ++k)

{

float px = v[k] * (1 - b) + b * v2[k];

if (px < 0.0)

px = 0.0;

if (px > 255.0)

px = 255.0;

r.c[k] = (npy_uint8)px;

}

return r.v;

}

};

#if 0

template<class T, class TR>

struct QuadInterpolation {

T operator()(const Array2D<T>& src, const TR& tr, const typename TR::point& p) {

int nx = p.ix();

int ny = p.iy();

double v = src.value(nx, ny);

double a=0;

if (nx==0||nx==src.nj-1) return src.value(nx,ny);

if (ny==0||ny==src.ni-1) return src.value(nx,ny);

v0 = interp(src.value(nx-1,ny-1));

if (nx<src.nj-1) {

a = p.x()-nx;

v = (1-a)*v+a*src.value(nx+1,ny);

}

if (ny>=src.ni-1) return v;

double v2 = src.value(nx,ny+1);

double b = p.y()-ny;

if (nx<src.nj-1) {

v2 = (1-a)*v2+a*src.value(nx+1,ny+1);

}

return v*(1-b)+b*v2;

}

};

#endif

template <class T>

struct LinearInterpolation<T, XYScale>

{

T operator()(const Array2D<T> &src, const XYScale &tr, const typename XYScale::point &p)

{

int nx = p.ix();

int ny = p.iy();

double v = src.value(nx, ny);

double a = 0;

if (nx == 0 || nx == src.nj - 1)

return (T)v;

if (ny == 0 || ny == src.ni - 1)

return (T)v;

if (nx < src.nj - 1)

{

double x0 = tr.ax.value(nx);

double x1 = tr.ax.value(nx + 1);

a = (p.x() - x0) / (x1 - x0);

v = (1 - a) * v + a * src.value(nx + 1, ny);

}

if (ny >= src.ni - 1)

return (T)v;

double v2 = src.value(nx, ny + 1);

double y0 = tr.ay.value(ny);

double y1 = tr.ay.value(ny + 1);

double b = (p.y() - y0) / (y1 - y0);

if (nx < src.nj - 1)

{

v2 = (1 - a) * v2 + a * src.value(nx + 1, ny + 1);

}

return (T)(v * (1 - b) + b * v2);

}

};

template <>

struct LinearInterpolation<npy_uint32, XYScale>

{

npy_uint32 operator()(const Array2D<npy_uint32> &src, const XYScale &tr, const XYScale::point &p)

{

return 0;

}

};

template <class T, class TR>

struct SubSampleInterpolation

{

SubSampleInterpolation(const Array2D<T> &_mask) : mask(_mask)

{

ki = 1. / (mask.ni - 1);

kj = 1. / (mask.nj - 1);

}

T operator()(const Array2D<T> &src, const TR &tr, const typename TR::point &p0)

{

int i, j;

typename TR::point p, p1;

typename num_trait<T>::large_type value = 0;

typename num_trait<T>::large_type count = 0, msk, val;

p1.copy(p0);

tr.incy(p1, -0.5);

tr.incx(p1, -0.5);

for (i = 0; i < mask.ni; ++i)

{

p.copy(p1);

for (j = 0; j < mask.nj; ++j)

{

if (p.inside())

{

msk = mask.value(j, i);

val = src.value(p.ix(), p.iy());

value += msk * val;

count += msk;

// printf("i,j=%d,%d : %f,%f / %d,%d = %lfx%lf\n", i, j, p.x(), p.y(), p.ix(), p.iy(), (double)val, (double)msk );

}

tr.incx(p, kj);

}

tr.incy(p1, ki);

}

// printf("%d/%d\n", (int)value, (int)count);

if (count)

return value / count;

else

return value;

}

typename TR::real ki, kj;

const Array2D<T> &mask;

};

template <class DEST, class ST, class Scale, class Trans, class Interpolation>

void _scale_rgb(DEST &dest,

Array2D<ST> &src, const Scale &scale, const Trans &tr,

int dx1, int dy1, int dx2, int dy2,

Interpolation &interpolate)

{

int i, j;

ST val;

int round = fegetround();

PixelIterator<DEST> it(dest);

typename Trans::point p, p0;

fesetround(FE_TOWARDZERO);

/*

printf("SRC: ni=%d nj=%d\n", src.ni, src.nj);

printf("TR: dx=%lf dy=%lf\n", tr.dx, tr.dy);

printf("DST: ni=%d nj=%d si=%d sj=%d\n", dest.ni, dest.nj, dest.si, dest.sj);

*/

tr.set(p0, dx1, dy1);

for (i = dy1; i < dy2; ++i)

{

it.moveto(dx1, i);

p = p0;

for (j = dx1; j < dx2; ++j)

{

if (!p.inside())

{

scale.set_bg(it());

}

else

{

val = interpolate(src, tr, p);

if (isnan(val))

{

scale.set_bg(it());

}

else

{

it() = scale.eval(val);

}

}

tr.incx(p);

it.move(1, 0);

}

tr.incy(p0);

}

fesetround(round);

}

static bool check_dispatch_type(const char *name, PyArrayObject *p_src)

{

if (PyArray_TYPE(p_src) != NPY_DOUBLE &&

PyArray_TYPE(p_src) != NPY_FLOAT &&

PyArray_TYPE(p_src) != NPY_UINT64 &&

PyArray_TYPE(p_src) != NPY_INT64 &&

PyArray_TYPE(p_src) != NPY_UINT32 &&

PyArray_TYPE(p_src) != NPY_INT32 &&

PyArray_TYPE(p_src) != NPY_UINT16 &&

PyArray_TYPE(p_src) != NPY_INT16 &&

PyArray_TYPE(p_src) != NPY_UINT8 &&

PyArray_TYPE(p_src) != NPY_INT8 &&

PyArray_TYPE(p_src) != NPY_BOOL)

{

PyErr_Format(PyExc_TypeError, "%s data type must be one of the following:"

" double, float, uint64, int64, uint32, int32, uint16, int16, uint8, int8, bool",

name);

return false;

}

return true;

}

static bool check_array_2d(const char *name, PyArrayObject *arr, int dtype)

{

if (!PyArray_Check(arr))

{

PyErr_Format(PyExc_TypeError, "%s must be a ndarray", name);

return false;

}

if (PyArray_NDIM(arr) != 2)

{

PyErr_Format(PyExc_TypeError, "%s must be 2-D array", name);

return false;

}

if (dtype >= 0 && PyArray_TYPE(arr) != dtype)

{

PyErr_Format(PyExc_TypeError, "%s data type must be %d", name, dtype);

return false;

}

return true;

}

bool check_arrays(PyArrayObject *p_src, PyArrayObject *p_dest)

{

if (!PyArray_Check(p_src) || !PyArray_Check(p_dest))

{

PyErr_SetString(PyExc_TypeError, "src and dst must be ndarrays");

return false;

}

if (PyArray_TYPE(p_dest) != NPY_UINT32 &&

PyArray_TYPE(p_dest) != NPY_FLOAT32 &&

PyArray_TYPE(p_dest) != NPY_FLOAT64)

{

PyErr_SetString(PyExc_TypeError, "dst data type must be uint32 or float");

return false;

}

if (PyArray_NDIM(p_src) != 2 || PyArray_NDIM(p_dest) != 2)

{

PyErr_SetString(PyExc_TypeError, "dst and src must be 2-D arrays");

return false;

}

return check_dispatch_type("src", p_src);

}

bool check_lut(PyArrayObject *p_lut)

{

if (!PyArray_Check(p_lut))

{

PyErr_SetString(PyExc_TypeError, "lut must be an ndarray");

return false;

}

if (PyArray_NDIM(p_lut) != 1)

{

PyErr_SetString(PyExc_TypeError, "lut must be a 1D array");

return false;

}

if (PyArray_TYPE(p_lut) != NPY_UINT32)

{

PyErr_SetString(PyExc_TypeError, "lut data type must be uint32");

return false;

}

return true;

}

static bool check_transform(PyArrayObject *p_tr)

{

if (!PyArray_Check(p_tr))

{

PyErr_SetString(PyExc_TypeError, "transform must be an ndarray");

return false;

}

if (PyArray_TYPE(p_tr) != NPY_DOUBLE)

{

PyErr_SetString(PyExc_TypeError, "transform data type must be float");

return false;

}

int ni = PyArray_DIM(p_tr, 0);

int nj = PyArray_DIM(p_tr, 1);

if (ni != 3 || nj != 3)

{

PyErr_SetString(PyExc_TypeError, "transform must be 3x3");

return false;

}

return true;

}

static void check_image_bounds(int ni, int nj, int &dx, int &dy)

{

if (dx < 0)

dx = 0;

if (dy < 0)

dy = 0;

if (dx > nj)

dx = nj;

if (dy > ni)

dy = ni;

}

template <class Params, class PixelScale, class Interp>

static bool scale_src_dst_interp(Params &p, PixelScale &pixel_scale, Interp &interp)

{

typedef typename PixelScale::source_type ST;

typedef typename PixelScale::dest_type DT;

Array2D<ST> src(p.p_src);

Array2D<DT> dst(p.p_dst);

_scale_rgb(dst, src, pixel_scale, p.trans,

p.dx1, p.dy1, p.dx2, p.dy2, interp);

return true;

}

template <class Params, class PixelScale>

static bool scale_src_dst(Params &p, PixelScale &pixel_scale)

{

typedef typename PixelScale::source_type ST;

typedef typename Params::transform_type TR;

typedef NearestInterpolation<ST, TR> Nearest;

typedef SubSampleInterpolation<ST, TR> SubAA;

typedef LinearInterpolation<ST, TR> Linear;

int interpolation;

PyArrayObject *p_mask = 0;

if (!PyArg_ParseTuple(p.p_interpolation, "i|O", &interpolation, &p_mask))

{

PyErr_SetString(PyExc_ValueError, "Can't interpret interpolation");

return false;

}

switch (interpolation)

{

case INTERP_NEAREST:

{

Nearest interp;

return scale_src_dst_interp<Params, PixelScale, Nearest>(p, pixel_scale, interp);

}

case INTERP_AA:

{

if (!check_array_2d("AA Mask", p_mask, PyArray_TYPE(p.p_src)))

return false;

Array2D<ST> mask(p_mask);

SubAA interp(mask);

return scale_src_dst_interp<Params, PixelScale, SubAA>(p, pixel_scale, interp);

}

case INTERP_LINEAR:

{

Linear interp;

return scale_src_dst_interp<Params, PixelScale, Linear>(p, pixel_scale, interp);

}

default:

PyErr_SetString(PyExc_ValueError, "Unknown interpolation type");

return false;

};

}

/* we know the transformation and source type, now we dispatch

on the destination type, which determines the LUT transformation

*/

template <class Params, class ST>

static bool scale_src_bw(Params &p)

{

typedef LutScale<ST, npy_uint32> color_scale;

typedef LinearScale<ST, npy_float32> bw32_scale;

typedef LinearScale<ST, npy_float64> bw64_scale;

double a, b;

PyObject *p_bg;

PyArrayObject *p_cmap = 0;

bool apply_bg = true;

if (!PyArg_ParseTuple(p.p_lut, "ddO|O", &a, &b, &p_bg, &p_cmap))

{

PyErr_SetString(PyExc_ValueError, "Can't interpret pixel transformation tuple");

return false;

}

if (p_bg == Py_None)

apply_bg = false;

switch (PyArray_TYPE(p.p_dst))

{

case NPY_UINT32:

{

/* Destination is RGB */

unsigned long bg = 0;

if (apply_bg)

{

#if PY_MAJOR_VERSION >= 3

bg = PyLong_AsUnsignedLongMask(p_bg);

#else

bg = PyInt_AsUnsignedLongMask(p_bg);

#endif

if (PyErr_Occurred())

return false;

}

if (!check_lut(p_cmap))

{

return false;

}

Array1D<npy_uint32> cmap(p_cmap);

color_scale scale(a, b, cmap, bg, apply_bg);

return scale_src_dst<Params, color_scale>(p, scale);

}

case NPY_FLOAT32:

{

double bg = 0.0;

if (apply_bg)

{

bg = PyFloat_AsDouble(p_bg);

if (PyErr_Occurred())

return false;

}

bw32_scale scale(a, b, bg, apply_bg);

return scale_src_dst<Params, bw32_scale>(p, scale);

}

case NPY_FLOAT64:

{

double bg = 0.0;

if (apply_bg)

{

bg = PyFloat_AsDouble(p_bg);

if (PyErr_Occurred())

return false;

}

bw64_scale scale(a, b, bg, apply_bg);

return scale_src_dst<Params, bw64_scale>(p, scale);

}

default:

PyErr_SetString(PyExc_TypeError, "Destination array must be uint32 (rgb) or float (BW)");

return false;

}

}

template <class Params, class ST>

static bool scale_src_rgb(Params &p)

{

// Special case p_lut = None

NoScale<ST, npy_uint32> scale(0, false);

return scale_src_dst<Params, NoScale<ST, npy_uint32>>(p, scale);

}

template <class Params>

static PyObject *dispatch_source(Params &p)

{

bool ok;

int dni = PyArray_DIM(p.p_dst, 0);

int dnj = PyArray_DIM(p.p_dst, 1);

if (!PyArg_ParseTuple(p.p_dst_data, "iiii",

&p.dx1, &p.dy1, &p.dx2, &p.dy2))

{

PyErr_SetString(PyExc_ValueError, "Invalid destination rectangle (expected tuple of 4 integers)");

return NULL;

}

if (p.dx2 < p.dx1)

swap(p.dx1, p.dx2);

if (p.dy2 < p.dy1)

swap(p.dy1, p.dy2);

check_image_bounds(dni, dnj, p.dx1, p.dy1);

check_image_bounds(dni, dnj, p.dx2, p.dy2);

switch (PyArray_TYPE(p.p_src))

{

case NPY_FLOAT32:

ok = scale_src_bw<Params, npy_float32>(p);

break;

case NPY_FLOAT64:

ok = scale_src_bw<Params, npy_float64>(p);

break;

case NPY_UINT64:

ok = scale_src_bw<Params, npy_uint64>(p);

break;

case NPY_INT64:

ok = scale_src_bw<Params, npy_int64>(p);

break;

case NPY_INT32:

ok = scale_src_bw<Params, npy_int32>(p);

break;

case NPY_UINT32: // RGBA

ok = scale_src_rgb<Params, npy_uint32>(p);

break;

case NPY_UINT16:

ok = scale_src_bw<Params, npy_uint16>(p);

break;

case NPY_INT16:

ok = scale_src_bw<Params, npy_int16>(p);

break;

case NPY_UINT8:

ok = scale_src_bw<Params, npy_uint8>(p);

break;

case NPY_INT8:

ok = scale_src_bw<Params, npy_int8>(p);

break;

case NPY_BOOL:

ok = scale_src_bw<Params, npy_bool>(p);

break;

default:

PyErr_SetString(PyExc_TypeError, "Unknown source data type");

return NULL;

}

if (!ok)

return NULL;

return Py_BuildValue("iiii", p.dx1, p.dy1, p.dx2, p.dy2);

}

/* Input data :

SRC, SRC_DATA, DST, DST_DATA, LUT_DATA

SRC : PyArrayObject (i8,u8,i16,u16,float32,float64)

DST : PyArrayObject : u32 -> rgb, float32 : bw

SRC_DATA : varies :

Scale : source rect (x1,y1,x2,y2)

Transform : transformation matrix

XY : source rect, X array, Y array

DST_DATA : dest rect (dx1,dy1,dx2,dy2)

LUT_DATA : (a,b,bg) if DST is bw or (a,b,bg,cmap) if DST is rgb

*/

static PyObject *py_scale_xy(PyObject *self, PyObject *args)

{

typedef params<XYScale> Params;

PyArrayObject *p_src = 0, *p_dst = 0, *p_ax = 0, *p_ay = 0;

PyObject *p_lut_data, *p_src_data, *p_dst_data, *p_interp_data;

double x1, y1, x2, y2;

if (!PyArg_ParseTuple(args, "OOOOOO:_scale_xy",

&p_src, &p_src_data,

&p_dst, &p_dst_data,

&p_lut_data, &p_interp_data))

{

return NULL;

}

if (!check_arrays(p_src, p_dst))

{

return NULL;

}

if (!PyArg_ParseTuple(p_src_data, "OO(dddd):_scale_xy",

&p_ax, &p_ay, &x1, &y1, &x2, &y2))

{

return NULL;

}

int ni = PyArray_DIM(p_src, 0);

int nj = PyArray_DIM(p_src, 1);

int dni = PyArray_DIM(p_dst, 0);

int dnj = PyArray_DIM(p_dst, 1);

double dx = (x2 - x1) / dnj;

double dy = (y2 - y1) / dni;

Array1D<double> ax(p_ax), ay(p_ay);

XYScale trans(nj, ni, ax, ay, x1, y1, dx, dy);

Params scale_params(p_src, p_dst, p_dst_data,

p_lut_data, p_interp_data, trans);

// examine source type

return dispatch_source<Params>(scale_params);

}

static PyObject *py_scale_tr(PyObject *self, PyObject *args)

{

typedef params<LinearTransform> Params;

PyArrayObject *p_src = 0, *p_dst = 0, *p_tr;

PyObject *p_lut_data, *p_dst_data, *p_interp_data;

if (!PyArg_ParseTuple(args, "OOOOOO:_scale_tr",

&p_src, &p_tr,

&p_dst, &p_dst_data,

&p_lut_data, &p_interp_data))

{

return NULL;

}

if (!check_arrays(p_src, p_dst))

{

return NULL;

}

if (!check_transform(p_tr))

{

return NULL;

}

int ni = PyArray_DIM(p_src, 0);

int nj = PyArray_DIM(p_src, 1);

Array2D<double> tr(p_tr);

LinearTransform trans(nj, ni,

tr.value(2, 0), tr.value(2, 1), // x0, y0

tr.value(0, 0), tr.value(1, 0), // xx, xy

tr.value(0, 1), tr.value(1, 1) // yx, yy

);

Params scale_params(p_src, p_dst, p_dst_data,

p_lut_data, p_interp_data, trans);

// examine source type

return dispatch_source<Params>(scale_params);

}

static PyObject *py_scale_rect(PyObject *self, PyObject *args)

{

typedef params<ScaleTransform> Params;

PyArrayObject *p_src = 0, *p_dst = 0;

PyObject *p_lut_data, *p_dst_data, *p_interp_data, *p_src_data;

double x1, x2, y1, y2;

if (!PyArg_ParseTuple(args, "OOOOOO:_scale_rect",

&p_src, &p_src_data,

&p_dst, &p_dst_data,

&p_lut_data, &p_interp_data))

{

return NULL;

}

if (!check_arrays(p_src, p_dst))

{

return NULL;

}

if (!PyArg_ParseTuple(p_src_data, "dddd:_scale_rect",

&x1, &y1, &x2, &y2))

{

return NULL;

}

int ni = PyArray_DIM(p_src, 0);

int nj = PyArray_DIM(p_src, 1);

int dni = PyArray_DIM(p_dst, 0);

int dnj = PyArray_DIM(p_dst, 1);

double dx = (x2 - x1) / dnj;

double dy = (y2 - y1) / dni;

ScaleTransform trans(nj, ni, x1, y1, dx, dy);

Params scale_params(p_src, p_dst, p_dst_data,

p_lut_data, p_interp_data, trans);

// examine source type

return dispatch_source<Params>(scale_params);

}

class Histogram

{

public:

Histogram(PyArrayObject *_data, PyArrayObject *_bins,

PyArrayObject *_res, int ignore_bounds) : p_data(_data),

p_bins(_bins),

p_res(_res),

mode(ignore_bounds)

{

}

template <class T>

void run()

{

Array1D<npy_uint32> res(p_res);

Array1D<T> data(p_data);

Array1D<T> bins(p_bins);

typename Array1D<T>::iterator it, end, bin, bin0, bin1;

end = data.end();

bin0 = bins.begin();

bin1 = bins.end();

for (it = data.begin(); it < end; ++it)

{

bin = std::lower_bound(bin0, bin1, *it);

res.value(bin - bin0)++;

}

}

PyArrayObject *p_data, *p_bins, *p_res;

int mode;

};

static PyObject *py_histogram(PyObject *self, PyObject *args)

{

PyArrayObject *p_data = 0, *p_bins = 0, *p_res = 0;

int ignore_bounds = 0;

if (!PyArg_ParseTuple(args, "OOO|i:_histogram", &p_data, &p_bins, &p_res,

&ignore_bounds))

{

return NULL;

}

if (!PyArray_Check(p_data) ||

!PyArray_Check(p_bins) ||

!PyArray_Check(p_res))

{

PyErr_SetString(PyExc_TypeError, "data, bins, dest must be ndarray");

return NULL;

}

if (!check_dispatch_type("data", p_data))

{

return NULL;

}

Histogram hist(p_data, p_bins, p_res, ignore_bounds);

dispatch_array(PyArray_TYPE(p_data), hist);

Py_INCREF(Py_None);

return Py_None;

}

PyObject *py_vert_line(PyObject *self, PyObject *args);

PyObject *py_scale_quads(PyObject *self, PyObject *args);

static PyMethodDef _meths[] = {

{"_scale_xy", py_scale_xy, METH_VARARGS,

"transform source to destination with arbitrary X and Y scale"},

{"_scale_tr", py_scale_tr, METH_VARARGS,

"Linear transformation of source to destination"},

{"_scale_rect", py_scale_rect, METH_VARARGS,

"Linear rescale of source to destination parallel to axes"},

{"_scale_quads", py_scale_quads, METH_VARARGS,

"Linear rescale of a structured grid to destination parallel to axes"},

{"_histogram", py_histogram, METH_VARARGS,

"Compute histogram of 1d data"},

{"_line_test", py_vert_line, METH_VARARGS,

"Rasterize lines"},

{NULL, NULL, 0, NULL} /* Sentinel */

};

#if PY_MAJOR_VERSION >= 3

static struct PyModuleDef moduledef = {

PyModuleDef_HEAD_INIT,

"_scaler", /* m_name */

"Scaler module", /* m_doc */

-1, /* m_size */

_meths, /* m_methods */

NULL, /* m_reload */

NULL, /* m_traverse */

NULL, /* m_clear */

NULL, /* m_free */

};

#endif

PyMODINIT_FUNC

#if PY_MAJOR_VERSION >= 3

PyInit__scaler(void)

#else

init_scaler()

#endif

{

PyObject *m;

#if PY_MAJOR_VERSION >= 3

m = PyModule_Create(&moduledef);

#else

m = Py_InitModule("_scaler", _meths);

#endif

import_array();

PyModule_AddIntConstant(m, "INTERP_NEAREST", INTERP_NEAREST);

PyModule_AddIntConstant(m, "INTERP_LINEAR", INTERP_LINEAR);

PyModule_AddIntConstant(m, "INTERP_AA", INTERP_AA);

#if PY_MAJOR_VERSION >= 3

return m;

#endif

}