//.........这里部分代码省略.........
        part = PyArray_EnsureAnyArray(part);
        round_part = PyArray_Round((PyArrayObject *)part,
                                   decimals, NULL);
        Py_DECREF(part);
        if (round_part == NULL) {
            Py_DECREF(arr);
            return NULL;
        }
        res = PyObject_SetAttrString(arr, "imag", round_part);
        Py_DECREF(round_part);
        if (res < 0) {
            Py_DECREF(arr);
            return NULL;
        }
        return arr;
    }
    /* do the most common case first */
    if (decimals >= 0) {
        if (PyArray_ISINTEGER(a)) {
            if (out) {
                if (PyArray_AssignArray(out, a,
                            NULL, NPY_DEFAULT_ASSIGN_CASTING) < 0) {
                    return NULL;
                }
                Py_INCREF(out);
                return (PyObject *)out;
            }
            else {
                Py_INCREF(a);
                return (PyObject *)a;
            }
        }
        if (decimals == 0) {
            if (out) {
                return PyObject_CallFunction(n_ops.rint, "OO", a, out);
            }
            return PyObject_CallFunction(n_ops.rint, "O", a);
        }
        op1 = n_ops.multiply;
        op2 = n_ops.true_divide;
    }
    else {
        op1 = n_ops.true_divide;
        op2 = n_ops.multiply;
        decimals = -decimals;
    }
    if (!out) {
        if (PyArray_ISINTEGER(a)) {
            ret_int = 1;
            my_descr = PyArray_DescrFromType(NPY_DOUBLE);
        }
        else {
            Py_INCREF(PyArray_DESCR(a));
            my_descr = PyArray_DESCR(a);
        }
        out = (PyArrayObject *)PyArray_Empty(PyArray_NDIM(a), PyArray_DIMS(a),
                                             my_descr,
                                             PyArray_ISFORTRAN(a));
        if (out == NULL) {
            return NULL;
        }
    }
    else {
        Py_INCREF(out);
    }
    f = PyFloat_FromDouble(power_of_ten(decimals));
    if (f == NULL) {
        return NULL;
    }
    ret = PyObject_CallFunction(op1, "OOO", a, f, out);
    if (ret == NULL) {
        goto finish;
    }
    tmp = PyObject_CallFunction(n_ops.rint, "OO", ret, ret);
    if (tmp == NULL) {
        Py_DECREF(ret);
        ret = NULL;
        goto finish;
    }
    Py_DECREF(tmp);
    tmp = PyObject_CallFunction(op2, "OOO", ret, f, ret);
    if (tmp == NULL) {
        Py_DECREF(ret);
        ret = NULL;
        goto finish;
    }
    Py_DECREF(tmp);

 finish:
    Py_DECREF(f);
    Py_DECREF(out);
    if (ret_int) {
        Py_INCREF(PyArray_DESCR(a));
        tmp = PyArray_CastToType((PyArrayObject *)ret,
                                 PyArray_DESCR(a), PyArray_ISFORTRAN(a));
        Py_DECREF(ret);
        return tmp;
    }
    return ret;
}

//.........这里部分代码省略.........
            Py_DECREF(dtype);
            return RETRY_WITH_UNICODE;
        }
        *out_dtype = dtype;
        return 0;
    }
    /* Do type promotion with 'out_dtype' */
    else {
        PyArray_Descr *res_dtype = PyArray_PromoteTypes(dtype, *out_dtype);
        Py_DECREF(dtype);
        if (res_dtype == NULL) {
            return -1;
        }
        if (!string_type &&
                res_dtype->type_num == NPY_UNICODE &&
                (*out_dtype)->type_num != NPY_UNICODE) {
            Py_DECREF(res_dtype);
            return RETRY_WITH_UNICODE;
        }
        if (!string_type &&
                res_dtype->type_num == NPY_STRING &&
                (*out_dtype)->type_num != NPY_STRING) {
            Py_DECREF(res_dtype);
            return RETRY_WITH_STRING;
        }
        Py_DECREF(*out_dtype);
        *out_dtype = res_dtype;
        return 0;
    }

fail:
    Py_XDECREF(*out_dtype);
    *out_dtype = NULL;
    return -1;
}

#undef RETRY_WITH_STRING
#undef RETRY_WITH_UNICODE

/* new reference */
NPY_NO_EXPORT PyArray_Descr *
_array_typedescr_fromstr(char *c_str)
{
    PyArray_Descr *descr = NULL;
    PyObject *stringobj = PyString_FromString(c_str);

    if (stringobj == NULL) {
        return NULL;
    }
    if (PyArray_DescrConverter(stringobj, &descr) != NPY_SUCCEED) {
        Py_DECREF(stringobj);
        return NULL;
    }
    Py_DECREF(stringobj);
    return descr;
}

NPY_NO_EXPORT int
check_and_adjust_index(npy_intp *index, npy_intp max_item, int axis)
{
    /* Check that index is valid, taking into account negative indices */
    if ((*index < -max_item) || (*index >= max_item)) {
        /* Try to be as clear as possible about what went wrong. */
        if (axis >= 0) {
            PyErr_Format(PyExc_IndexError,
                         "index %"NPY_INTP_FMT" is out of bounds "
                         "for axis %d with size %"NPY_INTP_FMT,
                         *index, axis, max_item);
        } else {
            PyErr_Format(PyExc_IndexError,
                         "index %"NPY_INTP_FMT" is out of bounds "
                         "for size %"NPY_INTP_FMT,
                         *index, max_item);
        }
        return -1;
    }
    /* adjust negative indices */
    if (*index < 0) {
        *index += max_item;
    }
    return 0;
}

