Navigation

← Previous Changeset
Next Changeset →

Changeset 660

Timestamp:

10/06/08 14:22:43 (2 months ago)

Author:

mdroe

Message:

Lots of fixes, mainly to SIP support

Files:

trunk/pywcs/pywcs/_docutil.py (modified) (1 diff)
trunk/pywcs/pywcs/pywcs.py (modified) (8 diffs)
trunk/pywcs/setup.py (modified) (3 diffs)
trunk/pywcs/src/distortion.c (modified) (5 diffs)
trunk/pywcs/src/pipeline.c (modified) (8 diffs)
trunk/pywcs/src/pipeline.h (modified) (2 diffs)
trunk/pywcs/src/pywcs.c (modified) (10 diffs)
trunk/pywcs/src/sip.c (modified) (14 diffs)
trunk/pywcs/src/sip.h (modified) (1 diff)
trunk/pywcs/src/sip_wrap.c (modified) (14 diffs)
trunk/pywcs/src/wcslib_wrap.c (modified) (35 diffs)
trunk/pywcs/test/distortion_test.py (deleted)
trunk/pywcs/test/sip_reference.fits (added)
trunk/pywcs/test/sip_test.py (added)
trunk/pywcs/test/test.py (modified) (4 diffs)

Legend:

: Unmodified
: Added
: Removed
: Modified
: Copied
: Moved

trunk/pywcs/pywcs/_docutil.py

r643	r660
	1	# Copyright (C) 2008 Association of Universities for Research in
	2	# Astronomy (AURA)
	3	#
	4	# Redistribution and use in source and binary forms, with or without
	5	# modification, are permitted provided that the following conditions
	6	# are met:
	7	#
	8	# 1. Redistributions of source code must retain the above
	9	# copyright notice, this list of conditions and the following
	10	# disclaimer.
	11	#
	12	# 2. Redistributions in binary form must reproduce the above
	13	# copyright notice, this list of conditions and the following
	14	# disclaimer in the documentation and/or other materials
	15	# provided with the distribution.
	16	#
	17	# 3. The name of AURA and its representatives may not be used to
	18	# endorse or promote products derived from this software without
	19	# specific prior written permission.
	20	#
	21	# THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR
	22	# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
	23	# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
	24	# ARE DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT,
	25	# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
	26	# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
	27	# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
	28	# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
	29	# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
	30	# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
	31	# OF THE POSSIBILITY OF SUCH DAMAGE.
	32
	33	"""
	34	pywcs-specific utilities for generating boilerplate in docstrings.
	35	"""
	36
1	37	def _fix(content, indent=0):
2	38	lines = content.split('\n')

trunk/pywcs/pywcs/pywcs.py

r652	r660
201	201
202	202	m = int(header["A_ORDER"])
203		a = numpy.zeros((m+1, m+1))
	203	a = numpy.zeros((m+1, m+1), numpy.double)
204	204	for i in range(m+1):
205	205	for j in range(m-i+1):
…	…
207	207
208	208	m = int(header["B_ORDER"])
209		b = numpy.zeros((m+1, m+1))
	209	b = numpy.zeros((m+1, m+1), numpy.double)
210	210	for i in range(m+1):
211	211	for j in range(m-i+1):
…	…
226	226
227	227	m = int(header["AP_ORDER"])
228		ap = numpy.zeros((m+1, m+1))
	228	ap = numpy.zeros((m+1, m+1), numpy.double)
229	229	for i in range(m+1):
230	230	for j in range(m-i+1):
…	…
232	232
233	233	m = int(header["BP_ORDER"])
234		bp = numpy.zeros((m+1, m+1))
	234	bp = numpy.zeros((m+1, m+1), numpy.double)
235	235	for i in range(m+1):
236	236	for j in range(m-i+1):
…	…
247	247	return None
248	248
249		return Sip(a, b, ap, bp, self.wcs.crpix)
	249	if not header.has_key("CRPIX1") or not header.has_key("CRPIX2"):
	250	raise ValueError(
	251	"Header has SIP keywords without CRPIX keywords")
	252
	253	crpix1 = header.get("CRPIX1")
	254	crpix2 = header.get("CRPIX2")
	255
	256	return Sip(a, b, ap, bp, (crpix1, crpix2))
250	257
251	258	def _array_converter(self, func, *args):
…	…
341	348
342	349	def wcs_pix2sky_fits(self, *args):
343		return self._array_converter(lambda x: self.wcs.p2s(x)['world'], *args)
	350	return self._array_converter(lambda x: self.wcs.p2s_fits(x)['world'], *args)
