Changeset cf8b630 in sasmodels for explore/jitter.py


Ignore:
Timestamp:
Nov 2, 2017 3:50:13 PM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
ff10479
Parents:
a06af5d (diff), bd36af0 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into ticket-776-orientation

File:
1 edited

Legend:

Unmodified
Added
Removed

  • explore/jitter.py

    • Property mode changed from 100644 to 100755
    r85190c2 rde71632  
     1#!/usr/bin/env python 
    12""" 
    23Application to explore the difference between sasview 3.x orientation 
    34dispersity and possible replacement algorithms. 
    45""" 
     6from __future__ import division, print_function 
     7 
     8import sys, os 
     9sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) 
     10 
     11import argparse 
     12 
    513import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot 
    614import matplotlib.pyplot as plt 
    715from matplotlib.widgets import Slider, CheckButtons 
    816from matplotlib import cm 
    9  
    1017import numpy as np 
    1118from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    1219 
    13 def draw_beam(ax): 
     20def draw_beam(ax, view=(0, 0)): 
     21    """ 
     22    Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1) 
     23    """ 
    1424    #ax.plot([0,0],[0,0],[1,-1]) 
    1525    #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
    … …  
    2232    x = r*np.outer(np.cos(u), np.ones_like(v)) 
    2333    y = r*np.outer(np.sin(u), np.ones_like(v)) 
    24     z = np.outer(np.ones_like(u), v) 
     34    z = 1.3*np.outer(np.ones_like(u), v) 
     35 
     36    theta, phi = view 
     37    shape = x.shape 
     38    points = np.matrix([x.flatten(), y.flatten(), z.flatten()]) 
     39    points = Rz(phi)*Ry(theta)*points 
     40    x, y, z = [v.reshape(shape) for v in points] 
    2541 
    2642    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
    2743 
    28 def draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi): 
    29     size=[0.1, 0.4, 1.0] 
    30     view=[theta, phi, psi] 
    31     shimmy=[0,0,0] 
    32     #draw_shape = draw_parallelepiped 
    33     draw_shape = draw_ellipsoid 
     44def draw_jitter(ax, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0)): 
     45    """ 
     46    Represent jitter as a set of shapes at different orientations. 
     47    """ 
     48    # set max diagonal to 0.95 
     49    scale = 0.95/sqrt(sum(v**2 for v in size)) 
     50    size = tuple(scale*v for v in size) 
     51    draw_shape = draw_parallelepiped 
     52    #draw_shape = draw_ellipsoid 
    3453 
    3554    #np.random.seed(10) 
    … …  
    6483        [ 1,  1,  1], 
    6584    ] 
     85    dtheta, dphi, dpsi = jitter 
    6686    if dtheta == 0: 
    6787        cloud = [v for v in cloud if v[0] == 0] 
    … …  
    7090    if dpsi == 0: 
    7191        cloud = [v for v in cloud if v[2] == 0] 
    72     draw_shape(ax, size, view, shimmy, steps=100, alpha=0.8) 
     92    draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8) 
     93    scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist] 
    7394    for point in cloud: 
    74         shimmy=[dtheta*point[0], dphi*point[1], dpsi*point[2]] 
    75         draw_shape(ax, size, view, shimmy, alpha=0.8) 
     95        delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]] 
     96        draw_shape(ax, size, view, delta, alpha=0.8) 
    7697    for v in 'xyz': 
    7798        a, b, c = size 
    … …  
    80101        getattr(ax, v+'axis').label.set_text(v) 
    81102 
    82 def draw_ellipsoid(ax, size, view, shimmy, steps=25, alpha=1): 
     103def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 
     104    """Draw an ellipsoid.""" 
    83105    a,b,c = size 
    84     theta, phi, psi = view 
    85     dtheta, dphi, dpsi = shimmy 
    86  
    87106    u = np.linspace(0, 2 * np.pi, steps) 
    88107    v = np.linspace(0, np.pi, steps) 
    … …  
    90109    y = b*np.outer(np.sin(u), np.sin(v)) 
    91110    z = c*np.outer(np.ones_like(u), np.cos(v)) 
    92  
    93     shape = x.shape 
    94     points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
    95     points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 
    96     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
    97     x,y,z = [v.reshape(shape) for v in points] 
     111    x, y, z = transform_xyz(view, jitter, x, y, z) 
    98112 
    99113    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
    100114 
    101 def draw_parallelepiped(ax, size, view, shimmy, alpha=1): 
     115    draw_labels(ax, view, jitter, [ 
     116         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     117         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     118         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     119         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     120         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     121         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     122    ]) 
     123 
     124def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 
     125    """Draw a parallelepiped.""" 
    102126    a,b,c = size 
    103     theta, phi, psi = view 
    104     dtheta, dphi, dpsi = shimmy 
    105  
    106127    x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    107128    y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
    … …  
    118139    ]) 
    119140 
    120     points = np.matrix([x,y,z]) 
    121     points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 
    122     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
    123  
    124     x,y,z = [np.array(v).flatten() for v in points] 
     141    x, y, z = transform_xyz(view, jitter, x, y, z) 
    125142    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    126143 
     144    draw_labels(ax, view, jitter, [ 
     145         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     146         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     147         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     148         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     149         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     150         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     151    ]) 
     152 
    127153def draw_sphere(ax, radius=10., steps=100): 
     154    """Draw a sphere""" 
    128155    u = np.linspace(0, 2 * np.pi, steps) 
    129156    v = np.linspace(0, np.pi, steps) 
    … …  
    134161    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    135162 
    136 def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'): 
    137     theta_center = radians(theta) 
    138     phi_center = radians(phi) 
    139     flow_center = radians(flow) 
    140     dtheta = radians(dtheta) 
    141     dphi = radians(dphi) 
    142  
    143     # 10 point 3-sigma gaussian weights 
    144     t = np.linspace(-3., 3., 11) 
    145     if dist == 'gauss': 
     163PROJECTIONS = [ 
     164    'equirectangular', 'azimuthal_equidistance', 'sinusoidal', 'guyou', 
     165] 
     166def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian', 
     167              projection='equirectangular'): 
     168    """ 
     169    Draw the dispersion mesh showing the theta-phi orientations at which 
     170    the model will be evaluated. 
     171 
     172    jitter projections 
     173    <https://en.wikipedia.org/wiki/List_of_map_projections> 
     174 
     175    equirectangular (standard latitude-longitude mesh) 
     176        <https://en.wikipedia.org/wiki/Equirectangular_projection> 
     177        Allows free movement in phi (around the equator), but theta is 
     178        limited to +/- 90, and points are cos-weighted. Jitter in phi is 
     179        uniform in weight along a line of latitude.  With small theta and 
     180        phi ranging over +/- 180 this forms a wobbling disk.  With small 
     181        phi and theta ranging over +/- 90 this forms a wedge like a slice 
     182        of an orange. 
     183    azimuthal_equidistance (Postel) 
     184        <https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection> 
     185        Preserves distance from center, and so is an excellent map for 
     186        representing a bivariate gaussian on the surface.  Theta and phi 
     187        operate identically, cutting wegdes from the antipode of the viewing 
     188        angle.  This unfortunately does not allow free movement in either 
     189        theta or phi since the orthogonal wobble decreases to 0 as the body 
     190        rotates through 180 degrees. 
     191    sinusoidal (Sanson-Flamsteed, Mercator equal-area) 
     192        <https://en.wikipedia.org/wiki/Sinusoidal_projection> 
     193        Preserves arc length with latitude, giving bad behaviour at 
     194        theta near +/- 90.  Theta and phi operate somewhat differently, 
     195        so a system with a-b-c dtheta-dphi-dpsi will not give the same 
     196        value as one with b-a-c dphi-dtheta-dpsi, as would be the case 
     197        for azimuthal equidistance.  Free movement using theta or phi 
     198        uniform over +/- 180 will work, but not as well as equirectangular 
     199        phi, with theta being slightly worse.  Computationally it is much 
     200        cheaper for wide theta-phi meshes since it excludes points which 
     201        lie outside the sinusoid near theta +/- 90 rather than packing 
     202        them close together as in equirectangle.  Note that the poles 
     203        will be slightly overweighted for theta > 90 with the circle 
     204        from theta at 90+dt winding backwards around the pole, overlapping 
     205        the circle from theta at 90-dt. 
     206    Guyou (hemisphere-in-a-square) **not weighted** 
     207        <https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection> 
     208        With tiling, allows rotation in phi or theta through +/- 180, with 
     209        uniform spacing.  Both theta and phi allow free rotation, with wobble 
     210        in the orthogonal direction reasonably well behaved (though not as 
     211        good as equirectangular phi). The forward/reverse transformations 
     212        relies on elliptic integrals that are somewhat expensive, so the 
     213        behaviour has to be very good to justify the cost and complexity. 
     214        The weighting function for each point has not yet been computed. 
     215        Note: run the module *guyou.py* directly and it will show the forward 
     216        and reverse mappings. 
     217    azimuthal_equal_area  **incomplete** 
     218        <https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection> 
     219        Preserves the relative density of the surface patches.  Not that 
     220        useful and not completely implemented 
     221    Gauss-Kreuger **not implemented** 
     222        <https://en.wikipedia.org/wiki/Transverse_Mercator_projection#Ellipsoidal_transverse_Mercator> 
     223        Should allow free movement in theta, but phi is distorted. 
     224    """ 
     225    theta, phi, psi = view 
     226    dtheta, dphi, dpsi = jitter 
     227 
     228    t = np.linspace(-1, 1, n) 
     229    weights = np.ones_like(t) 
     230    if dist == 'gaussian': 
     231        t *= 3 
    146232        weights = exp(-0.5*t**2) 
    147     elif dist == 'rect': 
    148         weights = np.ones_like(t) 
     233    elif dist == 'rectangle': 
     234        # Note: uses sasmodels ridiculous definition of rectangle width 
     235        t *= sqrt(3) 
     236    elif dist == 'uniform': 
     237        pass 
    149238    else: 
    150         raise ValueError("expected dist to be 'gauss' or 'rect'") 
    151     theta = dtheta*t 
    152     phi = dphi*t 
    153  
    154     x = radius * np.outer(cos(phi), cos(theta)) 
    155     y = radius * np.outer(sin(phi), cos(theta)) 
    156     z = radius * np.outer(np.ones_like(phi), sin(theta)) 
    157     #w = np.outer(weights, weights*abs(cos(dtheta*t))) 
    158     w = np.outer(weights, weights*abs(cos(theta))) 
    159  
    160     x, y, z, w = [v.flatten() for v in (x,y,z,w)] 
    161     x, y, z = rotate(x, y, z, phi_center, theta_center, flow_center) 
    162  
     239        raise ValueError("expected dist to be gaussian, rectangle or uniform") 
     240 
     241    if projection == 'equirectangular': 
     242        def rotate(theta_i, phi_j): 
     243            return Rx(phi_j)*Ry(theta_i) 
     244        def weight(theta_i, phi_j, wi, wj): 
     245            return wi*wj*abs(cos(radians(theta_i))) 
     246    elif projection == 'azimuthal_equidistance': 
     247        def rotate(theta_i, phi_j): 
     248            latitude = sqrt(theta_i**2 + phi_j**2) 
     249            longitude = degrees(np.arctan2(phi_j, theta_i)) 
     250            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
     251            return Rz(longitude)*Ry(latitude) 
     252        def weight(theta_i, phi_j, wi, wj): 
     253            # Weighting for each point comes from the integral: 
     254            #     \int\int I(q, lat, log) sin(lat) dlat dlog 
     255            # We are doing a conformal mapping from disk to sphere, so we need 
     256            # a change of variables g(theta, phi) -> (lat, long): 
     257            #     lat, long = sqrt(theta^2 + phi^2), arctan(phi/theta) 
     258            # giving: 
     259            #     dtheta dphi = det(J) dlat dlong 
     260            # where J is the jacobian from the partials of g. Using 
     261            #     R = sqrt(theta^2 + phi^2), 
     262            # then 
     263            #     J = [[x/R, Y/R], -y/R^2, x/R^2]] 
     264            # and 
     265            #     det(J) = 1/R 
     266            # with the final integral being: 
     267            #    \int\int I(q, theta, phi) sin(R)/R dtheta dphi 
     268            # 
     269            # This does approximately the right thing, decreasing the weight 
     270            # of each point as you go farther out on the disk, but it hasn't 
     271            # yet been checked against the 1D integral results. Prior 
     272            # to declaring this "good enough" and checking that integrals 
     273            # work in practice, we will examine alternative mappings. 
     274            # 
     275            # The issue is that the mapping does not support the case of free 
     276            # rotation about a single axis correctly, with a small deviation 
     277            # in the orthogonal axis independent of the first axis.  Like the 
     278            # usual polar coordiates integration, the integrated sections 
     279            # form wedges, though at least in this case the wedge cuts through 
     280            # the entire sphere, and treats theta and phi identically. 
     281            latitude = sqrt(theta_i**2 + phi_j**2) 
     282            w = sin(radians(latitude))/latitude if latitude != 0 else 1 
     283            return w*wi*wj if latitude < 180 else 0 
     284    elif projection == 'sinusoidal': 
     285        def rotate(theta_i, phi_j): 
     286            latitude = theta_i 
     287            scale = cos(radians(latitude)) 
     288            longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0 
     289            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
     290            return Rx(longitude)*Ry(latitude) 
     291        def weight(theta_i, phi_j, wi, wj): 
     292            latitude = theta_i 
     293            scale = cos(radians(latitude)) 
     294            w = 1 if abs(phi_j) < abs(scale)*180 else 0 
     295            return w*wi*wj 
     296    elif projection == 'azimuthal_equal_area': 
     297        def rotate(theta_i, phi_j): 
     298            R = min(1, sqrt(theta_i**2 + phi_j**2)/180) 
     299            latitude = 180-degrees(2*np.arccos(R)) 
     300            longitude = degrees(np.arctan2(phi_j, theta_i)) 
     301            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
     302            return Rz(longitude)*Ry(latitude) 
     303        def weight(theta_i, phi_j, wi, wj): 
     304            latitude = sqrt(theta_i**2 + phi_j**2) 
     305            w = sin(radians(latitude))/latitude if latitude != 0 else 1 
     306            return w*wi*wj if latitude < 180 else 0 
     307    elif projection == 'guyou': 
     308        def rotate(theta_i, phi_j): 
     309            from guyou import guyou_invert 
     310            #latitude, longitude = guyou_invert([theta_i], [phi_j]) 
     311            longitude, latitude = guyou_invert([phi_j], [theta_i]) 
     312            return Rx(longitude[0])*Ry(latitude[0]) 
     313        def weight(theta_i, phi_j, wi, wj): 
     314            return wi*wj 
     315    else: 
     316        raise ValueError("unknown projection %r"%projection) 
     317 
     318    # mesh in theta, phi formed by rotating z 
     319    z = np.matrix([[0], [0], [radius]]) 
     320    points = np.hstack([rotate(theta_i, phi_j)*z 
     321                        for theta_i in dtheta*t 
     322                        for phi_j in dphi*t]) 
     323    # select just the active points (i.e., those with phi < 180 
     324    w = np.array([weight(theta_i, phi_j, wi, wj) 
     325                  for wi, theta_i in zip(weights, dtheta*t) 
     326                  for wj, phi_j in zip(weights, dphi*t)]) 
     327    #print(max(w), min(w), min(w[w>0])) 
     328    points = points[:, w>0] 
     329    w = w[w>0] 
     330    w /= max(w) 
     331 
     332    if 0: # Kent distribution 
     333        points = np.hstack([Rx(phi_j)*Ry(theta_i)*z for theta_i in 30*t for phi_j in 60*t]) 
     334        xp, yp, zp = [np.array(v).flatten() for v in points] 
     335        kappa = max(1e6, radians(dtheta)/(2*pi)) 
     336        beta = 1/max(1e-6, radians(dphi)/(2*pi))/kappa 
     337        w = exp(kappa*zp) #+ beta*(xp**2 + yp**2) 
     338        print(kappa, dtheta, radians(dtheta), min(w), max(w), sum(w)) 
     339        #w /= abs(cos(radians( 
     340        #w /= sum(w) 
     341 
     342    # rotate relative to beam 
     343    points = orient_relative_to_beam(view, points) 
     344 
     345    x, y, z = [np.array(v).flatten() for v in points] 
     346    #plt.figure(2); plt.clf(); plt.hist(z, bins=np.linspace(-1, 1, 51)) 
    163347    ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.) 
    164348 
    165 def rotate(x, y, z, phi, theta, psi): 
    166     R = Rz(phi)*Ry(theta)*Rz(psi) 
    167     p = np.vstack([x,y,z]) 
    168     return R*p 
    169  
     349def draw_labels(ax, view, jitter, text): 
     350    """ 
     351    Draw text at a particular location. 
     352    """ 
     353    labels, locations, orientations = zip(*text) 
     354    px, py, pz = zip(*locations) 
     355    dx, dy, dz = zip(*orientations) 
     356 
     357    px, py, pz = transform_xyz(view, jitter, px, py, pz) 
     358    dx, dy, dz = transform_xyz(view, jitter, dx, dy, dz) 
     359 
     360    # TODO: zdir for labels is broken, and labels aren't appearing. 
     361    for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)): 
     362        zdir = np.asarray(zdir).flatten() 
     363        ax.text(p[0], p[1], p[2], label, zdir=zdir) 
     364 
     365# Definition of rotation matrices comes from wikipedia: 
     366#    https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations 
    170367def Rx(angle): 
     368    """Construct a matrix to rotate points about *x* by *angle* degrees.""" 
    171369    a = radians(angle) 
    172     R = [[1., 0., 0.], 
    173          [0.,  cos(a), sin(a)], 
    174          [0., -sin(a), cos(a)]] 
     370    R = [[1, 0, 0], 
     371         [0, +cos(a), -sin(a)], 
     372         [0, +sin(a), +cos(a)]] 
    175373    return np.matrix(R) 
    176374 
    177375def Ry(angle): 
     376    """Construct a matrix to rotate points about *y* by *angle* degrees.""" 
    178377    a = radians(angle) 
    179     R = [[cos(a), 0., -sin(a)], 
    180          [0., 1., 0.], 
    181          [sin(a), 0.,  cos(a)]] 
     378    R = [[+cos(a), 0, +sin(a)], 
     379         [0, 1, 0], 
     380         [-sin(a), 0, +cos(a)]] 
    182381    return np.matrix(R) 
    183382 
    184383def Rz(angle): 
     384    """Construct a matrix to rotate points about *z* by *angle* degrees.""" 
    185385    a = radians(angle) 
    186     R = [[cos(a), -sin(a), 0.], 
    187          [sin(a),  cos(a), 0.], 
    188          [0., 0., 1.]] 
     386    R = [[+cos(a), -sin(a), 0], 
     387         [+sin(a), +cos(a), 0], 
     388         [0, 0, 1]] 
    189389    return np.matrix(R) 
    190390 
    191 def main(): 
     391def transform_xyz(view, jitter, x, y, z): 
     392    """ 
     393    Send a set of (x,y,z) points through the jitter and view transforms. 
     394    """ 
     395    x, y, z = [np.asarray(v) for v in (x, y, z)] 
     396    shape = x.shape 
     397    points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
     398    points = apply_jitter(jitter, points) 
     399    points = orient_relative_to_beam(view, points) 
     400    x, y, z = [np.array(v).reshape(shape) for v in points] 
     401    return x, y, z 
     402 
     403def apply_jitter(jitter, points): 
     404    """ 
     405    Apply the jitter transform to a set of points. 
     406 
     407    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
     408    """ 
     409    dtheta, dphi, dpsi = jitter 
     410    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
     411    return points 
     412 
     413def orient_relative_to_beam(view, points): 
     414    """ 
     415    Apply the view transform to a set of points. 
     416 
     417    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
     418    """ 
     419    theta, phi, psi = view 
     420    points = Rz(phi)*Ry(theta)*Rz(psi)*points 
     421    return points 
     422 
     423# translate between number of dimension of dispersity and the number of 
     424# points along each dimension. 
     425PD_N_TABLE = { 
     426    (0, 0, 0): (0, 0, 0),     # 0 
     427    (1, 0, 0): (100, 0, 0),   # 100 
     428    (0, 1, 0): (0, 100, 0), 
     429    (0, 0, 1): (0, 0, 100), 
     430    (1, 1, 0): (30, 30, 0),   # 900 
     431    (1, 0, 1): (30, 0, 30), 
     432    (0, 1, 1): (0, 30, 30), 
     433    (1, 1, 1): (15, 15, 15),  # 3375 
     434} 
     435 
     436def clipped_range(data, portion=1.0, mode='central'): 
     437    """ 
     438    Determine range from data. 
     439 
     440    If *portion* is 1, use full range, otherwise use the center of the range 
     441    or the top of the range, depending on whether *mode* is 'central' or 'top'. 
     442    """ 
     443    if portion == 1.0: 
     444        return data.min(), data.max() 
     445    elif mode == 'central': 
     446        data = np.sort(data.flatten()) 
     447        offset = int(portion*len(data)/2 + 0.5) 
     448        return data[offset], data[-offset] 
     449    elif mode == 'top': 
     450        data = np.sort(data.flatten()) 
     451        offset = int(portion*len(data) + 0.5) 
     452        return data[offset], data[-1] 
     453 
     454def draw_scattering(calculator, ax, view, jitter, dist='gaussian'): 
     455    """ 
     456    Plot the scattering for the particular view. 
     457 
     458    *calculator* is returned from :func:`build_model`.  *ax* are the 3D axes 
     459    on which the data will be plotted.  *view* and *jitter* are the current 
     460    orientation and orientation dispersity.  *dist* is one of the sasmodels 
     461    weight distributions. 
     462    """ 
     463    if dist == 'uniform':  # uniform is not yet in this branch 
     464        dist, scale = 'rectangle', 1/sqrt(3) 
     465    else: 
     466        scale = 1 
     467 
     468    # add the orientation parameters to the model parameters 
     469    theta, phi, psi = view 
     470    theta_pd, phi_pd, psi_pd = [scale*v for v in jitter] 
     471    theta_pd_n, phi_pd_n, psi_pd_n = PD_N_TABLE[(theta_pd>0, phi_pd>0, psi_pd>0)] 
     472    ## increase pd_n for testing jitter integration rather than simple viz 
     473    #theta_pd_n, phi_pd_n, psi_pd_n = [5*v for v in (theta_pd_n, phi_pd_n, psi_pd_n)] 
     474 
     475    pars = dict( 
     476        theta=theta, theta_pd=theta_pd, theta_pd_type=dist, theta_pd_n=theta_pd_n, 
     477        phi=phi, phi_pd=phi_pd, phi_pd_type=dist, phi_pd_n=phi_pd_n, 
     478        psi=psi, psi_pd=psi_pd, psi_pd_type=dist, psi_pd_n=psi_pd_n, 
     479    ) 
     480    pars.update(calculator.pars) 
     481 
     482    # compute the pattern 
     483    qx, qy = calculator._data.x_bins, calculator._data.y_bins 
     484    Iqxy = calculator(**pars).reshape(len(qx), len(qy)) 
     485 
     486    # scale it and draw it 
     487    Iqxy = np.log(Iqxy) 
     488    if calculator.limits: 
     489        # use limits from orientation (0,0,0) 
     490        vmin, vmax = calculator.limits 
     491    else: 
     492        vmin, vmax = clipped_range(Iqxy, portion=0.95, mode='top') 
     493    #print("range",(vmin,vmax)) 
     494    #qx, qy = np.meshgrid(qx, qy) 
     495    if 0: 
     496        level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i') 
     497        level[level<0] = 0 
     498        colors = plt.get_cmap()(level) 
     499        ax.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) 
     500    elif 1: 
     501        ax.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, 
     502                    levels=np.linspace(vmin, vmax, 24)) 
     503    else: 
     504        ax.pcolormesh(qx, qy, Iqxy) 
     505 
     506def build_model(model_name, n=150, qmax=0.5, **pars): 
     507    """ 
     508    Build a calculator for the given shape. 
     509 
     510    *model_name* is any sasmodels model.  *n* and *qmax* define an n x n mesh 
     511    on which to evaluate the model.  The remaining parameters are stored in 
     512    the returned calculator as *calculator.pars*.  They are used by 
     513    :func:`draw_scattering` to set the non-orientation parameters in the 
     514    calculation. 
     515 
     516    Returns a *calculator* function which takes a dictionary or parameters and 
     517    produces Iqxy.  The Iqxy value needs to be reshaped to an n x n matrix 
     518    for plotting.  See the :class:`sasmodels.direct_model.DirectModel` class 
     519    for details. 
     520    """ 
     521    from sasmodels.core import load_model_info, build_model 
     522    from sasmodels.data import empty_data2D 
     523    from sasmodels.direct_model import DirectModel 
     524 
     525    model_info = load_model_info(model_name) 
     526    model = build_model(model_info) #, dtype='double!') 
     527    q = np.linspace(-qmax, qmax, n) 
     528    data = empty_data2D(q, q) 
     529    calculator = DirectModel(data, model) 
     530 
     531    # stuff the values for non-orientation parameters into the calculator 
     532    calculator.pars = pars.copy() 
     533    calculator.pars.setdefault('backgound', 1e-3) 
     534 
     535    # fix the data limits so that we can see if the pattern fades 
     536    # under rotation or angular dispersion 
     537    Iqxy = calculator(theta=0, phi=0, psi=0, **calculator.pars) 
     538    Iqxy = np.log(Iqxy) 
     539    vmin, vmax = clipped_range(Iqxy, 0.95, mode='top') 
     540    calculator.limits = vmin, vmax+1 
     541 
     542    return calculator 
     543 
     544def select_calculator(model_name, n=150, size=(10,40,100)): 
     545    """ 
     546    Create a model calculator for the given shape. 
     547 
     548    *model_name* is one of sphere, cylinder, ellipsoid, triaxial_ellipsoid, 
     549    parallelepiped or bcc_paracrystal. *n* is the number of points to use 
     550    in the q range.  *qmax* is chosen based on model parameters for the 
     551    given model to show something intersting. 
     552 
     553    Returns *calculator* and tuple *size* (a,b,c) giving minor and major 
     554    equitorial axes and polar axis respectively.  See :func:`build_model` 
     555    for details on the returned calculator. 
     556    """ 
     557    a, b, c = size 
     558    if model_name == 'sphere': 
     559        calculator = build_model('sphere', n=n, radius=c) 
     560        a = b = c 
     561    elif model_name == 'bcc_paracrystal': 
     562        calculator = build_model('bcc_paracrystal', n=n, dnn=c, 
     563                                  d_factor=0.06, radius=40) 
     564        a = b = c 
     565    elif model_name == 'cylinder': 
     566        calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c) 
     567        a = b 
     568    elif model_name == 'ellipsoid': 
     569        calculator = build_model('ellipsoid', n=n, qmax=1.0, 
     570                                 radius_polar=c, radius_equatorial=b) 
     571        a = b 
     572    elif model_name == 'triaxial_ellipsoid': 
     573        calculator = build_model('triaxial_ellipsoid', n=n, qmax=0.5, 
     574                                 radius_equat_minor=a, 
     575                                 radius_equat_major=b, 
     576                                 radius_polar=c) 
     577    elif model_name == 'parallelepiped': 
     578        calculator = build_model('parallelepiped', n=n, a=a, b=b, c=c) 
     579    else: 
     580        raise ValueError("unknown model %s"%model_name) 
     581 
     582    return calculator, (a, b, c) 
     583 
     584SHAPES = [ 
     585    'parallelepiped', 'triaxial_ellipsoid', 'bcc_paracrystal', 
     586    'cylinder', 'ellipsoid', 
     587    'sphere', 
     588 ] 
     589 
     590DISTRIBUTIONS = [ 
     591    'gaussian', 'rectangle', 'uniform', 
     592] 
     593DIST_LIMITS = { 
     594    'gaussian': 30, 
     595    'rectangle': 90/sqrt(3), 
     596    'uniform': 90, 
     597} 
     598 
     599def run(model_name='parallelepiped', size=(10, 40, 100), 
     600        dist='gaussian', mesh=30, 
     601        projection='equirectangular'): 
     602    """ 
     603    Show an interactive orientation and jitter demo. 
     604 
     605    *model_name* is one of the models available in :func:`select_model`. 
     606    """ 
     607    # set up calculator 
     608    calculator, size = select_calculator(model_name, n=150, size=size) 
     609 
     610    ## uncomment to set an independent the colour range for every view 
     611    ## If left commented, the colour range is fixed for all views 
     612    calculator.limits = None 
     613 
     614    ## initial view 
     615    #theta, dtheta = 70., 10. 
     616    #phi, dphi = -45., 3. 
     617    #psi, dpsi = -45., 3. 
     618    theta, phi, psi = 0, 0, 0 
     619    dtheta, dphi, dpsi = 0, 0, 0 
     620 
     621    ## create the plot window 
    192622    #plt.hold(True) 
     623    plt.subplots(num=None, figsize=(5.5, 5.5)) 
    193624    plt.set_cmap('gist_earth') 
    194625    plt.clf() 
     626    plt.gcf().canvas.set_window_title(projection) 
    195627    #gs = gridspec.GridSpec(2,1,height_ratios=[4,1]) 
    196628    #ax = plt.subplot(gs[0], projection='3d') 
    197629    ax = plt.axes([0.0, 0.2, 1.0, 0.8], projection='3d') 
    198  
    199     theta, dtheta = 70., 10. 
    200     phi, dphi = -45., 3. 
    201     psi, dpsi = -45., 3. 
    202     theta, phi, psi = 0, 0, 0 
    203     dtheta, dphi, dpsi = 0, 0, 0 
    204     #dist = 'rect' 
    205     dist = 'gauss' 
     630    try:  # CRUFT: not all versions of matplotlib accept 'square' 3d projection 
     631        ax.axis('square') 
     632    except Exception: 
     633        pass 
    206634 
    207635    axcolor = 'lightgoldenrodyellow' 
     636 
     637    ## add control widgets to plot 
    208638    axtheta  = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    209639    axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
    … …  
    212642    sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi) 
    213643    spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi) 
     644 
    214645    axdtheta  = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    215646    axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
    216647    axdpsi= plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) 
    217     sdtheta = Slider(axdtheta, 'dTheta', 0, 30, valinit=dtheta) 
    218     sdphi = Slider(axdphi, 'dPhi', 0, 30, valinit=dphi) 
    219     sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dphi) 
    220  
     648    # Note: using ridiculous definition of rectangle distribution, whose width 
     649    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 
     650    # the maximum width to 90. 
     651    dlimit = DIST_LIMITS[dist] 
     652    sdtheta = Slider(axdtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta) 
     653    sdphi = Slider(axdphi, 'dPhi', 0, 2*dlimit, valinit=dphi) 
     654    sdpsi = Slider(axdpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi) 
     655 
     656    ## callback to draw the new view 
    221657    def update(val, axis=None): 
    222         theta, phi, psi = stheta.val, sphi.val, spsi.val 
    223         dtheta, dphi, dpsi = sdtheta.val, sdphi.val, sdpsi.val 
     658        view = stheta.val, sphi.val, spsi.val 
     659        jitter = sdtheta.val, sdphi.val, sdpsi.val 
     660        # set small jitter as 0 if multiple pd dims 
     661        dims = sum(v > 0 for v in jitter) 
     662        limit = [0, 0, 0.5, 5][dims] 
     663        jitter = [0 if v < limit else v for v in jitter] 
    224664        ax.cla() 
    225         draw_beam(ax) 
    226         draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi) 
    227         #if not axis.startswith('d'): 
    228         #    ax.view_init(elev=theta, azim=phi) 
     665        draw_beam(ax, (0, 0)) 
     666        draw_jitter(ax, view, jitter, dist=dist, size=size) 
     667        #draw_jitter(ax, view, (0,0,0)) 
     668        draw_mesh(ax, view, jitter, dist=dist, n=mesh, projection=projection) 
     669        draw_scattering(calculator, ax, view, jitter, dist=dist) 
    229670        plt.gcf().canvas.draw() 
    230671 
     672    ## bind control widgets to view updater 
    231673    stheta.on_changed(lambda v: update(v,'theta')) 
    232674    sphi.on_changed(lambda v: update(v, 'phi')) 
    … …  
    236678    sdpsi.on_changed(lambda v: update(v, 'dpsi')) 
    237679 
     680    ## initialize view 
    238681    update(None, 'phi') 
    239682 
     683    ## go interactive 
    240684    plt.show() 
     685 
     686def main(): 
     687    parser = argparse.ArgumentParser( 
     688        description="Display jitter", 
     689        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
     690        ) 
     691    parser.add_argument('-p', '--projection', choices=PROJECTIONS, default=PROJECTIONS[0], help='coordinate projection') 
     692    parser.add_argument('-s', '--size', type=str, default='10,40,100', help='a,b,c lengths') 
     693    parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, default=DISTRIBUTIONS[0], help='jitter distribution') 
     694    parser.add_argument('-m', '--mesh', type=int, default=30, help='#points in theta-phi mesh') 
     695    parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], help='oriented shape') 
     696    opts = parser.parse_args() 
     697    size = tuple(int(v) for v in opts.size.split(',')) 
     698    run(opts.shape, size=size, 
     699        mesh=opts.mesh, dist=opts.distribution, 
     700        projection=opts.projection) 
    241701 
    242702if __name__ == "__main__": 
Note: See TracChangeset for help on using the changeset viewer.