NPY_NO_EXPORT char *
index2ptr(PyArrayObject *mp, npy_intp i)
{
    npy_intp dim0;

    if (PyArray_NDIM(mp) == 0) {
        PyErr_SetString(PyExc_IndexError, "0-d arrays can't be indexed");
        return NULL;
    }
    dim0 = PyArray_DIMS(mp)[0];
    if (check_and_adjust_index(&i, dim0, 0) < 0)
        return NULL;
    if (i == 0) {
        return PyArray_DATA(mp);
    }
    return PyArray_BYTES(mp)+i*PyArray_STRIDES(mp)[0];
}

/*
 * Conforms an output parameter 'out' to have 'ndim' dimensions
 * with dimensions of size one added in the appropriate places
 * indicated by 'axis_flags'.
 *
 * The return value is a view into 'out'.
 */
static PyArrayObject *
conform_reduce_result(int ndim, npy_bool *axis_flags,
                    PyArrayObject *out, int keepdims, const char *funcname)
{
    npy_intp strides[NPY_MAXDIMS], shape[NPY_MAXDIMS];
    npy_intp *strides_out = PyArray_STRIDES(out);
    npy_intp *shape_out = PyArray_DIMS(out);
    int idim, idim_out, ndim_out = PyArray_NDIM(out);
    PyArray_Descr *dtype;
    PyArrayObject_fields *ret;

    /*
     * If the 'keepdims' parameter is true, do a simpler validation and
     * return a new reference to 'out'.
     */
    if (keepdims) {
        if (PyArray_NDIM(out) != ndim) {
            PyErr_Format(PyExc_ValueError,
                    "output parameter for reduction operation %s "
                    "has the wrong number of dimensions (must match "
                    "the operand's when keepdims=True)", funcname);
            return NULL;
        }

        for (idim = 0; idim < ndim; ++idim) {
            if (axis_flags[idim]) {
                if (shape_out[idim] != 1) {
                    PyErr_Format(PyExc_ValueError,
                            "output parameter for reduction operation %s "
                            "has a reduction dimension not equal to one "
                            "(required when keepdims=True)", funcname);
                    return NULL;
                }
            }
        }

        Py_INCREF(out);
        return out;
    }

    /* Construct the strides and shape */
    idim_out = 0;
    for (idim = 0; idim < ndim; ++idim) {
        if (axis_flags[idim]) {
            strides[idim] = 0;
            shape[idim] = 1;
        }
        else {
            if (idim_out >= ndim_out) {
                PyErr_Format(PyExc_ValueError,
                        "output parameter for reduction operation %s "
                        "does not have enough dimensions", funcname);
                return NULL;
            }
            strides[idim] = strides_out[idim_out];
            shape[idim] = shape_out[idim_out];
            ++idim_out;
        }
    }

    if (idim_out != ndim_out) {
        PyErr_Format(PyExc_ValueError,
                "output parameter for reduction operation %s "
                "has too many dimensions", funcname);
        return NULL;
    }

    /* Allocate the view */
    dtype = PyArray_DESCR(out);
    Py_INCREF(dtype);
    ret = (PyArrayObject_fields *)PyArray_NewFromDescr(&PyArray_Type,
                               dtype,
                               ndim, shape,
                               strides,
                               PyArray_DATA(out),
                               PyArray_FLAGS(out),
                               NULL);
    if (ret == NULL) {
        return NULL;
    }
    Py_INCREF(out);
    if (PyArray_SetBaseObject((PyArrayObject *)ret, (PyObject *)out) < 0) {
        Py_DECREF(ret);
        return NULL;
    }

    return (PyArrayObject *)ret;
}

static PyObject* gist_extract(PyObject *self, PyObject *args)
{
	int nblocks=4;
	int n_scale=3;
	int orientations_per_scale[50]={8,8,4};
	PyArrayObject *image, *descriptor;

	if (!PyArg_ParseTuple(args, "O", &image))
	{
		return NULL;
	}

	if (PyArray_TYPE(image) != NPY_UINT8) {
		PyErr_SetString(PyExc_TypeError, "type of image must be uint8");
		return NULL;
	}

	if (PyArray_NDIM(image) != 3) {
		PyErr_SetString(PyExc_TypeError, "dimensions of image must be 3.");
		return NULL;
	}

	npy_intp *dims_image = PyArray_DIMS(image);


	const int w = (int) *(dims_image+1);
	const int h = (int) *(dims_image);

	// Read image to color_image_t structure
	color_image_t *im=color_image_new(w,h);

	for (int y=0, i=0 ; y<h ; ++y) {
		for (int x=0 ; x<w ; ++x, ++i) {
			im->c1[i] = *(unsigned char *)PyArray_GETPTR3(image, y, x, 0);
			im->c2[i] = *(unsigned char *)PyArray_GETPTR3(image, y, x, 1);
			im->c3[i] = *(unsigned char *)PyArray_GETPTR3(image, y, x, 2);
		}
	}

	// Extract descriptor
	float *desc=color_gist_scaletab(im,nblocks,n_scale,orientations_per_scale);

	int descsize=0;
	/* compute descriptor size */
	for(int i=0;i<n_scale;i++)
		descsize+=nblocks*nblocks*orientations_per_scale[i];

	descsize*=3; /* color */


	// Allocate output
	npy_intp dim_desc[1] = {descsize};
	descriptor = (PyArrayObject *) PyArray_SimpleNew(1, dim_desc, NPY_FLOAT);

	// Set val
	for (int i=0 ; i<descsize ; ++i) {
		*(float *)PyArray_GETPTR1(descriptor, i) = desc[i];
	}

	// Release memory
	color_image_delete(im);
	free(desc);

	return PyArray_Return(descriptor);
}

