#!/usr/bin/env python # ----------------------------------------------------------------------------------- # # Python/ROOT script # # Author: Florian Uekermann (NBI) # Email: # Date: 16th of September 2014 # # Additions: Troels Petersen (NBI) # Email: petersen@nbi.dk # Date: 6th of January 2016 # # ----------------------------------------------------------------------------------- # # # This exercise is a little different in style, as it mostly involves commenting in # calls to external functions, which are contained in the file "timeSeries_functions.py". # Occationally, a little coding is needed, also in the external function file, but # otherwise the exercise should be straight forward, and illustrate well a complex analysis. # # ----------------------------------------------------------------------------------- # import ROOT import array import math import time execfile('timeSeries_functions.py') # Load function definitions # You will analyse pictures, which I obtained by filming a small # slime mold plasmodium over a few hours. The pictures you will see # have already been processed in various ways and only contain pixels # with the values 0 and 1, where 1 means that in the original microscope # picture was part of the plamodium. # Slime molds contract in a cyclic manner to move stuff around inside # themselves (a slime mold consists of one huge cell). We want # to measure the typical duration of a contraction-expansion cycle. # TASK 1: Importing the data & understand the data format. # Use the ImportJSON function to store the content of the data file in a variable. # The filename is 'images.json.gz' and it contains a list of images. # Each image consists of a list of list of int (a 2D list because we have # a 2D coordinate system). Note that you can specify how many pictures are # imported (optional) in the ImportJSON function, to save time (or prevent # out of memory crashes). #YOUR CODE HERE import_result = ImportJSON("images.json.gz", 400) # Can be extended up to 2400 frames! images = import_result[0] times = import_result[1] # Draw the third image in the images list, using the DrawImage function #YOUR CODE HERE image3 = images[2] #for i in range ( 10, 45 ) : # for j in range ( 10, 45 ) : # image3[i][j] = 1 DrawImage(image3, False) # To check if you understood the data format (this will be important later), # change the value of the pixel with the coordinates (50,100) to 1 in the # picture you have drawn above and draw it again. #YOUR CODE HERE # See above! # Have a look at the whole dataset by calling the AnimateImages function # with the whole list of pictures (You can change the speed if you want) #YOUR CODE HERE #AnimateImages(images, 0.02) # TASK 2: Compute cell cross-section area. Plot an area time series. # Fill in the missing code for the SumOfPixelValues function. Test it by # printing the sum of values in picture 3. # Print the total number of pixels (not values) of picture 3. Check if # the fraction of the two numbers make sense (compare with area in drawings). # Hint: Use the len function for finding the dimensions of the picture. # Feel free to get rid of the animation in the code above (to save time). #YOUR CODE HERE npix3 = SumOfPixelValues(image3) # print npix3 # We can use the SumOfPixels function to obtain an estimate of the # Area of the cross-section of the cell in each picture. # Compute this area for each picture and store the results in # a list. #YOUR CODE HERE npix = [] for i in range (len(images)) : npix.append(SumOfPixelValues(images[i])) # print npix # TASK 3: Separate slow cell area changes (cell-growth) and fast # ones (contractions, expansions) using a moving average and FFT. # Fill in the missing code for the MovingAverage function. Test it with a window # of 22 data points in each direction on the time series of the area (45 in total). # Note that the resulting time series should have the same length as the orignal # one and should not display any artifacts at beginning and end. #YOUR CODE HERE canvas = DrawGraph(times, npix, "", "time", "area") # Use the DFT function to calculate the fourier coefficients # of the area time series. Use the Amplitude function to # transform the amplitudes into coefficients. # You can use the Frequencies function to generate a # list of frequencies that matches the coefficients # and amplitudes. # Plot the amplitudes over frequencies (you can use # the Log function for rescaling the amplitudes). #YOUR CODE HERE print " Number of time steps in total = %d"%(len(times)) print " Time step (times[1]) = %6.3f minutes"%(times[1]) fourier_result = DFT(npix) amp_result = Amplitudes(fourier_result) log_amp_result = Log(amp_result) freq_result = Frequencies(len(amp_result), times[1]) canvas2 = DrawGraph(freq_result, log_amp_result, "", "frequency", "amplitude") # Compare the log-amplitude-frequency plot with the # area time series. Which frequency range corresponds to # the contractions and which part to cell growth? Find # an appropriate threshold (it may help make a plot of # a small part of the time series for easier frequency estimation). # Use the LowPassFilter function to set all coefficients # with lower frequency than the threshold to zero. # Then transform the filtered coefficients back into a # time series using the InverseDFT. Plot the resulting time # series together with the original time series of the area. #YOUR CODE HERE fourier_filtered = LowPassFilter(freq_result, fourier_result, 0.5) npix_lowpassfiltered = InverseDFT(fourier_filtered) DrawGraph(times, npix, "", "time", "area") AddGraph(times, npix_lowpassfiltered) # If you carefully tune your threshold, the low-pass # filtered area curve should reflect the slow trends # better than the moving average. The moving average # always shows a small bump for each contraction cycle, # which is undesirable. The low pass filtered curve # also reflects quick changes in area better that are # only slightly above the timescale of individual # contraction cycles. Compare the moving average plot # and the low pass plot to find these features. Adjust # your threshold if necessary. # Can you identify disadvantages of the low pass filter # method, compared to the moving average? # TASK 4: Use an autocorrelation to fit the residuals # of the smoothened time series. # Calculate the difference of the original # time series and one of the two # smoothened time series (low pass or moving average). # Use the Residuals function. Plot the result over # time. #YOUR CODE HERE npix_highpassfiltered = Residuals(npix, npix_lowpassfiltered) DrawGraph(times, npix_highpassfiltered, "", "time", "filtered area") # Fill in the missing two lines in the autocorrelation function. # Calculate and plot the autocorrelation of the residuals. # Hint: Use the same list of times you use for your time # series plots as x-values in the autocorrelation plot. # 2nd Hint: Limit the autocorrelation plot to the first # 10min (200 data points). Then you can easily see the # length of a single contraction cycle. #YOUR CODE HERE autocorr = Autocorrelation(npix) autocorr_highpassfiltered = Autocorrelation(npix_highpassfiltered) DrawGraph(times, autocorr_highpassfiltered, "", "time", "AutoCorrelation") AddGraph(times, autocorr) print autocorr, autocorr_highpassfiltered # My lab report states an estimated contraction period # of ~62.5 seconds. Does that fit with your own estimate? raw_input("Press any key to exit")