Monday, February 15, 2021

Animated GIFs for Scientific Visualization - An Example Antenna Array Animation using Python, Matplotlib & GIMP

This animation shows how the radiation patterns from an array of antennas as more antennas are added. Using an animation its a lot easier to quickly see how the patterns changes with the number of antennas. Saving the results as a GIF file rather than an animation using Python allows the results to be used in a website or in a Powerpoint presentation. As of Matplotlib 3.3, there is not a gif export capability, but this king of animation can be generated using Python, Matplotlib, and GIMP. 


Discussion

Often scientific information and simulation results need to be visualization to be explained and understood. Live animations in a script are useful, but distributing scripts has a lot of issues. Creating movies (mp4, avi) sometime have issues with encoders. GIF animations are a reliable way to share a lot of information. They can be embedded in PowerPoint, send via email, or embedded in a website with minimal effort and the results just work. 

Once a script to generate a Matplotlib animation, it can usually be modified to generate a sequence of images and save them using the 'savefig' command. Then there are several tools that can convert the images to an animated GIF. With GIMP, the GIF can be generated in a a few minutes and there are a lot of options to control the speed, resolution, and resolution of the result.


Steps to Generate An Animated GIF

  1. Generate a sequence of images of identical size using Matplotlib and save using savefig.
  2. Open GIMP
    • File -> Load the plots as layers
      • Select all image files
    • Filters -> Animation -> Optimize (for GIF)
      • Image->Mode->Index
      • Generate optimum palette
      • Maximum number of colors: 256
      • Color dithering: Positioned
      • Enable dithering of transparency: checked
      • Select Convert
    • Filter -> Animation -> Optimize (GIF)
    • File -> Export
      • Name the file with a .gif extension
      • Select 'Export'
      • On new dialog
        • Select as animation
        • Set timing to 20 msec

    References

    Sample Script


    
    # -*- coding: utf-8 -*-
    
    import matplotlib.pyplot as plt
    from matplotlib.colors import BoundaryNorm
    from matplotlib.ticker import MaxNLocator
    import numpy as np
    
    
    
    # make these smaller to increase the resolution
    dx, dy = 0.005, 0.005
    
    # list of antenna locations
    wavelength = 0.1
    spacing = 0.25
    x_list = np.array([0,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,7,-7])*spacing*wavelength
    y_list = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
    phase_deg_list = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,] 
    phase_offset_deg_list = range(0,720,9)
    
    num_antennas = [1+int(pp*8/720)*2 for pp in phase_offset_deg_list]
    
    
    
    for idx,(phase_anim,antennas) in enumerate(zip(phase_offset_deg_list,num_antennas)):
        
        # generate 2 2d grids for the x & y bounds
        y, x = np.mgrid[slice(-1, 1 + dy, dy),
                        slice(-1, 1 + dx, dx)]
        
        z_list= []
        for xx,yy,phase_offset in zip(x_list[:antennas],y_list,phase_deg_list):
            phase_shift = 2*np.pi*(phase_offset + phase_anim)/360
            zz = np.real(np.exp(1j*phase_shift)*np.exp(np.pi*2/wavelength*-1j*np.sqrt((x-xx)**2 + (y-yy)**2)))
            z_list += [zz]
        
        
        z = sum(z_list)/len(z_list)
        
        # x and y are bounds, so z should be the value *inside* those bounds.
        # Therefore, remove the last value from the z array.
        z = z[:-1, :-1]
        levels = MaxNLocator(nbins=15).tick_values(z.min(), z.max())
        
        
        # pick the desired colormap, sensible levels, and define a normalization
        # instance which takes data values and translates those into levels.
        cmap = plt.get_cmap('bwr')
        norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)
        
        fig, (ax1) = plt.subplots(figsize=(12,10),dpi=75)
        
        
        # contours are *point* based plots, so convert our bound into point
        # centers
        ax1.plot(x_list,y_list,'o',linestyle='',markerfacecolor='black')
        ax1.plot(x_list[:antennas],y_list[:antennas],'ko',linestyle='',markerfacecolor='white')
        cf = ax1.contourf(x[:-1, :-1] + dx/2.,
                          y[:-1, :-1] + dy/2., z, levels=levels,
                          cmap=cmap,
                          vmin=-1, vmax=1.0)
        #plt.clim(-1,1)
        fig.colorbar(cf, ax=ax1)
        ax1.set_title('Number of Antennas = %i' % antennas)
        
        
        # adjust spacing between subplots so `ax1` title and `ax0` tick labels
        # don't overlap
        fig.tight_layout()
        plt.savefig('image-%03i.png' % idx)
        
        plt.show()
    
    

    Sunday, February 14, 2021

    Generating 3D Plots and Exporting the Geometry - An Example with Antenna Radiation in Matplotlib, PyVista and Blender

    This was made starting from a Python script using Matplotlib. This post shows how. 

    Discussion

    Matplotlib is a great tool for generating beautiful 2D and 3D images quickly. 


    However, sometimes its desirable to generate something more visual and appealing to help explain a concept more clearly. Often this is not possible in Matplotlib. Other tools must be used.  However, there is not facility in Matplotlib (as of 3.2.1) to export that plot surface for use in other tools. 
    Fortunately, using almost the same syntax as Matplotlib, a 3D plot can be created in PyVista and exported in one of the common 3D file formation like stl, ply, or fbx. Once this is done, its possible to import into a tool like Blender and create a high quality rendering or animation to help explain a concept. 

    This posting will cover the key steps in doing this and use the plotting of the radiation pattern of an antenna to help illustrate.


    Setup

    You'll need Python 3 with Matplotlib, PyVista, and Blender to replicate this. My personal setup uses the portable version of Python 3 and Blender.
    1. Install WinPython portable (https://winpython.github.io/)
    2. Use pypi to install PyVista (https://pypi.org/project/pyvista/)
    3. Install Blender portable (https://www.blender.org/download/ and chose the portable version)

    Example

    Sunday, March 18, 2012

    Installing sailfish-CFD under PythonXY

    image

    Discussion

    Sailfish CFD is an interesting python program. It uses OpenCL to solve fluid flow problems using the Lattice-Boltzmann method. Sailfish uses PyOpenCL or PyCuda to manage the simulations and pygame as one of its visualization packages. The following steps can be used to add Sailfish CFD to PythonXY and execute the examples using OpenCL.

    Installation Steps:

    1. Install the OpenCL driver for the computer from the CPU/GPU vendor. NVIDIA drivers can be found here. Intel Drivers can be found here.
    2. Download the source using Git following the links here. If you do not use Git or don’t want to add it to your machine, a portable version is available here.
    3. Since the developers do not offer packages (yet?), use git to clone their repository. I clone to a temp directory, then work from there.
      • git clone git://gitorious.org/sailfish/sailfish.git c:\temp\sailfish
      • The repository can be viewed through a GUI using the gitk command.
    4. I like to keep the python pieces together, so I copy the contents of the sailfish directory from my temp location to the a sailfish directory created under c:\Python27.
    5. There appears to be an issue with the location of the Mako files. So following the recommendation on the Sailfish google group in this exchange, copy all of the template files from “C:\Python27\sailfish\sailfish\templates” to “C:\Python27\sailfish\sailfish”.
    6. Create a file named “sailfish.pth” at “C:\Python27\Lib\site-packages”. Edit the file and add a single line with “C:\Python27\sailfish\”. This tells python to look in that directory for the sailfish modules.
    7. Test the installation by opening Python and typing “import sailfish” at the prompt. If it imports without errors, your installation is probably good.
    8. Since PythonXY works with PyOpenCL out of the box and Sailfish works with OpenCL, the examples can now be run by using the “—backend=opencl” option. Open up “C:\Python27\sailfish\examples\lbm_cylinder.py” and try to run using the opencl option.
    9. If your system does not have a GPU, you might encounter an error like “pyopencl.LogicError: Context failed: invalid value”. This is because Sailfish is configured to only use GPUs. I was able to get my installation to work on a system without a GPU, but with Intel OpenCL installed. This was done by modifying “C:\Python27\sailfish\sailfish\backend_opencl.py”  modifying line 31 from “devices = platform.get_devices(device_type=cl.device_type.GPU)” to “devices = platform.get_devices()”.
    The image at the top of the posting was generated from “C:\Python27\sailfish\examples\lbm_cylinder.py”.

    References:


    Test Environment:

    • Win 7 Professional
    • PythonXY 2.7.2.1
    • Sailfish git commit with SHA 1 of 29d7c934be8c02cf714ce2307796a4550428b512 from 2011-11-06 at 14:35:56
    • Intel OpenCL
    • i7 CPU, no GPU
    This work is licensed under a Creative Commons Attribution By license.

    Saturday, March 17, 2012

    Yet Another Countdown Timer Application built using .HTA (Windows HTml Applications)

    image

    Discussion:

    The time until something is due is always important. There are a lot of countdown timers which are available. There are a lot of interesting ways to implement a countdown timer. In restricted IT environments, where installing applications is frowned upon, one approach to creating a lightweight application is using HTML Applications (HTA).  The following example shows how to use an HTA to build a standalone countdown timer which is displayed as a window.

    Several techniques are used to make this countdown timer work. First, the name of the file is parsed to determine the target date and time. This way, the remaining time can be set by renaming the file. Next, javascript is used to rewrite a portion of the document by assigning a value to the ‘innerHTML’  of  div element. Then, the countdown function is called once when the body of the document is loaded. However, once the function is called, the last thing it does is reschedule itself to run again in 1 second. Finally, this file, which looks like an HTML file is saved as with the extension HTA.

    Code Example:

    <html>
    <head>
    	<TITLE>Countdown Timer</TITLE>
    	<HTA:APPLICATION 
    		ID="oMyApp" 
    		APPLICATIONNAME="CountDownTimer" 
    		BORDER="yes"
    		CAPTION="yes"
    		SHOWINTASKBAR="yes"
    		SINGLEINSTANCE="yes"
    		SYSMENU="yes"
    		WINDOWSTATE="normal"
    		MAXIMIZEBUTTON="no">
        <script>
    		// set the window size and place it on the screen in the upper right hand corner
    		var windowWidth = 600
    		var windowHeight = 150
    		window.resizeTo(600,150);  // Can only be done in HTA
    		//half the screen height minus half the new window height (plus title and status bars).
    		iMyHeight = (window.screen.height/4) - (windowHeight/2 + 50)/2;
            window.moveTo(window.screen.width-windowWidth-5,5);
        </script>
    	<style type="text/css">
    		body {
    			background-color: green;
    			color: black; 
    			text-align: center; 
    			font-family: arial; 
    			font-weight: bold; 
    			font-size: 20px;
    			font-variant: small-caps;
    			border: medium double rgb(0,255,15);
    			width:100%;
    			overflow: hidden;
    			filter:progid:DXImageTransform.Microsoft.Gradient(endColorstr='#00C000',
    						startColorstr='#00FFFF', gradientType='0'); 			
    			}
    	</style>
    </head>
    <body>
    	<!-- Define a division for displaying the countdown message -->
    	<div id="countdown" ></div>
    	<!-- Define the actions to take when the body loads -->
    	<script>
    	// parse the filename to get the target date and time
    	var filename = String(window.location)              // get filename
    	var iSubStrEnd = filename.lastIndexOf('.')       // get point in string where the date data may end
    	var iSubStrStart = filename.lastIndexOf('.',iSubStrEnd-1) // get point in string where the date data may start
    	if ((iSubStrEnd>0)&&(iSubStrStart>0)&&(iSubStrStart<iSubStrEnd)) {
    		// get the date string from the file name
    		endDateTime = filename.slice(iSubStrStart+1,iSubStrEnd)
    		}
    	else {
    		// enter the stop date here that will be used if the name is not formatted correctly
    		// format for this string is YYYY,MM,DD,HH,MM,SS
    		endDateTime = '2015,01,01,0,0,0'
    		}
    		
    	// parse the date
    	var params = endDateTime.split(',')
    	lyear  = parseInt(params[0])
    	lmonth = parseInt(params[1])
    	lday   = parseInt(params[2])
    	// call the countdown function which will hook itself into a timer and
    	//     continue to run every second.
    	countdown(lyear,lmonth,lday,params[3],params[4],params[5])	
    
    	//-------------------------------------------------------
    	function countdown(yr,m,d,hh,mm,ss){
    		var montharray=new Array("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")
    		hh = ("00"+hh).slice(-2);
    		mm = ("00"+mm).slice(-2); //if (mm<10) mm = "0"+(mm+1);
    		ss = ("00"+ss).slice(-2); //if (ss<10) ss = "0"+ss;
    		theyear=yr; themonth=m; theday=d;
    		thehh=parseInt(hh); themm=parseInt(mm); thess=parseInt(ss);
    		var today=new Date();
    		var todayy=today.getFullYear();
    		var todaym=today.getMonth()+1;
    		var todayd=today.getDate();
    		var todayh=today.getHours();
    		var todaymin=today.getMinutes();
    		var todaysec=today.getSeconds();
    		var todaystring=montharray[todaym]+" "+todayd+", "+todayy+" "+todayh+":"+todaymin+":"+todaysec
    		futurestring   =montharray[m-1]   +" "+d     +", "+yr    +" "+hh    +":"+mm      +":"+ss
    		var dd    = Date.parse(futurestring)-Date.parse(todaystring);
    		var dday  = Math.floor(dd/(60*60*1000*24)*1);
    		var dhour = Math.floor((dd%(60*60*1000*24))/(60*60*1000)*1);
    		var dmin  = Math.floor(((dd%(60*60*1000*24))%(60*60*1000))/(60*1000)*1);
    		var dsec  = Math.floor((((dd%(60*60*1000*24))%(60*60*1000))%(60*1000))/1000*1);
    		if(dday<=0&&dhour<=0&&dmin<=0&&dsec<=0){
    			// put text in the "countdown" division to show the countdown expired
    			document.getElementById("countdown").innerHTML = "<br> Time has run out! <br>";
    			return;
    			}
    		else {
    			dhour = ("00"+dhour).slice(-2);
    			dmin  = ("00"+dmin).slice(-2);
    			dsec  = ("00"+dsec).slice(-2);
    			// format the delta time and put in the countdown division to display
    			document.getElementById("countdown").innerHTML = 
    				"<br> Times out in "+dday+ 
    				" days, "+dhour+" hours, "+dmin+" minutes, and "+dsec+
    				" seconds  <br>";
    			}
    		// trigger this routine to run again in 1 second
    		setTimeout("countdown(theyear,themonth,theday,thehh,themm,thess)",1000);
    		} 
    
    	</script>
    	</body>
    </html>

     

    TEsting Environment:

    • IE9
    • Win7
    This work is licensed under a Creative Commons Attribution By license.

    Friday, February 24, 2012

    Lesson Learned: Managing Large Numbers of Plots, with an Example in Python

    28056200-5ef7-11e1-8741-08002700a056

    Problem Statement:

    A project generates hundreds, thousands, or more graphs over its life. These graphs are copied and pasted into e-mail, power point slides, etc. The plots become divorced from any of the documents they were originally distributed with. Invariably, at some point in the project, a plot is brought back and the question is what were the assumptions used to generate this graph. With only the graph available, it can be difficult or impossible to answer this question.

    To complicate matters, the plots are generated using legacy codes and modifying all of the existing code base is a detailed endeavor.

    How can this situation be improved?

    Discussion:

    There are two problems here. First, a given graph is not traceable to its origin. This can be remedied in one of many ways.  If the source data is well controlled and can be described using a short phase, then adding that phrase somewhere on the chart is helpful. If the source data is constantly changing or requires too much information to describe with a short phrase, then something else is needed. A hash of the input data can help identify and verify the source data set used to generate the graph. A universally unique ID (UUID) can be used to give a graph a unique name. If the source data, assumptions, etc are stored using that same UUID, then when a graph if brought back for review, all of the necessary parts can be found.

    The second problem is handling legacy code. There are at least three choices.

    • The first is probably the easiest. A function to add hash and uuid could be created and inserted at the appropriate location in each of the major pieces of plotting code. This is problematic because there are several interfaces and actions in the scripts which could make this work poorly. Also, every plotting routine would need to be modified.
    • A second choice is to wrap the plotting routines into function, then pass this function and its data to a wrapper function which would add a hash and uuid as the last thing done by the plotting routine. If the plotting functions already exist, then this can be done without changing any of the plotting code.   
    • A third choice is to create a decorator which wraps plotting routines, adding the hash and uuid. This has the same issues as using a wrapper call and requires changes to the plotting routines source. However, the changes consist of an import statement and application of the decorator at the correct location.

    For this problem, the decorator solution is used. The code that following implements a decorator that creates a hash of the plot data and a UUID which are added to the right side of the plot.  This way, no matter where the plot goes, there is a high likelihood that its pedigree can be preserved.

    '''
    Script to demonstrate the use of decorators to add a unique identifier to
       a plot. The identifier incudes a hash of input data, to help see if 
       version of a plot really have different data, and a UUID to uniquely 
       identify this plot independent of when or where it was generated.
    '''
    
    __author__  = 'Ed Tate'
    __email__   = 'edtate<AT>gmail-dot-com'
    __website__ = 'exnumerus.blogspot.com'
    __license__ = 'Creative Commons Attribute By - http://creativecommons.org/licenses/by/3.0/us/'
    
    from matplotlib.pylab import *
    import random
    import md5
    import uuid
    
    def identifiable_plot(fn):
        def newfun(*args,**kwargs):
            # do before decorated function call
            fn(*args,**kwargs)
            # do after decorated function call
            # create the tag string from a hash of the data and a 
            #    universially unique ID
            x = args[0]
            y = args[1]
            xy = zip(x,y)
            m = md5.new()
            m.update(str(xy))
            this_uuid = str(uuid.uuid1())
            this_tag = 'hash=' + m.hexdigest() + ',' + 'UUID=' + this_uuid
            # write the tag to the figure for future reference
            figtext(1,0.5,this_tag ,rotation='vertical',
                    horizontalalignment='right',
                    verticalalignment='center',
                    size = 'x-small',
                    )
            return this_uuid
    
        return newfun
    
    ###############################
        
    @identifiable_plot
    def my_plot(x,y):
        plot(x,y,'o')
        grid(True)
        
    ###############################
        
    x = [random.random() for i in range(100)]
    y = [random.random() for i in range(100)]
        
    plot_uuid = my_plot(x,y)
    savefig(plot_uuid+'.png')
    
    show()
    

     


    Test Configuration:
    • win7
    • PythonXY 2.7.2.1

    References:
    This work is licensed under a Creative Commons Attribution By license.