Mat fromNDArrayToMat(PyObject* o) {
	cv::Mat m;
	bool allowND = true;
	if (!PyArray_Check(o)) {
		failmsg("argument is not a numpy array");
		if (!m.data)
			m.allocator = &g_numpyAllocator;
	} else {
		PyArrayObject* oarr = (PyArrayObject*) o;

		bool needcopy = false, needcast = false;
		int typenum = PyArray_TYPE(oarr), new_typenum = typenum;
		int type = typenum == NPY_UBYTE ? CV_8U : typenum == NPY_BYTE ? CV_8S :
					typenum == NPY_USHORT ? CV_16U :
					typenum == NPY_SHORT ? CV_16S :
					typenum == NPY_INT ? CV_32S :
					typenum == NPY_INT32 ? CV_32S :
					typenum == NPY_FLOAT ? CV_32F :
					typenum == NPY_DOUBLE ? CV_64F : -1;

		if (type < 0) {
			if (typenum == NPY_INT64 || typenum == NPY_UINT64
					|| type == NPY_LONG) {
				needcopy = needcast = true;
				new_typenum = NPY_INT;
				type = CV_32S;
			} else {
				failmsg("Argument data type is not supported");
				m.allocator = &g_numpyAllocator;
				return m;
			}
		}

#ifndef CV_MAX_DIM
		const int CV_MAX_DIM = 32;
#endif

		int ndims = PyArray_NDIM(oarr);
		if (ndims >= CV_MAX_DIM) {
			failmsg("Dimensionality of argument is too high");
			if (!m.data)
				m.allocator = &g_numpyAllocator;
			return m;
		}

		int size[CV_MAX_DIM + 1];
		size_t step[CV_MAX_DIM + 1];
		size_t elemsize = CV_ELEM_SIZE1(type);
		const npy_intp* _sizes = PyArray_DIMS(oarr);
		const npy_intp* _strides = PyArray_STRIDES(oarr);
		bool ismultichannel = ndims == 3 && _sizes[2] <= CV_CN_MAX;

		for (int i = ndims - 1; i >= 0 && !needcopy; i--) {
			// these checks handle cases of
			//  a) multi-dimensional (ndims > 2) arrays, as well as simpler 1- and 2-dimensional cases
			//  b) transposed arrays, where _strides[] elements go in non-descending order
			//  c) flipped arrays, where some of _strides[] elements are negative
			if ((i == ndims - 1 && (size_t) _strides[i] != elemsize)
					|| (i < ndims - 1 && _strides[i] < _strides[i + 1]))
				needcopy = true;
		}

		if (ismultichannel && _strides[1] != (npy_intp) elemsize * _sizes[2])
			needcopy = true;

		if (needcopy) {

			if (needcast) {
				o = PyArray_Cast(oarr, new_typenum);
				oarr = (PyArrayObject*) o;
			} else {
				oarr = PyArray_GETCONTIGUOUS(oarr);
				o = (PyObject*) oarr;
			}

			_strides = PyArray_STRIDES(oarr);
		}

		for (int i = 0; i < ndims; i++) {
			size[i] = (int) _sizes[i];
			step[i] = (size_t) _strides[i];
		}

		// handle degenerate case
		if (ndims == 0) {
			size[ndims] = 1;
			step[ndims] = elemsize;
			ndims++;
		}

		if (ismultichannel) {
			ndims--;
			type |= CV_MAKETYPE(0, size[2]);
		}

		if (ndims > 2 && !allowND) {
			failmsg("%s has more than 2 dimensions");
		} else {

			m = Mat(ndims, size, type, PyArray_DATA(oarr), step);
//.........这里部分代码省略.........

/*@[email protected]*/ static PyObject*
Wcs_all_pix2world(
    Wcs* self,
    PyObject* args,
    PyObject* kwds) {

    int            naxis      = 2;
    PyObject*      pixcrd_obj = NULL;
    int            origin     = 1;
    PyArrayObject* pixcrd     = NULL;
    PyArrayObject* world      = NULL;
    int            status     = -1;
    const char*    keywords[] = {
        "pixcrd", "origin", NULL
    };

    if (!PyArg_ParseTupleAndKeywords(
                args, kwds, "Oi:all_pix2world", (char **)keywords,
                &pixcrd_obj, &origin)) {
        return NULL;
    }

    naxis = self->x.wcs->naxis;

    pixcrd = (PyArrayObject*)PyArray_ContiguousFromAny(pixcrd_obj, PyArray_DOUBLE, 2, 2);
    if (pixcrd == NULL) {
        return NULL;
    }

    if (PyArray_DIM(pixcrd, 1) < naxis) {
        PyErr_Format(
            PyExc_RuntimeError,
            "Input array must be 2-dimensional, where the second dimension >= %d",
            naxis);
        goto exit;
    }

    world = (PyArrayObject*)PyArray_SimpleNew(2, PyArray_DIMS(pixcrd), PyArray_DOUBLE);
    if (world == NULL) {
        goto exit;
    }

    /* Make the call */
    Py_BEGIN_ALLOW_THREADS
    preoffset_array(pixcrd, origin);
    wcsprm_python2c(self->x.wcs);
    status = pipeline_all_pixel2world(&self->x,
                                      (unsigned int)PyArray_DIM(pixcrd, 0),
                                      (unsigned int)PyArray_DIM(pixcrd, 1),
                                      (double*)PyArray_DATA(pixcrd),
                                      (double*)PyArray_DATA(world));
    wcsprm_c2python(self->x.wcs);
    unoffset_array(pixcrd, origin);
    Py_END_ALLOW_THREADS
    /* unoffset_array(world, origin); */

exit:
    Py_XDECREF(pixcrd);

    if (status == 0 || status == 8) {
        return (PyObject*)world;
    } else if (status == -1) {
        PyErr_SetString(
            PyExc_ValueError,
            "Wrong number of dimensions in input array.  Expected 2.");
        return NULL;
    } else {
        Py_DECREF(world);
        if (status == -1) {
            /* exception already set */
            return NULL;
        } else {
            wcserr_to_python_exc(self->x.err);
            return NULL;
        }
    }
}

/*@[email protected]*/ static PyObject*
Wcs_pix2foc(
    Wcs* self,
    PyObject* args,
    PyObject* kwds) {

    PyObject*      pixcrd_obj = NULL;
    int            origin     = 1;
    PyArrayObject* pixcrd     = NULL;
    PyArrayObject* foccrd     = NULL;
    int            status     = -1;
    const char*    keywords[] = {
        "pixcrd", "origin", NULL
    };

    if (!PyArg_ParseTupleAndKeywords(args, kwds, "Oi:pix2foc", (char **)keywords,
                                     &pixcrd_obj, &origin)) {
        return NULL;
    }

    pixcrd = (PyArrayObject*)PyArray_ContiguousFromAny(pixcrd_obj, PyArray_DOUBLE, 2, 2);
    if (pixcrd == NULL) {
        return NULL;
    }

    if (PyArray_DIM(pixcrd, 1) != NAXES) {
        PyErr_SetString(PyExc_ValueError, "Pixel array must be an Nx2 array");
        goto _exit;
    }

    foccrd = (PyArrayObject*)PyArray_SimpleNew(2, PyArray_DIMS(pixcrd), PyArray_DOUBLE);
    if (foccrd == NULL) {
        goto _exit;
    }

    Py_BEGIN_ALLOW_THREADS
    preoffset_array(pixcrd, origin);
    status = pipeline_pix2foc(&self->x,
                              (unsigned int)PyArray_DIM(pixcrd, 0),
                              (unsigned int)PyArray_DIM(pixcrd, 1),
                              (double*)PyArray_DATA(pixcrd),
                              (double*)PyArray_DATA(foccrd));
    unoffset_array(pixcrd, origin);
    unoffset_array(foccrd, origin);
    Py_END_ALLOW_THREADS

_exit:

    Py_XDECREF(pixcrd);

    if (status == 0) {
        return (PyObject*)foccrd;
    } else {
        Py_XDECREF(foccrd);
        if (status == -1) {
            /* Exception already set */
            return NULL;
        } else {
            wcserr_to_python_exc(self->x.err);
            return NULL;
        }
    }
}

//.........这里部分代码省略.........
        }

        /* here we have always intent(in) or intent(inout) or intent(inplace) */

        if (check_and_fix_dimensions(arr,rank,dims)) {
            return NULL; /*XXX: set exception */
        }
	/*
	printf("intent alignement=%d\n", F2PY_GET_ALIGNMENT(intent));
	printf("alignement check=%d\n", F2PY_CHECK_ALIGNMENT(arr, intent));
	int i;
	for (i=1;i<=16;i++)
	  printf("i=%d isaligned=%d\n", i, ARRAY_ISALIGNED(arr, i));
	*/
        if ((! (intent & F2PY_INTENT_COPY))
            && PyArray_ITEMSIZE(arr)==elsize
            && ARRAY_ISCOMPATIBLE(arr,type_num)
	    && F2PY_CHECK_ALIGNMENT(arr, intent)
            ) {
            if ((intent & F2PY_INTENT_C)?PyArray_ISCARRAY(arr):PyArray_ISFARRAY(arr)) {
                if ((intent & F2PY_INTENT_OUT)) {
                    Py_INCREF(arr);
                }
                /* Returning input array */
                return arr;
            }
        }

        if (intent & F2PY_INTENT_INOUT) {
            strcpy(mess, "failed to initialize intent(inout) array");
            if ((intent & F2PY_INTENT_C) && !PyArray_ISCARRAY(arr))
                strcat(mess, " -- input not contiguous");
            if (!(intent & F2PY_INTENT_C) && !PyArray_ISFARRAY(arr))
                strcat(mess, " -- input not fortran contiguous");
            if (PyArray_ITEMSIZE(arr)!=elsize)
                sprintf(mess+strlen(mess),
                        " -- expected elsize=%d but got %" NPY_INTP_FMT,
                        elsize,
                        (npy_intp)PyArray_ITEMSIZE(arr)
                        );
            if (!(ARRAY_ISCOMPATIBLE(arr,type_num)))
                sprintf(mess+strlen(mess)," -- input '%c' not compatible to '%c'",
                        PyArray_DESCR(arr)->type,typechar);
	    if (!(F2PY_CHECK_ALIGNMENT(arr, intent)))
	      sprintf(mess+strlen(mess)," -- input not %d-aligned", F2PY_GET_ALIGNMENT(intent));
            PyErr_SetString(PyExc_ValueError,mess);
            return NULL;
        }

        /* here we have always intent(in) or intent(inplace) */

        {
            PyArrayObject *retarr = (PyArrayObject *) \
                PyArray_New(&PyArray_Type, PyArray_NDIM(arr), PyArray_DIMS(arr), type_num,
                            NULL,NULL,0,
                            !(intent&F2PY_INTENT_C),
                            NULL);
            if (retarr==NULL)
                return NULL;
            F2PY_REPORT_ON_ARRAY_COPY_FROMARR;
            if (PyArray_CopyInto(retarr, arr)) {
                Py_DECREF(retarr);
                return NULL;
            }
            if (intent & F2PY_INTENT_INPLACE) {
                if (swap_arrays(arr,retarr))
                    return NULL; /* XXX: set exception */
                Py_XDECREF(retarr);
                if (intent & F2PY_INTENT_OUT)
                    Py_INCREF(arr);
            } else {
                arr = retarr;
            }
        }
        return arr;
    }

    if ((intent & F2PY_INTENT_INOUT) ||
            (intent & F2PY_INTENT_INPLACE) ||
            (intent & F2PY_INTENT_CACHE)) {
        PyErr_SetString(PyExc_TypeError,
                        "failed to initialize intent(inout|inplace|cache) "
                        "array, input not an array");
        return NULL;
    }

    {
        F2PY_REPORT_ON_ARRAY_COPY_FROMANY;
        arr = (PyArrayObject *) \
            PyArray_FromAny(obj,PyArray_DescrFromType(type_num), 0,0,
                            ((intent & F2PY_INTENT_C)?NPY_ARRAY_CARRAY:NPY_ARRAY_FARRAY) \
                            | NPY_ARRAY_FORCECAST, NULL);
        if (arr==NULL)
            return NULL;
        if (check_and_fix_dimensions(arr,rank,dims))
            return NULL; /*XXX: set exception */
        return arr;
    }

}

