Mesh making material

README.md

-# OVERVIEW
\ No newline at end of file
+# 3Dskymaps_utils
+Hosting basic tools to create 3D skymaps

cron

+#!/bin/bash
+set -e
+SCRIPTDIR=$(dirname `realpath "$0"`)
+echo $SCRIPTDIR
+LOGDIR=/dev/stdout
+LOCKF=/var/lock/3dmap.lock
+
+# Read config
+. $SCRIPTDIR/settings
+
+flock -n $LOCKF $SCRIPTDIR/wrapper &>> $LOGDIR
+

gw_mesh.py

+from warnings import warn
+import healpy as hp
+import numpy as np
+from scipy.stats import norm
+from astropy.utils.data import download_file
+import ligo.skymap.io.fits as iofits
+
+import argparse as op
+import os
+import sys
+import trimesh
+import trimesh.transformations as tt
+import trimesh.creation as tc
+import pickle
+import mcubes
+from skimage import measure
+import yaml
+# @profile (use with python3 -m memory_profiler gw_mesh.py
+
+def dp_dV_cart(x, y, z, event):
+    npix = len(event['prob'])
+    nside = hp.npix2nside(npix)
+    pixarea = hp.nside2pixarea(nside)
+    ipix = hp.vec2pix(nside, x, y, z)
+    r = np.sqrt(x**2 + y**2 + z**2)
+    return (event['prob'][ipix] * event['distnorm'][ipix] *
+            norm(event['distmu'][ipix], event['distsigma'][ipix]).pdf(r) /
+            pixarea)
+
+def logdp_dV_cart(x, y, z, event):
+    npix = len(event['prob'])
+    nside = hp.npix2nside(npix)
+    pixarea = hp.nside2pixarea(nside)
+    ipix = hp.vec2pix(nside, x, y, z)
+    r = np.sqrt(x**2 + y**2 + z**2)
+    return (np.log(event['prob'][ipix]) +
+            np.log(event['distnorm'][ipix]) -
+            (r - event['distmu'] / event['distsigma'])**2 -
+            np.log(pixarea))
+
+def FindHeightForLevel(inArr, adLevels):
+    # Flatten the array assuming it is 3D, containing the log-density
+    oldshape = np.shape(inArr)
+    adInput = np.reshape(inArr, oldshape[0] * oldshape[1] * oldshape[2])
+    # Get array specifics
+    nLength = np.size(adInput)
+    # Reverse sort the flattened adInput
+    adInput[::-1].sort()
+    # Create normalized log - cdf (using log-density)
+    # Assumes a linearly spaced grid (ignored being just a constant)
+    adCum = np.log(np.cumsum(np.exp(adInput)))
+    adCum -= adCum[-1]
+    # Find value closest to required density levels
+    adHeights = []
+    for item in adLevels:
+        idx = (np.abs(adCum - np.log(item))).argmin()
+        adHeights.append(adInput[idx])
+    adHeights = np.array(adHeights)
+    return adHeights
+
+
+# def FindExtent(event, sigma):
+#     distlim = event['distmean'] + sigma * event['diststd']
+#     return np.max(distlim[np.isfinite(distlim)])
+
+
+if __name__ == '__main__':
+    
+    parser = op.ArgumentParser()
+    parser.add_argument("-S", "--superevent", type=str, dest="superevent",
+                        help="Superevent identifier",
+                        default=None,
+                        required=False)
+    parser.add_argument("-b", "--box", type=int, dest="box",
+                        help="Halfsize of the box (in Mpc, centered to the SSB)",
+                        default=600,
+                        required=False)
+    parser.add_argument("--adaptive", action='store_true', dest="adaptive",
+                        help="Find the rendering box size adaptively",
+                        default=False)
+    parser.add_argument("-s", "--sigma", type=int, dest="sigma",
+                        help="""Set the number of sigmas for the distance probability (only used for the adaptive strategy)""",
+                        default=3)
+    parser.add_argument("-n", "--npoints", type=int, dest="npoints",
+                        help="Number of point per box size",
+                        default=400)
+    parser.add_argument("-l", "--levels", type=float, nargs='+', dest="levels",
+                        help="""Probability contour level (one mesh file for each value)""",
+                        default=[0.90])
+    parser.add_argument("--lal", action='store_true', dest="lal",
+                        help="""Use lalinference fits instead of bayestar""",
+                        default=False)
+    parser.add_argument("--local", type=str, dest="local", help="Use local fits file", default=None)
+    parser.add_argument("--url", type=str, dest="url", help="""Use direct url to download file""", default=None)
+    parser.add_argument("--name", type=str, dest="name", help="""Set a name for events outside the superevent catalogue""", default=None)
+    parser.add_argument("--lewiner", action='store_true', dest="lewiner",
+                        help="""Use lewiner algorithm from skimage.measure""",
+                        default=False)
+    parser.add_argument("--smooth", action='store_true', dest="smooth",
+                        help="""Smooth and export the mesh into a .gltf file.""",
+                        default=False)
+    parser.add_argument("--annulus", action='store_true', dest="annulus",
+                        help="""Add an horizontal annulus to identify the galactic plane""",
+                        default=False)
+    parser.add_argument("-c", "--config", dest="config",
+                        help="Path to YAML configuration file.\n"
+                        +"Can be used to override any arguments provided by the command line interface.",
+                        default=None)
+    parser.add_argument("-o", "--outputdir", dest="outputdir",
+                        help="Folder in which to output files.\n"+
+                        "Default behaviour is to output file into a folder with the event name.",
+                        default='')
+    options = vars(parser.parse_args())
+     
+    if options['config'] is not None:
+        with open(options['config']) as stream:
+            for key,value in yaml.safe_load(stream).items():
+                if key in options:
+                    options[key] = value
+    
+    # Load the input data (by default is bayestar, assuming knowledge of superevent
+    remote_file = "https://gracedb.ligo.org/api/superevents/{}/files/bayestar.fits.gz".format(options['superevent'])
+    
+    # Load the LALinference instead (still assuming its the default one)
+    if options['lal'] is not None:
+        remote_file= "https://gracedb.ligo.org/api/superevents/{}/files/LALInference.fits.gz".format(options['superevent'])
+
+    # Use a direct url (useful for testing or for manual jobs)
+    if options['url'] is not None:
+        remote_file = options['url']
+
+    # Use for locally stored file
+    if options['local'] is not None:
+        file = options['local']
+    else:
+        file = download_file(remote_file)
+
+    prob, distmu, distsigma, distnorm = hp.read_map(file, field=[0, 1, 2, 3])
+
+    # Assuming a conventional structure for the fits file
+    event = {'prob': prob, 'distmu': distmu,
+             'distsigma': distsigma, 'distnorm': distnorm}
+    
+    # Set the grid points
+    n = options['npoints']
+    
+    # Choose the bounding box, either adaptively or fixed by input
+    if options['adaptive']:
+        mapfits = iofits.read_sky_map(file)
+        extent = mapfits[1]['distmean'] + options['sigma'] * mapfits[1]['diststd']
+    else:
+        extent = options['box']
+
+    ext = [-extent, extent,
+           -extent, extent,
+           -extent, extent]
+
+    # Setting the grid (MEMORY DEMANDING)
+    x, y, z = np.mgrid[ext[0]:ext[1]:n * 1j,
+                       ext[2]:ext[3]:n * 1j,
+                       ext[4]:ext[5]:n * 1j]
+    
+    scalars = dp_dV_cart(x, y, z, event) # Could be parallelized
+    heights = FindHeightForLevel(np.log(scalars), options['levels']) # Could be parallelized
+   
+    # create output directory
+    if options['superevent'] is not None:
+        name = str(options['superevent']).upper()
+    else:
+        name = str(options['name'])
+    
+    if options['outputdir'] == '':
+        folder = name
+    else:
+        folder = options['outputdir']
+
+    try:
+        os.mkdir(folder)
+    except FileExistsError:
+        warn('Folder {} already exists.'.format(folder))
+
+    # Creating the mesh
+    for l,h in zip(options['levels'], heights):
+        if options['lewiner']:
+            sys.stdout.write("Using Lewiner algorithm...\n")
+            vertices, triangles, _ , __= measure.marching_cubes_lewiner(np.log(scalars), level=h, allow_degenerate=False)
+            # Add here cube for rendering
+            sys.stdout.write("Finished\n")
+        else:
+            vertices, triangles = mcubes.marching_cubes(np.log(scalars),h)
+        vertices -= int(n/2)
+        vertices *= (x[1][0][0]-x[0][0][0])/1000 # Output in Gpc
+        sys.stdout.write("Most distant point on the contour surface: {}\n".format(vertices[np.argmax(np.linalg.norm(vertices, axis=1))]))
+
+        filename = os.path.join(folder,"{}.dae".format(l))
+        sys.stdout.write("Saving {}...\n".format(filename))
+        mcubes.export_mesh(vertices, triangles, filename, name)
+
+        if options['smooth']:
+            sys.stdout.write("Smoothing {}...\n".format(filename))
+            with open(filename, 'rb') as f:
+    	        scene = trimesh.load(file_obj=f, file_type='dae')
+            mesh = scene.geometry['geometry0.0'] #Do not change this!
+            lap = trimesh.smoothing.laplacian_calculation(mesh, equal_weight=True)
+            __ = trimesh.smoothing.filter_laplacian(mesh, 0.5, 10, False, True, lap)
+
+            # Add annulus as reference system. Diameter equal
+            if options['annulus']:
+                box = extent
+                rmin=box/1000
+                rmax=rmin*1.01
+                height = rmax-rmin
+                sections = 256
+                transform = np.array([[1,0,0,0],[0,0,-1,0],[0,1,0,0],[0,0,0,1]])
+                my_annulus = tc.annulus(rmin, rmax, height, sections, transform)
+                scene.add_geometry(geometry=my_annulus, geom_name='annulus')
+
+            filename_smooth = os.path.splitext(filename)[0]+'.glb'
+            sys.stdout.write("Saving smooth {}...\n".format(filename_smooth))
+            scene.export(filename_smooth, file_type='glb')

