Changeset bcb5594 in sasmodels


Ignore:
Timestamp:
Oct 31, 2017 1:54:21 PM (6 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:
d73a5ac
Parents:
904cd9c
Message:

allow selection of jitter mesh behaviour from command line

File:
1 edited

Legend:

Unmodified
Added
Removed
  • explore/jitter.py

    r5b5ea20 rbcb5594  
    88import sys, os 
    99sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) 
     10 
     11import argparse 
    1012 
    1113import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot 
     
    8991        cloud = [v for v in cloud if v[2] == 0] 
    9092    draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8) 
    91     scale = 1/sqrt(3) if dist == 'rectangle' else 1 
     93    scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist] 
    9294    for point in cloud: 
    9395        delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]] 
     
    159161    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    160162 
    161 def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian'): 
     163PROJECTIONS = [ 
     164    'equirectangular', 'azimuthal_equidistance', 'sinusoidal', 
     165] 
     166def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian', 
     167              projection='equirectangular'): 
    162168    """ 
    163169    Draw the dispersion mesh showing the theta-phi orientations at which 
    164170    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. 
     203    Guyour (hemisphere-in-a-square)  **not implemented** 
     204        <https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection> 
     205        Promising.  With tiling should allow rotation in phi or theta 
     206        through +/- 180, preserving almost disk-like behaviour in either 
     207        direction (phi rotation will not be as uniform as it is in 
     208        equirectangular; not sure about theta).  Unfortunately, distortion 
     209        is not restricted to the corners of the theta-phi mesh, so this will 
     210        not be as good as the azimuthal equidistance project for gaussian 
     211        distributions. 
     212    azimuthal_equal_area  **incomplete** 
     213        <https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection> 
     214        Preserves the relative density of the surface patches.  Not that 
     215        useful and not completely implemented 
     216    Gauss-Kreuger **not implemented** 
     217        <https://en.wikipedia.org/wiki/Transverse_Mercator_projection#Ellipsoidal_transverse_Mercator> 
     218        Should allow free movement in theta, but phi is distorted. 
    165219    """ 
    166220    theta, phi, psi = view 
    167221    dtheta, dphi, dpsi = jitter 
    168222 
     223    t = np.linspace(-1, 1, n) 
     224    weights = np.ones_like(t) 
    169225    if dist == 'gaussian': 
    170         t = np.linspace(-3, 3, n) 
     226        t *= 3 
    171227        weights = exp(-0.5*t**2) 
    172228    elif dist == 'rectangle': 
    173229        # Note: uses sasmodels ridiculous definition of rectangle width 
    174         t = np.linspace(-1, 1, n)*sqrt(3) 
    175         weights = np.ones_like(t) 
     230        t *= sqrt(3) 
     231    elif dist == 'uniform': 
     232        pass 
    176233    else: 
    177         raise ValueError("expected dist to be 'gaussian' or 'rectangle'") 
    178  
    179     #PROJECTION = 'stretched_phi' 
    180     PROJECTION = 'azimuthal_equidistance' 
    181     #PROJECTION = 'azimuthal_equal_area' 
    182     if PROJECTION == 'stretced_phi': 
     234        raise ValueError("expected dist to be gaussian, rectangle or uniform") 
     235 
     236    if projection == 'equirectangular': 
    183237        def rotate(theta_i, phi_j): 
    184             if theta_i != 90: 
    185                 phi_j /= cos(radians(theta_i)) 
    186238            return Rx(phi_j)*Ry(theta_i) 
    187239        def weight(theta_i, phi_j, wi, wj): 
    188             if theta_i != 90: 
    189                 phi_j /= cos(radians(theta_i)) 
    190             return wi*wj if abs(phi_j) < 180 else 0 
    191     elif PROJECTION == 'azimuthal_equidistance': 
    192         # https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection 
     240            return wi*wj*cos(radians(theta_i)) 
     241    elif projection == 'azimuthal_equidistance': 
    193242        def rotate(theta_i, phi_j): 
    194243            latitude = sqrt(theta_i**2 + phi_j**2) 
     
    228277            w = sin(radians(latitude))/latitude if latitude != 0 else 1 
    229278            return w*wi*wj if latitude < 180 else 0 
    230     elif PROJECTION == 'azimuthal_equal_area': 
    231         # https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection 
     279    elif projection == 'sinusoidal': 
     280        def rotate(theta_i, phi_j): 
     281            latitude = theta_i 
     282            scale = cos(radians(latitude)) 
     283            longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0 
     284            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
     285            return Rx(longitude)*Ry(latitude) 
     286        def weight(theta_i, phi_j, wi, wj): 
     287            latitude = theta_i 
     288            scale = cos(radians(latitude)) 
     289            w = 1 if abs(phi_j) < abs(scale)*180 else 0 
     290            return w*wi*wj 
     291    elif projection == 'azimuthal_equal_area': 
    232292        def rotate(theta_i, phi_j): 
    233293            R = min(1, sqrt(theta_i**2 + phi_j**2)/180) 
     
    240300            w = sin(radians(latitude))/latitude if latitude != 0 else 1 
    241301            return w*wi*wj if latitude < 180 else 0 
    242     elif SCALED_PHI == 10:  # random thrashing 
    243         def rotate(theta_i, phi_j): 
    244             theta_i, phi_j = 2*theta_i/abs(cos(radians(phi_j))), 2*phi_j/cos(radians(theta_i)) 
    245             return Rx(phi_j)*Ry(theta_i) 
    246         def weight(theta_i, phi_j, wi, wj): 
    247             theta_i, phi_j = 2*theta_i/abs(cos(radians(phi_j))), 2*phi_j/cos(radians(theta_i)) 
    248             return wi*wj if abs(phi_j) < 180 else 0 
    249302    else: 
    250         def rotate(theta_i, phi_j): 
    251             return Rx(phi_j)*Ry(theta_i) 
    252         def weight(theta_i, phi_j, wi, wj): 
    253             return wi*wj*cos(radians(theta_i)) 
     303        raise ValueError("unknown projection %r"%projection) 
    254304 
    255305    # mesh in theta, phi formed by rotating z 
     
    398448    weight distributions. 
    399449    """ 
    400     ## Sasmodels use sqrt(3)*width for the rectangle range; scale to the 
    401     ## proper width for comparison. Commented out since now using the 
    402     ## sasmodels definition of width for rectangle. 
    403     #scale = 1/sqrt(3) if dist == 'rectangle' else 1 
    404     scale = 1 
     450    if dist == 'uniform':  # uniform is not yet in this branch 
     451        dist, scale = 'rectangle', 1/sqrt(3) 
     452    else: 
     453        scale = 1 
    405454 
    406455    # add the orientation parameters to the model parameters 
     
    520569    return calculator, (a, b, c) 
    521570 
    522 def main(model_name='parallelepiped', size=(10, 40, 100)): 
     571SHAPES = [ 
     572    'parallelepiped', 'triaxial_ellipsoid', 'bcc_paracrystal', 
     573    'cylinder', 'ellipsoid', 
     574    'sphere', 
     575 ] 
     576 
     577DISTRIBUTIONS = [ 
     578    'gaussian', 'rectangle', 'uniform', 
     579] 
     580DIST_LIMITS = { 
     581    'gaussian': 30, 
     582    'rectangle': 90/sqrt(3), 
     583    'uniform': 90, 
     584} 
     585 
     586def run(model_name='parallelepiped', size=(10, 40, 100), 
     587        dist='gaussian', mesh=30, 
     588        projection='equirectangular'): 
    523589    """ 
    524590    Show an interactive orientation and jitter demo. 
     
    532598    ## If left commented, the colour range is fixed for all views 
    533599    calculator.limits = None 
    534  
    535     ## use gaussian distribution unless testing integration 
    536     dist = 'rectangle' 
    537     #dist = 'gaussian' 
    538600 
    539601    ## initial view 
     
    572634    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 
    573635    # the maximum width to 90. 
    574     dlimit = 30 if dist == 'gaussian' else 90/sqrt(3) 
     636    dlimit = DIST_LIMITS[dist] 
    575637    sdtheta = Slider(axdtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta) 
    576638    sdphi = Slider(axdphi, 'dPhi', 0, 2*dlimit, valinit=dphi) 
     
    583645        # set small jitter as 0 if multiple pd dims 
    584646        dims = sum(v > 0 for v in jitter) 
    585         limit = [0, 0, 2, 5][dims] 
     647        limit = [0, 0, 0.5, 5][dims] 
    586648        jitter = [0 if v < limit else v for v in jitter] 
    587649        ax.cla() 
     
    589651        draw_jitter(ax, view, jitter, dist=dist, size=size) 
    590652        #draw_jitter(ax, view, (0,0,0)) 
    591         draw_mesh(ax, view, jitter, dist=dist, n=30) 
     653        draw_mesh(ax, view, jitter, dist=dist, n=mesh, projection=projection) 
    592654        draw_scattering(calculator, ax, view, jitter, dist=dist) 
    593655        plt.gcf().canvas.draw() 
     
    607669    plt.show() 
    608670 
     671def main(): 
     672    parser = argparse.ArgumentParser( 
     673        description="Display jitter", 
     674        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
     675        ) 
     676    parser.add_argument('-p', '--projection', choices=PROJECTIONS, default=PROJECTIONS[0], help='coordinate projection') 
     677    parser.add_argument('-s', '--size', type=str, default='10,40,100', help='a,b,c lengths') 
     678    parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, default=DISTRIBUTIONS[0], help='jitter distribution') 
     679    parser.add_argument('-m', '--mesh', type=int, default=30, help='#points in theta-phi mesh') 
     680    parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], help='oriented shape') 
     681    opts = parser.parse_args() 
     682    size = tuple(int(v) for v in opts.size.split(',')) 
     683    run(opts.shape, size=size, 
     684        mesh=opts.mesh, dist=opts.distribution, 
     685        projection=opts.projection) 
     686 
    609687if __name__ == "__main__": 
    610     model_name = sys.argv[1] if len(sys.argv) > 1 else 'parallelepiped' 
    611     size = tuple(int(v) for v in sys.argv[2].split(',')) if len(sys.argv) > 2 else (10, 40, 100) 
    612     main(model_name, size) 
     688    main() 
Note: See TracChangeset for help on using the changeset viewer.