/*@[email protected]*/ static PyObject*
PySip_foc2pix(
    PySip* self,
    PyObject* args,
    PyObject* kwds) {

  PyObject*      foccrd_obj = NULL;
  int            origin     = 1;
  PyArrayObject* foccrd     = NULL;
  PyArrayObject* pixcrd     = NULL;
  int            status     = -1;
  const char*    keywords[] = {
    "foccrd", "origin", NULL };

  if (!PyArg_ParseTupleAndKeywords(args, kwds, "Oi:foc2pix", (char **)keywords,
                                   &foccrd_obj, &origin)) {
    return NULL;
  }

  if (self->x.ap == NULL || self->x.bp == NULL) {
    PyErr_SetString(
        PyExc_ValueError,
        "SIP object does not have coefficients for foc2pix transformation (AP and BP)");
    return NULL;
  }

  foccrd = (PyArrayObject*)PyArray_ContiguousFromAny(foccrd_obj, PyArray_DOUBLE, 2, 2);
  if (foccrd == NULL) {
    goto exit;
  }

  if (PyArray_DIM(foccrd, 1) != 2) {
    PyErr_SetString(PyExc_ValueError, "Pixel array must be an Nx2 array");
    goto exit;
  }

  pixcrd = (PyArrayObject*)PyArray_SimpleNew(2, PyArray_DIMS(foccrd),
                                             PyArray_DOUBLE);
  if (pixcrd == NULL) {
    status = 2;
    goto exit;
  }

  Py_BEGIN_ALLOW_THREADS
  preoffset_array(foccrd, origin);
  status = sip_foc2pix(&self->x,
                       (unsigned int)PyArray_DIM(pixcrd, 1),
                       (unsigned int)PyArray_DIM(pixcrd, 0),
                       (double*)PyArray_DATA(foccrd),
                       (double*)PyArray_DATA(pixcrd));
  unoffset_array(foccrd, origin);
  unoffset_array(pixcrd, origin);
  Py_END_ALLOW_THREADS

 exit:
  Py_XDECREF(foccrd);

  if (status == 0) {
    return (PyObject*)pixcrd;
  } else {
    Py_XDECREF(pixcrd);
    if (status == -1) {
      /* Exception already set */
      return NULL;
    } else {
      wcserr_to_python_exc(self->x.err);
      return NULL;
    }
  }
}

