Commits

Sam Skillman committed b911ed3

Restoring an interpolated off-axis-projection. This is now controlled by the keyword interpolated=(True/False), with False as the default. Additionally, make the off_axis_projection call simply a wrapper around the much more capable ProjectionCamera class.

Comments (0)

Files changed (2)

yt/utilities/lib/grid_traversal.pyx

     for i in range(imin(3, vc.n_fields)):
         im.rgba[i] += vc.data[i][di] * dl
 
-cdef class ProjectionSampler(ImageSampler):
-    cdef void setup(self, PartitionedGrid pg):
-        self.sampler = projection_sampler
-
 cdef struct VolumeRenderAccumulator:
     int n_fits
     int n_samples
     np.float64_t *light_dir
     np.float64_t *light_rgba
 
+
+cdef class ProjectionSampler(ImageSampler):
+    cdef void setup(self, PartitionedGrid pg):
+        self.sampler = projection_sampler
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.cdivision(True)
+cdef void interpolated_projection_sampler(
+                 VolumeContainer *vc, 
+                 np.float64_t v_pos[3],
+                 np.float64_t v_dir[3],
+                 np.float64_t enter_t,
+                 np.float64_t exit_t,
+                 int index[3],
+                 void *data) nogil:
+    cdef ImageAccumulator *im = <ImageAccumulator *> data
+    cdef VolumeRenderAccumulator *vri = <VolumeRenderAccumulator *> \
+            im.supp_data
+    # we assume this has vertex-centered data.
+    cdef int offset = index[0] * (vc.dims[1] + 1) * (vc.dims[2] + 1) \
+                    + index[1] * (vc.dims[2] + 1) + index[2]
+    cdef np.float64_t slopes[6], dp[3], ds[3]
+    cdef np.float64_t dt = (exit_t - enter_t) / vri.n_samples
+    cdef np.float64_t dvs[6]
+    for i in range(3):
+        dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
+        dp[i] -= index[i] * vc.dds[i] + vc.left_edge[i]
+        dp[i] *= vc.idds[i]
+        ds[i] = v_dir[i] * vc.idds[i] * dt
+    for i in range(vri.n_samples):
+        for j in range(vc.n_fields):
+            dvs[j] = offset_interpolate(vc.dims, dp,
+                    vc.data[j] + offset)
+        for j in range(imin(3, vc.n_fields)):
+            im.rgba[j] += dvs[j] * dt
+        for j in range(3):
+            dp[j] += ds[j]
+
+cdef class InterpolatedProjectionSampler(ImageSampler):
+    cdef VolumeRenderAccumulator *vra
+    cdef public object tf_obj
+    cdef public object my_field_tables
+    def __cinit__(self, 
+                  np.ndarray vp_pos,
+                  np.ndarray vp_dir,
+                  np.ndarray[np.float64_t, ndim=1] center,
+                  bounds,
+                  np.ndarray[np.float64_t, ndim=3] image,
+                  np.ndarray[np.float64_t, ndim=1] x_vec,
+                  np.ndarray[np.float64_t, ndim=1] y_vec,
+                  np.ndarray[np.float64_t, ndim=1] width,
+                  n_samples = 10):
+        ImageSampler.__init__(self, vp_pos, vp_dir, center, bounds, image,
+                               x_vec, y_vec, width)
+        cdef int i
+        # Now we handle tf_obj
+        self.vra = <VolumeRenderAccumulator *> \
+            malloc(sizeof(VolumeRenderAccumulator))
+        self.vra.n_samples = n_samples
+        self.supp_data = <void *> self.vra
+
+    cdef void setup(self, PartitionedGrid pg):
+        self.sampler = interpolated_projection_sampler
+
 @cython.boundscheck(False)
 @cython.wraparound(False)
 @cython.cdivision(True)

yt/visualization/volume_rendering/camera.py

 
 from yt.utilities.lib import \
     PartitionedGrid, ProjectionSampler, VolumeRenderSampler, \
-    LightSourceRenderSampler, \
+    LightSourceRenderSampler, InterpolatedProjectionSampler, \
     arr_vec2pix_nest, arr_pix2vec_nest, arr_ang2pix_nest, \
     pixelize_healpix, arr_fisheye_vectors
 
     def finalize_image(self, image):
         pass
 