require.txt

+appdirs==1.4.3
+astropy==3.2.3
+backcall==0.1.0
+certifi==2019.11.28
+cycler==0.10.0
+decorator==4.4.1
+distlib==0.3.0
+filelock==3.0.12
+healpy==1.13.0
+imageio==2.8.0
+importlib-metadata==1.6.0
+importlib-resources==1.4.0
+ipython==7.9.0
+ipython-genutils==0.2.0
+jedi==0.16.0
+kiwisolver==1.1.0
+matplotlib==3.0.3
+meshtool==0.3
+networkx==2.4
+numpy==1.18.1
+parso==0.6.2
+pexpect==4.8.0
+pickleshare==0.7.5
+Pillow==7.0.0
+pipenv==2018.11.26
+pkg-resources==0.0.0
+prompt-toolkit==2.0.10
+ptyprocess==0.6.0
+pycollada==0.7.1
+Pygments==2.6.1
+PyMCubes==0.0.12
+pyparsing==2.4.6
+python-dateutil==2.8.1
+PyWavelets==1.1.1
+scikit-image==0.15.0
+scipy==1.4.1
+six==1.14.0
+traitlets==4.3.3
+trimesh==3.6.17
+virtualenv==20.0.15
+virtualenv-clone==0.5.4
+wcwidth==0.1.9
+zipp==1.2.0