int APPLY_SPECIFIC(ctc_cost_cpu)(PyArrayObject *  in_activations,
                                 PyArrayObject *  in_labels,
                                 PyArrayObject *  in_input_lengths,
                                 PyArrayObject ** out_costs,
                                 PyArrayObject ** out_gradients)
{
    ctc_context_t ctc_object;
    ctc_context_t * context = &ctc_object;
    ctc_context_init( context );

    if ( !PyArray_IS_C_CONTIGUOUS( in_activations ) )
    {
        PyErr_SetString( PyExc_RuntimeError,
            "ConnectionistTemporalClassification: activations array must be C-contiguous." );
        return 1;
    }

    npy_float32 * activations = (npy_float32 *) PyArray_DATA( in_activations );

    create_contiguous_input_lengths( in_input_lengths, &(context->input_lengths) );

    if ( NULL == context->input_lengths )
    {
        // Destroy previous CTC context before returning exception
        ctc_context_destroy( context );

        PyErr_Format( PyExc_MemoryError,
            "ConnectionistTemporalClassification: Could not allocate memory for input lengths" );
        return 1;
    }

    // flatten labels to conform with library memory layout
    create_flat_labels( in_labels, &(context->flat_labels), &(context->label_lengths) );

    if ( ( NULL == context->label_lengths ) || ( NULL == context->flat_labels ) )
    {
        // Destroy previous CTC context before returning exception
        ctc_context_destroy( context );

        PyErr_Format( PyExc_MemoryError,
            "ConnectionistTemporalClassification: Could not allocate memory for labels and their lengths" );
        return 1;
    }

    npy_int minibatch_size = PyArray_DIMS( in_activations )[1];
    npy_int alphabet_size = PyArray_DIMS( in_activations )[2];

    npy_float32 * costs = NULL;
    npy_intp cost_size = minibatch_size;

    if ( (*out_costs) == NULL ||                       // Symbolic variable has no memory backing
         PyArray_NDIM( *out_costs ) != 1 ||            // or, matrix has the wrong size
         PyArray_DIMS( *out_costs )[0] != cost_size )
    {
        Py_XDECREF( *out_costs );
        // Allocate new matrix
        *out_costs = (PyArrayObject *) PyArray_ZEROS( 1, &cost_size, NPY_FLOAT32, 0 );

        if ( NULL == (*out_costs) )
        {
            // Destroy previous CTC context before returning exception
            ctc_context_destroy( context );

            PyErr_Format( PyExc_MemoryError,
                "ConnectionistTemporalClassification: Could not allocate memory for CTC costs" );
            return 1;
        }
    }

    costs = (npy_float32 *) PyArray_DATA( *out_costs );

    npy_float32 * gradients = NULL;

    if ( NULL != out_gradients )  // If gradient computation is not disabled
    {
        if ( NULL == (*out_gradients) ||  // Symbolic variable has no real backing
            PyArray_NDIM( *out_gradients ) != 3 ||
            PyArray_DIMS( *out_gradients )[0] != PyArray_DIMS( in_activations )[0] ||
            PyArray_DIMS( *out_gradients )[1] != PyArray_DIMS( in_activations )[1] ||
            PyArray_DIMS( *out_gradients )[2] != PyArray_DIMS( in_activations )[2] )
        {
            // Existing matrix is the wrong size. Make a new one.
            // Decrement ref counter to existing array
            Py_XDECREF( *out_gradients );
            // Allocate new array
            *out_gradients = (PyArrayObject *) PyArray_ZEROS(3, PyArray_DIMS( in_activations ),
                NPY_FLOAT32, 0);

            if ( NULL == (*out_gradients) )
            {
                // Destroy previous CTC context before returning exception
                ctc_context_destroy( context );

                PyErr_Format( PyExc_MemoryError,
                    "ConnectionistTemporalClassification: Could not allocate memory for CTC gradients!" );
                return 1;
            }
        }
        gradients = (npy_float32 *) PyArray_DATA( *out_gradients );
    }
//.........这里部分代码省略.........

static int
fortran_setattr(PyFortranObject *fp, char *name, PyObject *v) {
    int i,j,flag;
    PyArrayObject *arr = NULL;
    for (i=0,j=1;i<fp->len && (j=strcmp(name,fp->defs[i].name));i++);
    if (j==0) {
        if (fp->defs[i].rank==-1) {
            PyErr_SetString(PyExc_AttributeError,"over-writing fortran routine");
            return -1;
        }
        if (fp->defs[i].func!=NULL) { /* is allocatable array */
            npy_intp dims[F2PY_MAX_DIMS];
            int k;
            save_def = &fp->defs[i];
            if (v!=Py_None) {     /* set new value (reallocate if needed --
                                     see f2py generated code for more
                                     details ) */
                for(k=0;k<fp->defs[i].rank;k++) dims[k]=-1;
                if ((arr = array_from_pyobj(fp->defs[i].type,dims,fp->defs[i].rank,F2PY_INTENT_IN,v))==NULL)
                    return -1;
                (*(fp->defs[i].func))(&fp->defs[i].rank,PyArray_DIMS(arr),set_data,&flag);
            } else {             /* deallocate */
                for(k=0;k<fp->defs[i].rank;k++) dims[k]=0;
                (*(fp->defs[i].func))(&fp->defs[i].rank,dims,set_data,&flag);
                for(k=0;k<fp->defs[i].rank;k++) dims[k]=-1;
            }
            memcpy(fp->defs[i].dims.d,dims,fp->defs[i].rank*sizeof(npy_intp));
        } else {                     /* not allocatable array */
            if ((arr = array_from_pyobj(fp->defs[i].type,fp->defs[i].dims.d,fp->defs[i].rank,F2PY_INTENT_IN,v))==NULL)
                return -1;
        }
        if (fp->defs[i].data!=NULL) { /* copy Python object to Fortran array */
            npy_intp s = PyArray_MultiplyList(fp->defs[i].dims.d,PyArray_NDIM(arr));
            if (s==-1)
                s = PyArray_MultiplyList(PyArray_DIMS(arr),PyArray_NDIM(arr));
            if (s<0 ||
                (memcpy(fp->defs[i].data,PyArray_DATA(arr),s*PyArray_ITEMSIZE(arr)))==NULL) {
                if ((PyObject*)arr!=v) {
                    Py_DECREF(arr);
                }
                return -1;
            }
            if ((PyObject*)arr!=v) {
                Py_DECREF(arr);
            }
        } else return (fp->defs[i].func==NULL?-1:0);
        return 0; /* succesful */
    }
    if (fp->dict == NULL) {
        fp->dict = PyDict_New();
        if (fp->dict == NULL)
            return -1;
    }
    if (v == NULL) {
        int rv = PyDict_DelItemString(fp->dict, name);
        if (rv < 0)
            PyErr_SetString(PyExc_AttributeError,"delete non-existing fortran attribute");
        return rv;
    }
    else
        return PyDict_SetItemString(fp->dict, name, v);
}

Py::Object _path_module::affine_transform(const Py::Tuple& args)
{
    args.verify_length(2);

    Py::Object vertices_obj = args[0];
    Py::Object transform_obj = args[1];

    PyArrayObject* vertices = NULL;
    PyArrayObject* transform = NULL;
    PyArrayObject* result = NULL;

    try
    {
        vertices = (PyArrayObject*)PyArray_FromObject
                   (vertices_obj.ptr(), PyArray_DOUBLE, 1, 2);
        if (!vertices ||
	    (PyArray_NDIM(vertices) == 2 && PyArray_DIM(vertices, 1) != 2) ||
	    (PyArray_NDIM(vertices) == 1 && PyArray_DIM(vertices, 0) != 2))
            throw Py::ValueError("Invalid vertices array.");

        transform = (PyArrayObject*) PyArray_FromObject
                    (transform_obj.ptr(), PyArray_DOUBLE, 2, 2);
        if (!transform ||
	    PyArray_DIM(transform, 0) != 3 ||
            PyArray_DIM(transform, 1) != 3)
            throw Py::ValueError("Invalid transform.");

        double a, b, c, d, e, f;
        {
            size_t stride0 = PyArray_STRIDE(transform, 0);
            size_t stride1 = PyArray_STRIDE(transform, 1);
            char* row0 = PyArray_BYTES(transform);
            char* row1 = row0 + stride0;

            a = *(double*)(row0);
            row0 += stride1;
            c = *(double*)(row0);
            row0 += stride1;
            e = *(double*)(row0);

            b = *(double*)(row1);
            row1 += stride1;
            d = *(double*)(row1);
            row1 += stride1;
            f = *(double*)(row1);
        }

        result = (PyArrayObject*)PyArray_SimpleNew
          (PyArray_NDIM(vertices), PyArray_DIMS(vertices), PyArray_DOUBLE);
        if (PyArray_NDIM(vertices) == 2)
        {
            size_t n = PyArray_DIM(vertices, 0);
            char* vertex_in = PyArray_BYTES(vertices);
            double* vertex_out = (double*)PyArray_DATA(result);
            size_t stride0 = PyArray_STRIDE(vertices, 0);
            size_t stride1 = PyArray_STRIDE(vertices, 1);
            double x;
            double y;

            for (size_t i = 0; i < n; ++i)
            {
                x = *(double*)(vertex_in);
                y = *(double*)(vertex_in + stride1);

                *vertex_out++ = a*x + c*y + e;
                *vertex_out++ = b*x + d*y + f;

                vertex_in += stride0;
            }
        }
        else
        {
            char* vertex_in = PyArray_BYTES(vertices);
            double* vertex_out = (double*)PyArray_DATA(result);
            size_t stride0 = PyArray_STRIDE(vertices, 0);
            double x;
            double y;
            x = *(double*)(vertex_in);
            y = *(double*)(vertex_in + stride0);
            *vertex_out++ = a*x + c*y + e;
            *vertex_out++ = b*x + d*y + f;
        }
    }
    catch (...)
    {
        Py_XDECREF(vertices);
        Py_XDECREF(transform);
        Py_XDECREF(result);
        throw;
    }

    Py_XDECREF(vertices);
    Py_XDECREF(transform);

    return Py::Object((PyObject*)result, true);
}

//.........这里部分代码省略.........
        int res = PyArray_DTypeFromObjectHelper(objects[i], maxdims - 1,
                                                out_dtype, string_type);
        if (res < 0) {
            Py_DECREF(seq);
            goto fail;
        }
        else if (res > 0) {
            Py_DECREF(seq);
            return res;
        }
    }

    Py_DECREF(seq);

    return 0;