-    def _render(self, double_check, num_threads, image, na, sampler):
+    def _render(self, double_check, num_threads, image, sampler):
         pbar = get_pbar("Ray casting", (self.volume.brick_dimensions + 1).prod(axis=-1).sum())
         total_cells = 0
         if double_check:
 
         self.initialize_source()
 
-        image = self._render(double_check, num_threads, image, na, sampler)
+        image = self._render(double_check, num_threads, image, sampler)
 
         self.save_image(fn, clip_ratio, image)
 
         return args
  
 
-    def _render(self, double_check, num_threads, image, na, sampler):
+    def _render(self, double_check, num_threads, image, sampler):
         pbar = get_pbar("Ray casting", (self.volume.brick_dimensions + 1).prod(axis=-1).sum())
         total_cells = 0
         if double_check:
 
         self.volume.initialize_source()
 
-        image = self._render(double_check, num_threads, image, na, sampler)
+        image = self._render(double_check, num_threads, image, sampler)
 
         self.save_image(fn, clim, image)
 
     def finalize_image(self, image):
         image.shape = self.resolution, self.resolution, 3
 
-    def _render(self, double_check, num_threads, image, na, sampler):
+    def _render(self, double_check, num_threads, image, sampler):
         pbar = get_pbar("Ray casting", (self.volume.brick_dimensions + 1).prod(axis=-1).sum())
         total_cells = 0
         if double_check:
                 self.center += dx
             yield self.snapshot()
 