344	351	wcs_pix2sky_fits.__doc__ = """
345	352	wcs_pix2sky_fits(*args) -> sky
…	…
351	358	if self.wcs is None:
352	359	raise ValueError("No basic WCS settings were created.")
353		return self._array_converter(lambda x: self.wcs.s2p(x)['pixel'], *args)
	360	return self._array_converter(lambda x: self.wcs.s2p(x)['pixcrd'], *args)
354	361	wcs_sky2pix.__doc__ = """
355	362	wcs_sky2pix(*args) -> pixel
…	…
379	386	if self.wcs is None:
380	387	raise ValueError("No basic WCS settings were created.")
381		return self._array_converter(lambda x: self.wcs.s2p_fits(x)['pixel'],
	388	return self._array_converter(lambda x: self.wcs.s2p_fits(x)['pixcrd'],
382	389	*args)
383	390	wcs_sky2pix_fits.__doc__ = """

trunk/pywcs/setup.py

r654	r660
123	123
124	124	######################################################################
125		# DISTUTILS SETUP
126		if USE_ASSERTS:
127		define_macros = [('DEBUG', None)]
128		undef_macros = [('NDEBUG')]
129		else:
130		define_macros = [('NDEBUG', None)]
131		undef_macros = [('DEBUG')]
132
	125	# PYWCS-SPECIFIC AND WRAPPER SOURCE FILES
133	126	PYWCS_SOURCES = [ # List of pywcs files to compile
134	127	'distortion.c',
…	…
142	135	'wcslib_wrap.c']
143	136	PYWCS_SOURCES = [join('src', x) for x in PYWCS_SOURCES]
	137
	138	######################################################################
	139	# DISTUTILS SETUP
	140	if USE_ASSERTS:
	141	define_macros = [('DEBUG', None)]
	142	undef_macros = [('NDEBUG')]
	143	extra_compile_args = ["-fno-inline"]
	144	else:
	145	# Define ECHO as nothing to prevent spurious newlines from
	146	# printing within the libwcs parser
	147	define_macros = [('NDEBUG', None), ('ECHO', None)]
	148	undef_macros = [('DEBUG')]
	149	extra_compile_args = []
144	150
145	151	setup(name="pywcs",
…	…
159	165	],
160	166	define_macros=define_macros,
161		undef_macros=undef_macros
	167	undef_macros=undef_macros,
	168	extra_compile_args=extra_compile_args
162	169	)
163	170	]

trunk/pywcs/src/distortion.c

r643	r660
71	71	static inline float
72	72	get_dist(
73		const distortion_lookup_t *lookup,
	73	const distortion_lookup_t *lookup,
74	74	const unsigned int coord[NAXES]) {
75	75	unsigned int cropped[NAXES];
…	…
89	89	static inline double
90	90	image_coord_to_distortion_coord(
91		const distortion_lookup_t *lookup,
	91	const distortion_lookup_t *lookup,
92	92	const unsigned int axis,
93	93	const double img) {
…	…
101	101	lookup->crpix[axis]);
102	102
103		if (result < 0.0)
	103	if (result < 0.0) {
104	104	result = 0.0;
105		~~else if (result >= lookup->naxis[axis])~~
	105	} else if (result >= lookup->naxis[axis]) {
106	106	result = lookup->naxis[axis] - 1.0;
	107	}
107	108
108	109	return result;
…	…
115	116	static inline void
116	117	image_coords_to_distortion_coords(
117		const distortion_lookup_t *lookup,
	118	const distortion_lookup_t *lookup,
118	119	const double img[NAXES],
119	120	/* Output parameters */
…	…
172	173	coord[1] = dist_ifloor[1] + k;
173	174	result += (double)get_dist(lookup, coord) * \
174		calculate_weight(dist_weight[0], l, dist_weight[1], k);
	175	calculate_weight(dist_weight[0], l, dist_weight[1], k);
175	176	}
176	177	}

trunk/pywcs/src/pipeline.c

r644	r660
46	46	pipeline->cpdis[1] = NULL;
47	47	pipeline->wcs = NULL;
	48
	49	/* Temporary buffers */
	50	pipeline->alloc_ncoord = 0;
	51	pipeline->alloc_nelem = 0;
	52	pipeline->tmp = NULL;
	53	pipeline->imgcrd = NULL;
	54	pipeline->phi = NULL;
	55	pipeline->theta = NULL;
	56	pipeline->stat = NULL;
48	57	}
49	58
…	…
60	69	}
61	70
	71	static void
	72	pipeline_free_tmp(pipeline_t* pipeline) {
	73	/* Free all temporary buffers and reset pointers to NULL */
	74	free(pipeline->tmp);
	75	pipeline->tmp = NULL;
	76	free(pipeline->imgcrd);
	77	pipeline->imgcrd = NULL;
	78	free(pipeline->phi);
	79	pipeline->phi = NULL;
	80	free(pipeline->theta);
	81	pipeline->theta = NULL;
	82	free(pipeline->stat);
	83	pipeline->stat = NULL;
	84	}
	85
	86	void
	87	pipeline_free(pipeline_t* pipeline) {
	88	pipeline_free_tmp(pipeline);
	89	}
	90
	91	static int
	92	pipeline_realloc(
	93	pipeline_t* pipeline,
	94	unsigned int ncoord,
	95	unsigned int nelem) {
	96	if (pipeline->alloc_ncoord < ncoord \|\| pipeline->alloc_nelem < nelem) {
	97	pipeline_free_tmp(pipeline);
	98
	99	pipeline->alloc_ncoord = ncoord;
	100	pipeline->alloc_nelem = nelem;
	101
	102	pipeline->imgcrd = malloc(ncoord * nelem * sizeof(double));
	103	if (pipeline->imgcrd == NULL) {
	104	goto out_of_memory;
	105	}
	106
	107	pipeline->phi = malloc(ncoord * sizeof(double));
	108	if (pipeline->phi == NULL) {
	109	goto out_of_memory;
	110	}
	111
	112	pipeline->theta = malloc(ncoord * sizeof(double));
	113	if (pipeline->theta == NULL) {
	114	goto out_of_memory;
	115	}
	116
	117	pipeline->stat = malloc(ncoord * nelem * sizeof(int));
	118	if (pipeline->stat == NULL) {
	119	goto out_of_memory;
	120	}
	121
	122	pipeline->tmp = malloc(ncoord * nelem * sizeof(double));
	123	if (pipeline->tmp == NULL) {
	124	goto out_of_memory;
	125	}
	126	}
	127
	128	return 0;
	129
	130	out_of_memory:
	131	pipeline->alloc_nelem = 0;
	132	pipeline->alloc_ncoord = 0;
	133	pipeline_free_tmp(pipeline);
	134	return 2;
	135	}
	136
62	137	int
63	138	pipeline_all_pixel2world(
…	…
67	142	const double* pixcrd /* [ncoord][nelem] */,
68	143	double* world /* [ncoord][nelem] */) {
69		~~double* tmp = NULL;~~
70		~~double* imgcrd = NULL;~~
71		~~double* phi = NULL;~~
72		~~double* theta = NULL;~~
73		~~int* stat = NULL;~~
74	144	const double* wcs_input = NULL;
75	145	double* wcs_output = NULL;
…	…
90	160
91	161	if (has_wcs) {
92		imgcrd = malloc(ncoord * nelem * sizeof(double));
93		if (imgcrd == NULL) {
94		status = 2;
95		goto pipeline_pixel2world_exit;
96		}
97
98		phi = malloc(ncoord * sizeof(double));
99		if (phi == NULL) {
100		status = 2;
101		goto pipeline_pixel2world_exit;
102		}
103
104		theta = malloc(ncoord * sizeof(double));
105		if (theta == NULL) {
106		status = 2;
107		goto pipeline_pixel2world_exit;
108		}
109
110		stat = malloc(ncoord * nelem * sizeof(int));
111		if (stat == NULL) {
112		status = 2;
113		goto pipeline_pixel2world_exit;
	162	status = pipeline_realloc((pipeline_t*)pipeline, ncoord, nelem);
	163	if (status != 0) {
	164	goto exit;
114	165	}
115	166
116	167	if (has_sip \|\| has_p4) {
117		tmp = malloc(ncoord * nelem * sizeof(double));
118		if (tmp == NULL) {
119		status = 2;
120		goto pipeline_pixel2world_exit;
	168	status = pipeline_pix2foc(pipeline, ncoord, nelem, pixcrd, pipeline->tmp);
	169	if (status != 0) {
	170	goto exit;
121	171	}
122	172
123		status = pipeline_pix2foc(pipeline, ncoord, nelem, pixcrd, tmp);
124		if (status != 0) {
125		goto pipeline_pixel2world_exit;
126		}
127
128		wcs_input = tmp;
	173	wcs_input = pipeline->tmp;
129	174	wcs_output = world;
130	175	} else {
…	…
134	179
135	180	status = wcsp2s(pipeline->wcs, ncoord, nelem,
136		wcs_input, imgcrd, phi, theta, wcs_output, stat);
	181	wcs_input, pipeline->imgcrd, pipeline->phi,
	182	pipeline->theta, wcs_output, pipeline->stat);
137	183	} else {
138	184	if (has_sip \|\| has_p4) {
…	…
141	187	}
142	188
143		pipeline_pixel2world_exit:
144		free(tmp);
145		free(imgcrd);
146		free(phi);
147		free(theta);
148		free(stat);
149
	189	exit:
	190	/* We don't have any cleanup at the moment */
150	191	return status;
151	192	}
…	…
197	238	status = sip_pix2foc(pipeline->sip, 2, ncoord, sip_input, sip_output);
198	239	if (status) {
199		goto ~~pipeline_pix2foc_~~exit;
	240	goto exit;
200	241	}
201	242	}
…	…
204	245	status = p4_pix2foc(2, (void*)pipeline->cpdis, ncoord, p4_input, p4_output);
205	246	if (status) {
206		goto ~~pipeline_pix2foc_~~exit;
207		}
208		}
209
210		~~pipeline_pix2foc_~~exit:
	247	goto exit;
	248	}
	249	}
	250
	251	exit:
211	252	/* We don't have any cleanup at the moment */
212
213	253	return status;
214	254	}

trunk/pywcs/src/pipeline.h

r644	r660
46	46	distortion_lookup_t* cpdis[2];
47	47	struct wcsprm* wcs;
	48
	49	/* Temporary buffers for performing calculations */
	50	unsigned int alloc_ncoord;
	51	unsigned int alloc_nelem;
	52	double* tmp;
	53	double* imgcrd;
	54	double* phi;
	55	double* theta;
	56	int* stat;
48	57	} pipeline_t;
49	58
…	…
64	73	distortion_lookup_t* cpdis[2],
65	74	struct wcsprm* wcs);
	75
	76	/**
	77	Free all the temporary buffers of a pipeline_t. It does not free
	78	the underlying sip_t, distortion_lookup_t or wcsprm objects.
	79	*/
	80	void
	81	pipeline_free(
	82	pipeline_t* pipeline);
66	83
67	84	/**

trunk/pywcs/src/pywcs.c

r653	r660
66	66	*/
67	67
	68	static int
	69	PyWcs_traverse(PyWcs* self, visitproc visit, void* arg) {
	70	Py_VISIT(self->py_sip);
	71	Py_VISIT(self->py_distortion_lookup[0]);
	72	Py_VISIT(self->py_distortion_lookup[1]);
	73	Py_VISIT(self->py_wcsprm);
	74
	75	return 0;
	76	}
	77
68	78	static void
69	79	PyWcs_dealloc(PyWcs* self) {
…	…
72	82	Py_XDECREF(self->py_distortion_lookup[1]);
73	83	Py_XDECREF(self->py_wcsprm);
	84	pipeline_free(&self->x);
74	85	self->ob_type->tp_free((PyObject*)self);
75	86	}
…	…
169	180	world = (PyArrayObject*)PyArray_SimpleNew(2, PyArray_DIMS(pixcrd), PyArray_DOUBLE);
170	181	if (world == NULL) {
171		goto _exit;
	182	goto exit;
172	183	}
173	184
174	185	/* Adjust pixel coordinates to be 1-based */
175		if (do_shift)
	186	if (do_shift) {
176	187	offset_array(pixcrd, 1.0);
	188	}
177	189
178	190	/* Make the call */
…	…
185	197	wcsprm_c2python(self->x.wcs);
186	198	/* Adjust pixel coordinates back to 0-based */
187		if (do_shift)
	199	if (do_shift) {
188	200	offset_array(pixcrd, -1.0);
189
190		_exit:
	201	}
	202
	203	exit:
191	204	Py_XDECREF(pixcrd);
192	205
193	206	if (status == 0 \|\| status == 8) {
194	207	return (PyObject*)world;
195		~~} else if (status > 0 && status < WCS_ERRMSG_MAX) {~~
196		~~PyErr_SetString(*wcs_errexc[status], wcsp2s_errmsg[status]);~~
197		~~return NULL;~~
198		~~} else if (status == -1) {~~
199		~~return NULL;~~
200	208	} else {
201		PyErr_SetString(PyExc_RuntimeError,
202		"Unknown error occurred. Something is seriously wrong.");
203		return NULL;
	209	Py_DECREF(world);
	210	if (status > 0 && status < WCS_ERRMSG_MAX) {
	211	PyErr_SetString(*wcs_errexc[status], wcsp2s_errmsg[status]);
	212	return NULL;
	213	} else if (status == -1) {
	214	return NULL;
	215	} else {
	216	PyErr_SetString(PyExc_RuntimeError,
	217	"Unknown error occurred. Something is seriously wrong.");
	218	return NULL;
	219	}
204	220	}
205	221	}
…	…
233	249	if (PyArray_DIM(pixcrd, 1) != NAXES) {
234	250	PyErr_SetString(PyExc_ValueError, "Pixel array must be an Nx2 array");
235		goto _exit;
	251	goto exit;
236	252	}
237	253
238	254	foccrd = (PyArrayObject*)PyArray_SimpleNew(2, PyArray_DIMS(pixcrd), PyArray_DOUBLE);
239	255	if (foccrd == NULL) {
240		goto _exit;
	256	status = 2;
	257	goto exit;
241	258	}
242	259
…	…
248	265	PyArray_DATA(pixcrd), PyArray_DATA(foccrd));
249	266
250		~~if (status) {~~
251		~~Py_XDECREF(foccrd);~~
252		~~goto _exit;~~
253		}
254
255	267	if (do_shift) {
256	268	offset_array(pixcrd, -1.0);
257	269	}
258	270
259		_exit:
	271	exit:
260	272
261	273	Py_XDECREF(pixcrd);
262	274
263		if (status == 0 ~~\|\| status == 8~~) {
	275	if (status == 0) {
264	276	return (PyObject*)foccrd;
265		~~} else if (status > 0 && status < WCS_ERRMSG_MAX) {~~
266		~~PyErr_SetString(*wcs_errexc[status], wcsp2s_errmsg[status]);~~
267		~~return NULL;~~
268		~~} else if (status == -1) {~~
269		~~return NULL;~~
270	277	} else {
271		PyErr_SetString(PyExc_RuntimeError,
272		"Unknown error occurred. Something is seriously wrong.");
273		return NULL;
	278	Py_XDECREF(foccrd);
	279	if (status > 0 && status < WCS_ERRMSG_MAX) {
	280	PyErr_SetString(*wcs_errexc[status], wcsp2s_errmsg[status]);
	281	return NULL;
	282	} else if (status == -1) {
	283	return NULL;
	284	} else {
	285	PyErr_SetString(PyExc_RuntimeError,
	286	"Unknown error occurred. Something is seriously wrong.");
	287	return NULL;
	288	}
274	289	}
275	290	}
…	…
282	297	static PyObject*
283	298	PyWcs_p4_pix2foc_fits(PyWcs* self, PyObject* arg, PyObject* kwds) {
284		return PyWcs_p4_pix2foc_generic(self, arg, 1);
	299	return PyWcs_p4_pix2foc_generic(self, arg, 0);
285	300	}
286	301
…	…
316	331	(double*)PyArray_DATA(foccrd));
317	332
318		~~if (status) {~~
319		~~Py_DECREF(foccrd);~~
320		~~goto _exit;~~
321		}
322
323	333	if (do_shift) {
324	334	offset_array(pixcrd, -1.0);
…	…
329	339	Py_XDECREF(pixcrd);
330	340
331		if (status == 0 ~~\|\| status == 8~~) {
	341	if (status == 0) {
332	342	return (PyObject*)foccrd;
333		~~} else if (status > 0 && status < WCS_ERRMSG_MAX) {~~
334		~~PyErr_SetString(*wcs_errexc[status], wcsp2s_errmsg[status]);~~
335		~~return NULL;~~
336		~~} else if (status == -1) {~~
337		~~return NULL;~~
338	343	} else {
339		PyErr_SetString(PyExc_RuntimeError,
340		"Unknown error occurred. Something is seriously wrong.");
341		return NULL;
	344	Py_XDECREF(foccrd);
	345	if (status > 0 && status < WCS_ERRMSG_MAX) {
	346	PyErr_SetString(*wcs_errexc[status], wcsp2s_errmsg[status]);
	347	return NULL;
	348	} else if (status == -1) {
	349	return NULL;
	350	} else {
	351	PyErr_SetString(PyExc_RuntimeError,
	352	"Unknown error occurred. Something is seriously wrong.");
	353	return NULL;
	354	}
342	355	}
343	356	}
…	…
350	363	static PyObject*
351	364	PyWcs_pix2foc_fits(PyWcs* self, PyObject* arg, PyObject* kwds) {
352		return PyWcs_pix2foc_generic(self, arg, 1);
	365	return PyWcs_pix2foc_generic(self, arg, 0);
353	366	}
354	367

trunk/pywcs/src/sip.c

r643	r660
41	41	#include <string.h>
42	42
43		void sip_clear(sip_t* sip) {
	43	void
	44	sip_clear(sip_t* sip) {
44	45	assert(sip);
45	46
…	…
52	53	sip->bp_order = 0;
53	54	sip->bp = NULL;
	55	sip->crpix[0] = 0.0;
	56	sip->crpix[1] = 0.0;
	57	sip->scratch = NULL;
54	58	}
55	59
…	…
62	66	const unsigned int bp_order, const double* bp,
63	67	const double crpix[2]) {
64		int status = 0;
65		size_t a_size = 0;
66		size_t b_size = 0;
67		size_t ap_size = 0;
68		size_t bp_size = 0;
69		unsigned int scratch_sizes[4];
70		size_t i;
	68	unsigned int a_size = 0;
	69	unsigned int b_size = 0;
	70	unsigned int ap_size = 0;
	71	unsigned int bp_size = 0;
71	72	unsigned int scratch_size = 0;
	73	int status = 0;
72	74
73	75	assert(sip);
74		assert(sip->a == NULL);
75		assert(sip->b == NULL);
76		assert(sip->ap == NULL);
77		assert(sip->bp == NULL);
78
79		if (!((a == NULL) ^ (b == NULL)) \|\|
80		!((ap == NULL) ^ (bp == NULL))) {
	76	sip_clear(sip);
	77
	78	/* We we have one of A/B or AP/BP, we must have both. */
	79	if (((a == NULL) ^ (b == NULL)) \|\|
	80	((ap == NULL) ^ (bp == NULL))) {
81	81	return 6;
82	82	}
…	…
88	88	if (sip->a == NULL) {
89	89	status = 2;
90		goto ~~sip_init_~~exit;
	90	goto exit;
91	91	}
92	92	memcpy(sip->a, a, a_size);
	93	if (a_order > scratch_size) {
	94	scratch_size = a_order;
	95	}
93	96
94	97	sip->b_order = b_order;
…	…
97	100	if (sip->b == NULL) {
98	101	status = 2;
99		goto ~~sip_init_~~exit;
	102	goto exit;
100	103	}
101	104	memcpy(sip->b, b, b_size);
	105	if (b_order > scratch_size) {
	106	scratch_size = b_order;
	107	}
102	108	}
103	109
…	…
108	114	if (sip->ap == NULL) {
109	115	status = 2;
110		goto ~~sip_init_~~exit;
	116	goto exit;
111	117	}
112	118	memcpy(sip->ap, ap, ap_size);
	119	if (ap_order > scratch_size) {
	120	scratch_size = ap_order;
	121	}
113	122
114	123	sip->bp_order = bp_order;
…	…
117	126	if (sip->bp == NULL) {
118	127	status = 2;
119		goto ~~sip_init_~~exit;
	128	goto exit;
120	129	}
121	130	memcpy(sip->bp, bp, bp_size);
	131	if (bp_order > scratch_size) {
	132	scratch_size = bp_order;
	133	}
	134	}
	135
	136	if (scratch_size > 0) {
	137	scratch_size = (scratch_size + 1) * sizeof(double);
	138	sip->scratch = malloc(scratch_size);
	139	if (sip->scratch == NULL) {
	140	status = 2;
	141	goto exit;
	142	}
122	143	}
123	144
…	…
125	146	sip->crpix[1] = crpix[1];
126	147
127		/* Find the maximum order to create for our scratch memory space */
128		scratch_sizes[0] = a_order + 1;
129		scratch_sizes[1] = b_order + 1;
130		scratch_sizes[2] = ap_order + 1;
131		scratch_sizes[3] = bp_order + 1;
132		for (i = 0; i < 4; ++i) {
133		if (scratch_size > scratch_sizes[i]) {
134		scratch_size = scratch_sizes[i];
135		}
136		}
137
138		sip->scratch = malloc(scratch_size * sizeof(double));
139		if (sip->scratch == NULL) {
140		status = 2;
141		goto sip_init_exit;
142		}
143
144		sip_init_exit:
	148	exit:
145	149	if (status != 0) {
146	150	sip_free(sip);
…	…
168	172	const unsigned int order,
169	173	const double* matrix,
170		const unsigned int x,
171		const unsigned int y) {
172		assert(x <= order);
173		assert(y <= order);
174
175		return matrix[y * (order + 1) + x];
	174	const int x,
	175	const int y) {
	176	int index;
	177	assert(x >= 0 && x <= order);
	178	assert(y >= 0 && y <= order);
	179	index = x * (order + 1) + y;
	180	assert(index >= 0 && index < (order + 1) * (order + 1));
	181
	182	return matrix[index];
176	183	}
177	184
…	…
180	187	const unsigned int naxes,
181	188	const unsigned int nelem,
182		const unsigned int ~~a_order~~,
	189	const unsigned int m,
183	190	const double* a,
184		const unsigned int ~~b_order~~,
	191	const unsigned int n,
185	192	const double* b,
186	193	const double crpix[2],
…	…
188	195	const double* input /* [NAXES][nelem] */,
189	196	double* output /* [NAXES][nelem] */) {
190		~~size_~~t i, j, k;
	197	int i, j, k;
191	198	double x, y, tmp_x, tmp_y;
192	199	double sum;
…	…
199	206	assert(output);
200	207
201		if (a == NULL \|\| b == NULL) {
	208	/* Avoid segfaults */
	209	if (input == NULL \|\| output == NULL \|\| tmp == NULL \|\| crpix == NULL) {
	210	return 1;
	211	}
	212
	213	/* If we have one, we must have both... */
	214	if ((a == NULL) ^ (b == NULL)) {
202	215	return 6;
203	216	}
204	217
205		if (input == NULL \|\| output == NULL) {
206		return 1;
	218	/* If no distortion, just copy values */
	219	if (a == NULL /* && b == NULL ... implied */) {
	220	for (i = 0; i < nelem; ++i) {
	221	output[i << 1] = input[i << 1];
	222	output[(i << 1) + 1] = input[(i << 1) + 1];
	223	}
	224	return 0;
207	225	}
208	226
209	227	for (i = 0; i < nelem; ++i) {
210		x = input[i << 1] ~~+ 1.0~~;
211		y = input[(i << 1) + 1] ~~+ 1.0~~;
212
213		tmp_x = ~~crpix[0] - x~~;
214		tmp_y = ~~crpix[1] - y~~;
215
216		for (j = 0; j <= ~~a_order~~; ++j) {
217		tmp[j] = lu(~~a_order, a, a_order -~~ j, j);
218		for (k = j - 1; k >= 0; --k) {
219		tmp[j] = ~~tmp_y * tmp[j] + lu(a_order, a, a_order -~~ j, k);
	228	x = input[i << 1];
	229	y = input[(i << 1) + 1];
	230
	231	tmp_x = x - crpix[0];
	232	tmp_y = y - crpix[1];
	233
	234	for (j = 0; j <= m; ++j) {
	235	tmp[j] = lu(m, a, m-j, j);
	236	for (k = j-1; k >= 0; --k) {
	237	tmp[j] = (tmp_y * tmp[j]) + lu(m, a, m-j, k);
220	238	}
221	239	}
222	240
	241	/* Don't know why this loop is convoluted */
223	242	sum = tmp[0];
224		for (~~i = a_order; i > 0; --i~~) {
225		sum = tmp_x * sum + tmp[~~a_order - i~~ + 1];
	243	for (j = m; j > 0; --j) {
	244	sum = tmp_x * sum + tmp[m - j + 1];
226	245	}
227	246	output[i << 1] = sum + x;
228	247
229		for (j = 0; j <= ~~b_order~~; ++j) {
230		tmp[j] = lu(~~b_order, b, b_order -~~ j, j);
231		for (k = j - 1; k >= 0; --k) {
232		tmp[j] = tmp_y * tmp[j] + lu(~~b_order, b, b_order -~~ j, k);
	248	for (j = 0; j <= n; ++j) {
	249	tmp[j] = lu(n, b, n-j, j);
	250	for (k = j-1; k >= 0; --k) {
	251	tmp[j] = tmp_y * tmp[j] + lu(n, b, n-j, k);
233	252	}
234	253	}
235	254
	255	/* Don't know why this loop is convoluted */
236	256	sum = tmp[0];
237		for (~~i = b_order; i > 0; --i~~) {
238		sum = tmp_x * sum + tmp[~~b_order - i~~ + 1];
	257	for (j = n; j > 0; --j) {
	258	sum = tmp_x * sum + tmp[n - j + 1];
239	259	}
240	260	output[(i << 1) + 1] = sum + y;
…	…
259	279	sip->b_order, sip->b,
260	280	sip->crpix,
261		sip->scratch,
	281	(double *)sip->scratch,
262	282	pix, foc);
263	283	}
…	…
278	298	sip->bp_order, sip->bp,
279	299	sip->crpix,
280		sip->scratch,
	300	(double *)sip->scratch,
281	301	foc, pix);
282	302	}

trunk/pywcs/src/sip.h

r643	r660
40	40	typedef struct {
41	41	unsigned int a_order;
42		double* a; ~~/* [a_order+1][a_order+1] */~~
	42	double* a;
43	43	unsigned int b_order;
44		double* b; ~~/* [b_order+1][b_order+1] */~~
	44	double* b;
45	45	unsigned int ap_order;
46		double* ap; ~~/* [ap_order+1][ap_order+1] */~~
	46	double* ap;
47	47	unsigned int bp_order;
48		double* bp; ~~/* [bp_order+1][bp_order+1] */~~
	48	double* bp;
49	49	double crpix[2];
50	50	double* scratch;

trunk/pywcs/src/sip_wrap.c

r652	r660
37	37	#include "sip_wrap.h"
38	38	#include "docstrings.h"
	39	#include "wcs.h"
39	40
40	41	static void
…	…
56	57
57	58	static int
58		convert_matrix(PyObject* pyobj, PyArrayObject** array) {
59		array = (PyArrayObject)PyArray_ContiguousFromAny(pyobj, PyArray_DOUBLE, 2, 2);
	59	convert_matrix(
	60	PyObject* pyobj,
	61	PyArrayObject** array,
	62	double** data,
	63	size_t* order) {
	64	if (pyobj == Py_None) {
	65	*array = NULL;
	66	*data = NULL;
	67	*order = 0;
	68	return 0;
	69	}
	70
	71	array = (PyArrayObject)PyArray_ContiguousFromAny(
	72	pyobj, PyArray_DOUBLE, 2, 2);
60	73	if (*array == NULL)
61	74	return -1;
…	…
66	79	return -1;
67	80	}
	81
	82	data = (double)PyArray_DATA(*array);
	83	order = PyArray_DIM(array, 0) - 1;
68	84
69	85	return 0;
…	…
82	98	PyArrayObject* bp = NULL;
83	99	PyArrayObject* crpix = NULL;
84		int result = 0;
	100	double* a_data = NULL;
	101	double* b_data = NULL;
	102	double* ap_data = NULL;
	103	double* bp_data = NULL;
	104	size_t a_order = 0;
	105	size_t b_order = 0;
	106	size_t ap_order = 0;
	107	size_t bp_order = 0;
	108	int status = -1;
85	109
86	110	if (!PyArg_ParseTuple(args, "OOOOO:Sip.__init__",
…	…
89	113	}
90	114
91		if (convert_matrix(py_a, &a) \|\|
92		convert_matrix(py_b, &b) \|\|
93		convert_matrix(py_ap, &ap) \|\|
94		convert_matrix(py_bp, &bp)) {
95		result = -1;
96		goto PySip_init_exit;
	115	if (convert_matrix(py_a, &a, &a_data, &a_order) \|\|
	116	convert_matrix(py_b, &b, &b_data, &b_order) \|\|
	117	convert_matrix(py_ap, &ap, &ap_data, &ap_order) \|\|
	118	convert_matrix(py_bp, &bp, &bp_data, &bp_order)) {
	119	goto exit;
97	120	}
98	121
…	…
100	123	1, 1);
101	124	if (crpix == NULL) {
102		result = -1;
103		goto PySip_init_exit;
	125	goto exit;
104	126	}
105	127
106	128	if (PyArray_DIM(crpix, 0) != 2) {
107		PyErr_SetString(PyExc_ValueError,
108		"CRPIX wrong length");
109		result = -1;
110		goto PySip_init_exit;
111		}
112
113		if (sip_init(&self->x,
114		PyArray_DIM(a, 0) - 1, PyArray_DATA(a),
115		PyArray_DIM(b, 0) - 1, PyArray_DATA(b),
116		PyArray_DIM(ap, 0) - 1, PyArray_DATA(ap),
117		PyArray_DIM(bp, 0) - 1, PyArray_DATA(bp),
118		PyArray_DATA(crpix))) {
119		result = -1;
120		goto PySip_init_exit;
121		}
122
123		PySip_init_exit:
	129	PyErr_SetString(PyExc_ValueError, "CRPIX wrong length");
	130	goto exit;
	131	}
	132
	133	status = sip_init(&self->x,
	134	a_order, a_data,
	135	b_order, b_data,
	136	ap_order, ap_data,
	137	bp_order, bp_data,
	138	PyArray_DATA(crpix));
	139
	140	exit:
124	141	Py_XDECREF(a);
125	142	Py_XDECREF(b);
…	…
128	145	Py_XDECREF(crpix);
129	146
130		return result;
	147	if (status == 0) {
	148	return 0;
	149	} else if (status > 0 && status < WCS_ERRMSG_MAX) {
	150	PyErr_SetString(*wcs_errexc[status], wcsp2s_errmsg[status]);
	151	return -1;
	152	} else if (status == -1) {
	153	return -1;
	154	} else {
	155	PyErr_SetString(PyExc_RuntimeError,
	156	"Unknown error occurred. Something is seriously wrong.");
	157	return -1;
	158	}
131	159	}
132	160
…	…
135	163	PyArrayObject* pixcrd = NULL;
136	164	PyArrayObject* foccrd = NULL;
137		PyObject* result = NULL;
	165	int status = -1;
	166
	167	if (self->x.a == NULL \|\| self->x.b == NULL) {
	168	PyErr_SetString(
	169	PyExc_ValueError,
	170	"SIP object does not have coefficients for pix2foc transformation (A and B)");
	171	return NULL;
	172	}
138	173
139	174	pixcrd = (PyArrayObject*)PyArray_ContiguousFromAny(arg, PyArray_DOUBLE, 2, 2);
140	175	if (pixcrd == NULL) {
141		goto _exit;
	176	goto exit;
142	177	}
143	178
144	179	if (PyArray_DIM(pixcrd, 1) != 2) {
145	180	PyErr_SetString(PyExc_ValueError, "Pixel array must be an Nx2 array");
146		goto _exit;
	181	goto exit;
147	182	}
148	183
…	…
150	185	PyArray_DOUBLE);
151	186	if (foccrd == NULL) {
152		goto _exit;
153		}
154
155		if (do_shift)
	187	status = 2;
	188	goto exit;
	189	}
	190
	191	if (do_shift) {
156	192	offset_array(pixcrd, 1.0);
157
158		if (sip_pix2foc(&self->x,
159		PyArray_DIM(pixcrd, 0),
160		PyArray_DIM(pixcrd, 1),
161		(double*)PyArray_DATA(pixcrd),
162		(double*)PyArray_DATA(foccrd))) {
163		PyErr_SetString(PyExc_RuntimeError, "Error correcting distortion");
	193	}
	194