promote_types:
    /* Set 'out_dtype' if it's NULL */
    if (*out_dtype == NULL) {
        if (!string_type && dtype->type_num == NPY_STRING) {
            Py_DECREF(dtype);
            return RETRY_WITH_STRING;
        }
        if (!string_type && dtype->type_num == NPY_UNICODE) {
            Py_DECREF(dtype);
            return RETRY_WITH_UNICODE;
        }
        *out_dtype = dtype;
        return 0;
    }
    /* Do type promotion with 'out_dtype' */
    else {
        PyArray_Descr *res_dtype = PyArray_PromoteTypes(dtype, *out_dtype);
        Py_DECREF(dtype);
        if (res_dtype == NULL) {
            return -1;
        }
        if (!string_type &&
                res_dtype->type_num == NPY_UNICODE &&
                (*out_dtype)->type_num != NPY_UNICODE) {
            Py_DECREF(res_dtype);
            return RETRY_WITH_UNICODE;
        }
        if (!string_type &&
                res_dtype->type_num == NPY_STRING &&
                (*out_dtype)->type_num != NPY_STRING) {
            Py_DECREF(res_dtype);
            return RETRY_WITH_STRING;
        }
        Py_DECREF(*out_dtype);
        *out_dtype = res_dtype;
        return 0;
    }