-def off_axis_projection(pf, center, normal_vector, width, resolution,
-                        field, weight = None, num_threads = 0):
-    r"""Project through a parameter file, off-axis, and return the image plane.
-
-    This function will accept the necessary items to integrate through a volume
-    at an arbitrary angle and return the integrated field of view to the user.
-    Note that if a weight is supplied, it will multiply the pre-interpolated
-    values together, then create cell-centered values, then interpolate within
-    the cell to conduct the integration.
-
-    Parameters
-    ----------
-    pf : `~yt.data_objects.api.StaticOutput`
-        This is the parameter file to volume render.
-    center : array_like
-        The current 'center' of the view port -- the focal point for the
-        camera.
-    normal_vector : array_like
-        The vector between the camera position and the center.
-    width : float or list of floats
-        The current width of the image.  If a single float, the volume is
-        cubical, but if not, it is left/right, top/bottom, front/back
-    resolution : int or list of ints
-        The number of pixels in each direction.
-    field : string
-        The field to project through the volume
-    weight : optional, default None
-        If supplied, the field will be pre-multiplied by this, then divided by
-        the integrated value of this field.  This returns an average rather
-        than a sum.
-
-    Returns
-    -------
-    image : array
-        An (N,N) array of the final integrated values, in float64 form.
-
-    Examples
-    --------
-
-    >>> image = off_axis_projection(pf, [0.5, 0.5, 0.5], [0.2,0.3,0.4],
-                      0.2, N, "Temperature", "Density")
-    >>> write_image(na.log10(image), "offaxis.png")
-
-    """
-    # We manually modify the ProjectionTransferFunction to get it to work the
-    # way we want, with a second field that's also passed through.
-    fields = [field]
-
-    if weight is not None:
-        # This is a temporary field, which we will remove at the end.
-        def _wf(f1, w1):
-            def WeightField(field, data):
-                return data[f1].astype("float64") * \
-                       data[w1].astype("float64")
-            return WeightField
-        pf.field_info.add_field("temp_weightfield",
-                    function=_wf(field, weight))
-        fields = ["temp_weightfield", weight]
-    image = na.zeros((resolution, resolution, 3), dtype='float64',
-                      order='C')
-    normal_vector, north_vector, east_vector = ortho_find(normal_vector)
-    unit_vectors = [north_vector, east_vector, normal_vector]
-    back_center= center - 0.5*width * normal_vector
-    rotp = na.concatenate([na.linalg.pinv(unit_vectors).ravel('F'),
-                           back_center])
-    sampler = ProjectionSampler(
-        rotp, normal_vector * width, back_center,
-        (-width/2, width/2, -width/2, width/2),
-        image, north_vector, east_vector,
-        na.array([width, width, width], dtype='float64'))
-    # Calculate the eight corners of the box
-    # Back corners ...
-    mi = pf.domain_right_edge.copy()
-    ma = pf.domain_left_edge.copy()
-    for off1 in [-1, 1]:
-        for off2 in [-1, 1]:
-            for off3 in [-1, 1]:
-                this_point = (center + width/2.0 * off1 * north_vector
-                                     + width/2.0 * off2 * east_vector
-                                     + width/2.0 * off3 * normal_vector)
-                na.minimum(mi, this_point, mi)
-                na.maximum(ma, this_point, ma)
-    # Now we have a bounding box.
-    grids = pf.h.region(center, mi, ma)._grids
-    pb = get_pbar("Sampling ", len(grids))
-    for i,grid in enumerate(grids):
-        data = [(grid[field] * grid.child_mask).astype("float64")
-                for field in fields]
-        pg = PartitionedGrid(
-            grid.id, data,
-            grid.LeftEdge, grid.RightEdge, grid.ActiveDimensions.astype("int64"))
-        grid.clear_data()
-        sampler(pg)
-        pb.update(i)
-    pb.finish()
-    image = sampler.aimage
-    if weight is None:
-        dl = width * pf.units[pf.field_info[field].projection_conversion]
-        image *= dl
-    else:
-        image[:,:,0] /= image[:,:,1]
-        pf.field_info.pop("temp_weightfield")
-    return image[:,:,0]
-
 def allsky_projection(pf, center, radius, nside, field, weight = None,
                       inner_radius = 10, rotation = None):
     r"""Project through a parameter file, through an allsky-method
 
 class ProjectionCamera(Camera):
     def __init__(self, center, normal_vector, width, resolution,
-            field, weight=None, volume=None, le=None, re=None,
-            north_vector=None, pf=None):
-        Camera.__init__(self, center, normal_vector, width, resolution, None,
-                fields = field, pf=pf, volume=1,
-                le=le, re=re, north_vector=north_vector)
+            field, weight=None, volume=None, no_ghost = False, 
+            le=None, re=None,
+            north_vector=None, pf=None, interpolated=False):
+
+        if not interpolated:
+            volume = 1
+
+        self.interpolated = interpolated
         self.field = field
         self.weight = weight
         self.resolution = resolution
 
+        fields = [field]
+        if self.weight is not None:
+            # This is a temporary field, which we will remove at the end.
+            def _make_wf(f, w):
+                def temp_weightfield(a, b):
+                    tr = b[f].astype("float64") * b[w]
+                    return tr
+                return temp_weightfield
+            pf.field_info.add_field("temp_weightfield",
+                function=_make_wf(self.field, self.weight))
+            fields = ["temp_weightfield", self.weight]
+        
+        self.fields = fields
+        self.log_fields = [False]*len(self.fields)
+        Camera.__init__(self, center, normal_vector, width, resolution, None,
+                fields = fields, pf=pf, volume=volume,
+                log_fields=self.log_fields, 
+                le=le, re=re, north_vector=north_vector,
+                no_ghost=no_ghost)
+
     def get_sampler(self, args):
-        sampler = ProjectionSampler(*args)
+        if self.interpolated:
+            sampler = InterpolatedProjectionSampler(*args)
+        else:
+            sampler = ProjectionSampler(*args)
         return sampler
 
     def initialize_source(self):
-        pass
-
+        if self.interpolated:
+            Camera.initialize_source(self)
+        else:
+            pass
 
     def get_sampler_args(self, image):
         width = self.width[2]
         args = (rotp, normal_vector * width, back_center,
             (-width/2, width/2, -width/2, width/2),
             image, north_vector, east_vector,
-            na.array([width, width, width], dtype='float64'))
+            na.array([width, width, width], dtype='float64'),
+            self.sub_samples)
         return args
 
     def finalize_image(self,image):
             image *= dl
         else:
             image[:,:,0] /= image[:,:,1]
-            pf.field_info.pop("temp_weightfield")
         return image[:,:,0]
 
 
-    def _render(self, double_check, num_threads, image, na, sampler):
+    def _render(self, double_check, num_threads, image, sampler):
         # Calculate the eight corners of the box
         # Back corners ...
+        if self.interpolated:
+            return Camera._render(self, double_check, num_threads, image,
+                    sampler)
         pf = self.pf
         width = self.width[2]
         north_vector = self.orienter.unit_vectors[0]
             sampler(pg, num_threads = num_threads)
             pb.update(i)
         pb.finish()
-        
+
         image = sampler.aimage
         self.finalize_image(image)
         return image
 
         fields = [self.field]
         resolution = self.resolution
-        pf = self.pf
-        if self.weight is not None:
-            # This is a temporary field, which we will remove at the end.
-            def _make_wf(f, w):
-                def temp_weightfield(a, b):
-                    tr = b[f].astype("float64") * b[w]
-                    return tr
-                return temp_weightfield
-            pf.field_info.add_field("temp_weightfield",
-                function=_make_wf(self.field, self.weight))
-            fields = ["temp_weightfield", self.weight]
-        self.fields = fields
-        return Camera.snapshot(self, fn = fn, clip_ratio = clip_ratio, double_check = double_check,
-                 num_threads = num_threads)
 
+        image = self.new_image()
+
+        args = self.get_sampler_args(image)
+
+        sampler = self.get_sampler(args)
+
+        self.initialize_source()
+        
+        image = self._render(double_check, num_threads, image, sampler)
+
+        self.save_image(fn, clip_ratio, image)
+
+        return image
 
 data_object_registry["projection_camera"] = ProjectionCamera
 
+def off_axis_projection(pf, center, normal_vector, width, resolution,
+                        field, weight = None, num_threads = 0, 
+                        volume = None, no_ghost = False, interpolated = False):
+    r"""Project through a parameter file, off-axis, and return the image plane.
+
+    This function will accept the necessary items to integrate through a volume
+    at an arbitrary angle and return the integrated field of view to the user.
+    Note that if a weight is supplied, it will multiply the pre-interpolated
+    values together, then create cell-centered values, then interpolate within
+    the cell to conduct the integration.
+
+    Parameters
+    ----------
+    pf : `~yt.data_objects.api.StaticOutput`
+        This is the parameter file to volume render.
+    center : array_like
+        The current 'center' of the view port -- the focal point for the
+        camera.
+    normal_vector : array_like
+        The vector between the camera position and the center.
+    width : float or list of floats
+        The current width of the image.  If a single float, the volume is
+        cubical, but if not, it is left/right, top/bottom, front/back
+    resolution : int or list of ints
+        The number of pixels in each direction.
+    field : string
+        The field to project through the volume
+    weight : optional, default None
+        If supplied, the field will be pre-multiplied by this, then divided by
+        the integrated value of this field.  This returns an average rather
+        than a sum.
+    volume : `yt.extensions.volume_rendering.HomogenizedVolume`, optional
+        The volume to ray cast through.  Can be specified for finer-grained
+        control, but otherwise will be automatically generated.
+    no_ghost: bool, optional
+        Optimization option.  If True, homogenized bricks will
+        extrapolate out from grid instead of interpolating from
+        ghost zones that have to first be calculated.  This can
+        lead to large speed improvements, but at a loss of
+        accuracy/smoothness in resulting image.  The effects are
+        less notable when the transfer function is smooth and
+        broad. Default: True
+    interpolated : optional, default False
+        If True, the data is first interpolated to vertex-centered data, 
+        then tri-linearly interpolated along the ray. Not suggested for 
+        quantitative studies.
+
+    Returns
+    -------
+    image : array
+        An (N,N) array of the final integrated values, in float64 form.
+
+    Examples
+    --------
+
+    >>> image = off_axis_projection(pf, [0.5, 0.5, 0.5], [0.2,0.3,0.4],
+                      0.2, N, "Temperature", "Density")
+    >>> write_image(na.log10(image), "offaxis.png")
+
+    """
+    projcam = ProjectionCamera(center, normal_vector, width, resolution,
+            field, weight=weight, pf=pf, volume=volume,
+            no_ghost=no_ghost, interpolated=interpolated)
+    image = projcam.snapshot(num_threads=num_threads)
+    if weight is not None:
+        pf.field_info.pop("temp_weightfield")
+    del projcam
+    return image
+
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.