wrapper

+#!/bin/bash
+echo "Unless you know what you are doing"
+echo "This script should be run with"
+echo "flock -n /var/lock/3dmap.lock wrapper"
+
+set -e
+SCRIPTDIR=$(dirname `realpath "$0"`)
+echo $SCRIPTDIR
+
+# Arguments (can be overriden in config file)
+ulimit -m 32768000
+ulimit -v 32768000
+ulimit -H -v
+export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
+export PATH=/usr/local/bin:$PATH
+VENV=~/3d/bin/activate
+OUTPUTDIR=~/3DOUT
+
+# Read config
+. $SCRIPTDIR/settings
+echo "User: $USER"
+echo "Server: $SERVER"
+echo "Skymapdir: $SKYMAPDIR"
+echo "Queuedir: $QUEUEDIR"
+
+# canonicalize output path
+OUTPUTDIR=$(realpath $OUTPUTDIR)
+VENV=$(realpath $VENV)
+mkdir -p $OUTPUTDIR
+cd $OUTPUTDIR
+
+# Get list of outstanding jobs
+file=$(echo "ls" | sftp -q $USER@$SERVER:$QUEUEDIR )
+echo $file
+
+# Get first job
+file=$(echo "$file" | grep -v '^sftp>' | head -n1 | awk '{print $1;}')
+if [$file == '']; then
+    echo "No new files"
+    exit
+fi
+echo "Getting $file"
+echo "GET $file" | sftp -q $USER@$SERVER:$QUEUEDIR
+echo "rm $file" | sftp -q $USER@$SERVER:$QUEUEDIR
+
+# Overrite outputdir
+echo "Patching security"
+sed -i '/outputdir/d' $file
+echo "outputdir: $OUTPUTDIR/3D" >> $file
+
+# Delete existing skymaps (if they exist)
+echo "Removing old skymap files"
+echo "rmdir 3D" | sftp -q $USER@$SERVER:$SKYMAPDIR/$file/ || true
+
+# Make new skymaps
+echo "Making new files"
+source $VENV
+python $SCRIPTDIR/gw_mesh.py --config $file
+deactivate
+
+# Copy to server and clean up
+echo "Cleaning up"
+echo "PUT -r $OUTPUTDIR/3D" | sftp -q $USER@$SERVER:$SKYMAPDIR/$file/
+rm -r $OUTPUTDIR/*
+
+echo "Done"
+