fail:
    Py_XDECREF(*out_dtype);
    *out_dtype = NULL;
    return -1;
}

#undef RETRY_WITH_STRING
#undef RETRY_WITH_UNICODE

/* new reference */
NPY_NO_EXPORT PyArray_Descr *
_array_typedescr_fromstr(char *c_str)
{
    PyArray_Descr *descr = NULL;
    PyObject *stringobj = PyString_FromString(c_str);

    if (stringobj == NULL) {
        return NULL;
    }
    if (PyArray_DescrConverter(stringobj, &descr) != NPY_SUCCEED) {
        Py_DECREF(stringobj);
        return NULL;
    }
    Py_DECREF(stringobj);
    return descr;
}


NPY_NO_EXPORT char *
index2ptr(PyArrayObject *mp, npy_intp i)
{
    npy_intp dim0;

    if (PyArray_NDIM(mp) == 0) {
        PyErr_SetString(PyExc_IndexError, "0-d arrays can't be indexed");
        return NULL;
    }
    dim0 = PyArray_DIMS(mp)[0];
    if (check_and_adjust_index(&i, dim0, 0, NULL) < 0)
        return NULL;
    if (i == 0) {
        return PyArray_DATA(mp);
    }
    return PyArray_BYTES(mp)+i*PyArray_STRIDES(mp)[0];
}

