[Yt-svn] yt: 3 new changesets
hg at spacepope.org
hg at spacepope.org
Wed Feb 9 14:54:32 PST 2011
hg Repository: yt
details: yt/rev/15dbdcd32e69
changeset: 3726:15dbdcd32e69
user: Matthew Turk <matthewturk at gmail.com>
date:
Wed Feb 09 15:48:38 2011 -0500
description:
Adding in mk_pix2xy and making libconfig extension a bit clearer.
hg Repository: yt
details: yt/rev/3b2b8eb96702
changeset: 3727:3b2b8eb96702
user: Matthew Turk <matthewturk at gmail.com>
date:
Wed Feb 09 16:53:40 2011 -0500
description:
Adding in mk_xy2pix, vec2pix_nest, and adding on array interfaces to the
vec2pix and pix2vec routines.
hg Repository: yt
details: yt/rev/489e024f45fc
changeset: 3728:489e024f45fc
user: Matthew Turk <matthewturk at gmail.com>
date:
Wed Feb 09 17:54:21 2011 -0500
description:
Adding initial pass at HEALpixCamera class. This returns an image of
dimensionality (12*nside**2, ...). The mollview function in a modified version
of healpy's projection routines will plot this in matplotlib.
Right now, this only supports un-ordered projections -- for instance, column
density, weighted average, etc. To move beyond this I thought that I could do
corner/radius calculation, but now it's clear that it will require a
calculation of the minimum intersection radius of grids with spheres of
expanding radii, which might require a bit more care. Fortunately bricking
ensures we only have to touch each grid once, so we should only have to sort
one extra brick to get them in their final order.
diffstat:
yt/utilities/_amr_utils/VolumeIntegrator.pyx | 51 +++++++++
yt/utilities/_amr_utils/healpix_interface.pxd | 8 +-
yt/utilities/_amr_utils/healpix_mk_pix2xy.c | 74 +++++++++++++
yt/utilities/_amr_utils/healpix_mk_xy2pix.c | 66 ++++++++++++
yt/utilities/_amr_utils/healpix_vec2pix_nest.c | 137 +++++++++++++++++++++++++
yt/utilities/setup.py | 18 +--
yt/visualization/volume_rendering/camera.py | 65 +++++++++++-
7 files changed, 399 insertions(+), 20 deletions(-)
diffs (truncated from 489 to 300 lines):
diff -r 43ac8722363e -r 489e024f45fc yt/utilities/_amr_utils/VolumeIntegrator.pyx
--- a/yt/utilities/_amr_utils/VolumeIntegrator.pyx Wed Feb 09 15:40:01 2011 -0500
+++ b/yt/utilities/_amr_utils/VolumeIntegrator.pyx Wed Feb 09 17:54:21 2011 -0500
@@ -73,6 +73,57 @@
np.float64_t eval_gradient(int *ds, int *ci, np.float64_t *dp,
np.float64_t *data, np.float64_t *grad)
+def hp_pix2vec_nest(long nside, long ipix):
+ cdef double v[3]
+ healpix_interface.pix2vec_nest(nside, ipix, v)
+ cdef np.ndarray[np.float64_t, ndim=1] tr = np.empty((3,), dtype='float64')
+ tr[0] = v[0]
+ tr[1] = v[1]
+ tr[2] = v[2]
+ return tr
+
+def arr_pix2vec_nest(long nside,
+ np.ndarray[np.int64_t, ndim=1] aipix):
+ cdef int n = aipix.shape[0]
+ cdef int i
+ cdef double v[3]
+ cdef long ipix
+ cdef np.ndarray[np.float64_t, ndim=2] tr = np.zeros((n, 3), dtype='float64')
+ for i in range(n):
+ ipix = aipix[i]
+ healpix_interface.pix2vec_nest(nside, ipix, v)
+ tr[i,0] = v[0]
+ tr[i,1] = v[1]
+ tr[i,2] = v[2]
+ return tr
+
+def hp_vec2pix_nest(long nside, double x, double y, double z):
+ cdef double v[3]
+ v[0] = x
+ v[1] = y
+ v[2] = z
+ cdef long ipix
+ healpix_interface.vec2pix_nest(nside, v, &ipix)
+ return ipix
+
+
+def arr_vec2pix_nest(long nside,
+ np.ndarray[np.float64_t, ndim=1] x,
+ np.ndarray[np.float64_t, ndim=1] y,
+ np.ndarray[np.float64_t, ndim=1] z):
+ cdef int n = x.shape[0]
+ cdef int i
+ cdef double v[3]
+ cdef long ipix
+ cdef np.ndarray[np.int64_t, ndim=1] tr = np.zeros(n, dtype='int64')
+ for i in range(n):
+ v[0] = x[i]
+ v[1] = y[i]
+ v[2] = z[i]
+ healpix_interface.vec2pix_nest(nside, v, &ipix)
+ tr[i] = ipix
+ return tr
+
cdef class star_kdtree_container:
cdef kdtree_utils.kdtree *tree
cdef public np.float64_t sigma
diff -r 43ac8722363e -r 489e024f45fc yt/utilities/_amr_utils/healpix_interface.pxd
--- a/yt/utilities/_amr_utils/healpix_interface.pxd Wed Feb 09 15:40:01 2011 -0500
+++ b/yt/utilities/_amr_utils/healpix_interface.pxd Wed Feb 09 17:54:21 2011 -0500
@@ -30,9 +30,5 @@
from stdio cimport fopen, fclose, FILE
cdef extern from "healpix_vectors.h":
- #int mkPix2xy(long *ipix2x, long *ipix2y)
- #int mk_xy2pix(int *x2pix, int *y2pix)
- int pix2vec_nest(long nside, long ipix, double *v, double (*vertex)[3])
- #int vec2pix_nest( long nside, double *vec, long *ipix)
- #int pix2coord_nest( long nside, long ipix, int xsize, int ysize, int &npts,
- # int * &xpolygon, int * &ypolygon, int &draw_poly)
+ int pix2vec_nest(long nside, long ipix, double *v)
+ void vec2pix_nest(long nside, double *vec, long *ipix)
diff -r 43ac8722363e -r 489e024f45fc yt/utilities/_amr_utils/healpix_mk_pix2xy.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/utilities/_amr_utils/healpix_mk_pix2xy.c Wed Feb 09 17:54:21 2011 -0500
@@ -0,0 +1,74 @@
+/* -----------------------------------------------------------------------------
+ *
+ * Copyright (C) 1997-2010 Krzysztof M. Gorski, Eric Hivon,
+ * Benjamin D. Wandelt, Anthony J. Banday,
+ * Matthias Bartelmann,
+ * Reza Ansari & Kenneth M. Ganga
+ *
+ *
+ * This file is part of HEALPix.
+ *
+ * HEALPix is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * HEALPix is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with HEALPix; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+ *
+ * For more information about HEALPix see http://healpix.jpl.nasa.gov
+ *
+ *----------------------------------------------------------------------------- */
+/* mk_pix2xy.c
+ *
+ */
+
+/* Standard Includes */
+#include <math.h>
+
+void mk_pix2xy(int *pix2x, int *pix2y) {
+
+ /* =======================================================================
+ * subroutine mk_pix2xy
+ * =======================================================================
+ * constructs the array giving x and y in the face from pixel number
+ * for the nested (quad-cube like) ordering of pixels
+ *
+ * the bits corresponding to x and y are interleaved in the pixel number
+ * one breaks up the pixel number by even and odd bits
+ * =======================================================================
+ */
+
+ int i, kpix, jpix, IX, IY, IP, ID;
+ for (i = 0; i < 1023; i++) pix2x[i]=0;
+
+ for( kpix=0;kpix<1024;kpix++ ) {
+ jpix = kpix;
+ IX = 0;
+ IY = 0;
+ IP = 1 ;// ! bit position (in x and y)
+ while( jpix!=0 ){// ! go through all the bits
+ ID = (int)fmod(jpix,2);// ! bit value (in kpix), goes in ix
+ jpix = jpix/2;
+ IX = ID*IP+IX;
+
+ ID = (int)fmod(jpix,2);// ! bit value (in kpix), goes in iy
+ jpix = jpix/2;
+ IY = ID*IP+IY;
+
+ IP = 2*IP;// ! next bit (in x and y)
+ }
+
+ pix2x[kpix] = IX;// ! in 0,31
+ pix2y[kpix] = IY;// ! in 0,31
+ }
+
+ /* Later */
+ return;
+}
diff -r 43ac8722363e -r 489e024f45fc yt/utilities/_amr_utils/healpix_mk_xy2pix.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/utilities/_amr_utils/healpix_mk_xy2pix.c Wed Feb 09 17:54:21 2011 -0500
@@ -0,0 +1,66 @@
+/* -----------------------------------------------------------------------------
+ *
+ * Copyright (C) 1997-2010 Krzysztof M. Gorski, Eric Hivon,
+ * Benjamin D. Wandelt, Anthony J. Banday,
+ * Matthias Bartelmann,
+ * Reza Ansari & Kenneth M. Ganga
+ *
+ *
+ * This file is part of HEALPix.
+ *
+ * HEALPix is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * HEALPix is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with HEALPix; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+ *
+ * For more information about HEALPix see http://healpix.jpl.nasa.gov
+ *
+ *----------------------------------------------------------------------------- */
+/* Standard Includes */
+#include <math.h>
+
+void mk_xy2pix(int *x2pix, int *y2pix) {
+ /* =======================================================================
+ * subroutine mk_xy2pix
+ * =======================================================================
+ * sets the array giving the number of the pixel lying in (x,y)
+ * x and y are in {1,128}
+ * the pixel number is in {0,128**2-1}
+ *
+ * if i-1 = sum_p=0 b_p * 2^p
+ * then ix = sum_p=0 b_p * 4^p
+ * iy = 2*ix
+ * ix + iy in {0, 128**2 -1}
+ * =======================================================================
+ */
+ int i, K,IP,I,J,ID;
+
+ for(i = 0; i < 127; i++) x2pix[i] = 0;
+ for( I=1;I<=128;I++ ) {
+ J = I-1;// !pixel numbers
+ K = 0;//
+ IP = 1;//
+ truc : if( J==0 ) {
+ x2pix[I-1] = K;
+ y2pix[I-1] = 2*K;
+ }
+ else {
+ ID = (int)fmod(J,2);
+ J = J/2;
+ K = IP*ID+K;
+ IP = IP*4;
+ goto truc;
+ }
+ }
+
+}
+
diff -r 43ac8722363e -r 489e024f45fc yt/utilities/_amr_utils/healpix_vec2pix_nest.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/utilities/_amr_utils/healpix_vec2pix_nest.c Wed Feb 09 17:54:21 2011 -0500
@@ -0,0 +1,137 @@
+/* -----------------------------------------------------------------------------
+ *
+ * Copyright (C) 1997-2010 Krzysztof M. Gorski, Eric Hivon,
+ * Benjamin D. Wandelt, Anthony J. Banday,
+ * Matthias Bartelmann,
+ * Reza Ansari & Kenneth M. Ganga
+ *
+ *
+ * This file is part of HEALPix.
+ *
+ * HEALPix is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * HEALPix is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with HEALPix; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+ *
+ * For more information about HEALPix see http://healpix.jpl.nasa.gov
+ *
+ *----------------------------------------------------------------------------- */
+/* vec2pix_nest.c */
+
+/* Standard Includes */
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+/* Local Includes */
+#include "healpix_vectors.h"
+
+void vec2pix_nest( const long nside, double *vec, long *ipix) {
+
+ /* =======================================================================
+ * subroutine vec2pix_nest(nside, vec, ipix)
+ * =======================================================================
+ * gives the pixel number ipix (NESTED) corresponding to vector vec
+ *
+ * the computation is made to the highest resolution available (nside=8192)
+ * and then degraded to that required (by integer division)
+ * this doesn't cost more, and it makes sure that the treatement of round-off
+ * will be consistent for every resolution
+ * =======================================================================
+ */
+
+ double z, za, z0, tt, tp, tmp, phi;
+ int face_num,jp,jm;
+ long ifp, ifm;
+ int ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf, ntt;
+ double piover2 = 0.5*M_PI, pi = M_PI, twopi = 2.0*M_PI;
+ int ns_max = 8192;
+ static int x2pix[128], y2pix[128];
+ static char setup_done = 0;
+
+ if( nside<1 || nside>ns_max ) {
+ fprintf(stderr, "%s (%d): nside out of range: %ld\n", __FILE__, __LINE__, nside);
+ exit(0);
+ }
+ if( !setup_done ) {
+ mk_xy2pix(x2pix,y2pix);
+ setup_done = 1;
+ }
+
+ z = vec[2]/sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]);
+ phi = 0.0;
+ if (vec[0] != 0.0 || vec[1] != 0.0) {
More information about the yt-svn
mailing list