//.........这里部分代码省略.........
            goto fail;
        }
    }
    else {
        newin = self;
        Py_INCREF(newin);
    }

    /*
     * At this point, newin is a single-segment, aligned, and correct
     * byte-order array of the correct type
     *
     * if ingood == 0, then it is a copy, otherwise,
     * it is the original input.
     */

    /*
     * If we have already made a copy of the data, then use
     * that as the output array
     */
    if (out == NULL && !ingood) {
        out = newin;
    }

    /*
     * Now, we know newin is a usable array for fastclip,
     * we need to make sure the output array is available
     * and usable
     */
    if (out == NULL) {
        Py_INCREF(indescr);
        out = (PyArrayObject*)PyArray_NewFromDescr(Py_TYPE(self),
                                            indescr, PyArray_NDIM(self),
                                            PyArray_DIMS(self),
                                            NULL, NULL,
                                            PyArray_ISFORTRAN(self),
                                            (PyObject *)self);
        if (out == NULL) {
            goto fail;
        }

        outgood = 1;
    }
    else Py_INCREF(out);
    /* Input is good at this point */
    if (out == newin) {
        outgood = 1;
    }


    /* make sure the shape of the output array is the same */
    if (!PyArray_SAMESHAPE(newin, out)) {
        PyErr_SetString(PyExc_ValueError, "clip: Output array must have the"
                        "same shape as the input.");
        goto fail;
    }

    if (!outgood && PyArray_EQUIVALENTLY_ITERABLE(
                            self, out, PyArray_TRIVIALLY_ITERABLE_OP_READ,
                            PyArray_TRIVIALLY_ITERABLE_OP_NOREAD) &&
                        PyArray_CHKFLAGS(out, NPY_ARRAY_ALIGNED) &&
                        PyArray_ISNOTSWAPPED(out) &&
                        PyArray_EquivTypes(PyArray_DESCR(out), indescr)) {
        outgood = 1;
    }

static PyObject * anglesToGVec(PyObject * self, PyObject * args)
{
  PyArrayObject *angs, *bHat_l, *eHat_l, *rMat_c;
  PyArrayObject *gVec_c;
  double chi;
  npy_intp nvecs, rdims[2];

  int nangs, nbhat, nehat, nrmat;
  int da1, db1, de1, dr1, dr2;

  double *angs_ptr, *bHat_l_ptr, *eHat_l_ptr, *rMat_c_ptr;
  double *gVec_c_ptr;

  /* Parse arguments */
  if ( !PyArg_ParseTuple(args,"OOOdO",
			 &angs,
			 &bHat_l, &eHat_l,
			 &chi, &rMat_c)) return(NULL);
  if ( angs == NULL ) return(NULL);

  /* Verify shape of input arrays */
  nangs = PyArray_NDIM(angs);
  nbhat = PyArray_NDIM(bHat_l);
  nehat = PyArray_NDIM(eHat_l);
  nrmat = PyArray_NDIM(rMat_c);

  assert( nangs==2 && nbhat==1 && nehat==1 && nrmat==2 );

  /* Verify dimensions of input arrays */
  nvecs = PyArray_DIMS(angs)[0]; //rows
  da1   = PyArray_DIMS(angs)[1]; //cols

  db1   = PyArray_DIMS(bHat_l)[0];
  de1   = PyArray_DIMS(eHat_l)[0];
  dr1   = PyArray_DIMS(rMat_c)[0];
  dr2   = PyArray_DIMS(rMat_c)[1];

  assert( da1 == 3 );
  assert( db1 == 3 && de1 == 3);
  assert( dr1 == 3 && dr2 == 3);

  /* Allocate C-style array for return data */
  rdims[0] = nvecs; rdims[1] = 3;
  gVec_c = (PyArrayObject*)PyArray_EMPTY(2,rdims,NPY_DOUBLE,0);

  /* Grab pointers to the various data arrays */
  angs_ptr   = (double*)PyArray_DATA(angs);
  bHat_l_ptr = (double*)PyArray_DATA(bHat_l);
  eHat_l_ptr = (double*)PyArray_DATA(eHat_l);
  rMat_c_ptr = (double*)PyArray_DATA(rMat_c);
  gVec_c_ptr = (double*)PyArray_DATA(gVec_c);

  /* Call the actual function */
  anglesToGvec_cfunc(nvecs, angs_ptr,
		     bHat_l_ptr, eHat_l_ptr,
		     chi, rMat_c_ptr,
		     gVec_c_ptr);

  /* Build and return the nested data structure */
  return((PyObject*)gVec_c);
}

mitk::Image::Pointer mitk::PythonService::CopySimpleItkImageFromPython(const std::string &stdvarName)
{
  double*ds = nullptr;
  // access python module
  PyObject *pyMod = PyImport_AddModule((char*)"__main__");
  // global dictionarry
  PyObject *pyDict = PyModule_GetDict(pyMod);
  mitk::Image::Pointer mitkImage = mitk::Image::New();
  mitk::Vector3D spacing;
  mitk::Point3D origin;
  QString command;
  QString varName = QString::fromStdString( stdvarName );

  command.append( QString("%1_numpy_array = sitk.GetArrayFromImage(%1)\n").arg(varName) );
  command.append( QString("%1_spacing = numpy.asarray(%1.GetSpacing())\n").arg(varName) );
  command.append( QString("%1_origin = numpy.asarray(%1.GetOrigin())\n").arg(varName) );
  command.append( QString("%1_dtype = %1_numpy_array.dtype.name\n").arg(varName) );
  command.append( QString("%1_direction = numpy.asarray(%1.GetDirection())\n").arg(varName) );
  command.append( QString("%1_nrComponents = numpy.asarray(%1.GetNumberOfComponentsPerPixel())\n").arg(varName));
  command.append( QString("%1_dtype = %1_numpy_array.dtype.name\n").arg(varName) );


  MITK_DEBUG("PythonService") << "Issuing python command " << command.toStdString();
  this->Execute(command.toStdString(), IPythonService::MULTI_LINE_COMMAND );

  PyObject* py_dtype = PyDict_GetItemString(pyDict,QString("%1_dtype").arg(varName).toStdString().c_str() );
  std::string dtype = PyString_AsString(py_dtype);
  PyArrayObject* py_data = (PyArrayObject*) PyDict_GetItemString(pyDict,QString("%1_numpy_array").arg(varName).toStdString().c_str() );
  PyArrayObject* py_spacing = (PyArrayObject*) PyDict_GetItemString(pyDict,QString("%1_spacing").arg(varName).toStdString().c_str() );
  PyArrayObject* py_origin = (PyArrayObject*) PyDict_GetItemString(pyDict,QString("%1_origin").arg(varName).to