2

I have been bouncing some ideas around about this problem, but thought I would consult with the online community to see if someone has a better option.

So I have step-like function graphs like this:

step-1 graph

step-2 graph

step-3 graph

And with this I want to compute the y-displacement between steps.

As one can see the steps are not perfectly horizontal, but rather take a small range of y-values before stepping up.

So the question is:

(1) What is the 'proper' way (if there is one) to take the average y-value of each 'level'? I am not sure where I should be using as my left most point and my right most point on each level- so I can take the values between these points and average them, and then attain the 'average' of each level, hope this makes sense. As one can see they do not all span the same displacement in x. The ultimate goal is to get the y-displacement between levels, once I have the 'average' value of each 'level' it is trivial to take the difference.

I was maybe thinking of taking the derivative of the curve and seeing where it is equal to zero for my left and right most points on each level but i'm not sure it will work as each level also contains points where (dy/dx = 0) -- so I could use some insight.

Thank you :) Oh - and this must be done in Python -- and it is not just these graphs, but a lot of them with a similar style, so the code must be generic enough to handle other step-like graphs.

Date file for graph 1: http://textuploader.com/5nwsh

Data file for graph 2: http://textuploader.com/5nwsv

Data file for graph 3: http://textuploader.com/5nwsj

Scatter Plot Python code: 

import numpy as np
import matplotlib.pyplot as plt
import pylab as pl


data=np.loadtxt('data-file')
x= data[:,0]
y=data[:,1]


pl.plot(data[:,0],data[:,1],'r')

pl.xlabel('x')
pl.ylabel('y')
plt.show()
Improve this question
8
  • @roadrunner66 Would it matter if we tilted or not? The difference between levels will be the same y-displacement? Commented Mar 30, 2016 at 23:34
  • I believe your solution should lie in a model for the physical reality you are trying to understand. E.g. would you want to rotate the graphs first to align with x and y or is the current coordinate system the one you need? If the steps are very periodic, an FFT might be a way to test periodicity. The derivative would allow you to remove slant by putting an offset in such that most points would average zero as you mentioned. Commented Mar 30, 2016 at 23:42
  • No. If you tilted the graphs to be approximately horizontal, there would be no step in the y-axis. So step needs to be defined, I think you mean it in the sense where the graphs are tilted in a way to resemble stairs, like in the 2nd graph. Commented Mar 30, 2016 at 23:45
  • Yes I see what you mean. In that case, yes the coordinate system I am using right now is very important. The y-displacement on these 'slanted' graphs between step levels is what I am after. Commented Mar 30, 2016 at 23:49
  • I like the derivative, it gives you a function that (for the last data set) that looks like an inverted comb-function (with a positive offset for the slant). Now the data in the last set are unusually regular. To model this comb, you could make a sequence of pulses with irregular spacing and irregular width OR height (the integral over the pulse is your step height, so only one variable is needed for the 'pulse'), then fit the whole thing to the derivative of your data. Commented Mar 31, 2016 at 0:03

2 Answers 2

Reset to default
2

Let's assume your data are AFM measurements of a crystal surface (all units in m) and you want to get the step height of the crystal. The following would get you towards that. 

from __future__ import division
from ipywidgets import *
import numpy as np
import matplotlib.pyplot as p
%matplotlib inline

def rotatedata(x,y,a):
    cosa=np.cos(a)
    sina=np.sin(a)
    x = x*cosa-y*sina
    y = x*sina+y*cosa
    return x,y

data=np.loadtxt('plot3.txt')
data=data.T
x,y=data[0],data[1]

def workit(a2):
    fig=p.figure(num=None, figsize=(18, 16), dpi= 80, facecolor='w', edgecolor='k')
    p.subplot(511) # , aspect='equal')
    p.plot(x,y)

    #what is the slope? 
    m,b = np.polyfit(x, y, 1)

    x1,y1=rotatedata(x,y, -np.arctan(m) )  # rotate data to have flat surface, 
                                           # interesting for surface roughness
    p.subplot(512 )
    p.plot(x1,y1)

    x2,y2=rotatedata(x,y, a2)    # rotate data to 'sharpen' histogram
    p.subplot(513 )
    p.plot(x2,y2)

    p.subplot(514)
    p.hist(y2,bins=130)        

    y3=np.diff(y2)
    p.subplot(515 )
    p.plot(y3)

    return HTML()

interact(workit,a2=[-0.002,0.002,0.00001])

enter image description here

The first plot is your raw data. In the 2nd plot I removed the slope of the data to show data that you would use if you cared about calculating the surface roughness.

The third plot shows the same data if they are rotated such that all slopes are horizontal.

How this is done (interactive slider) is shown in the 4th plot, which is a histogram of the rotated data. You simple move the slider (rotate the data) until the histogram has maximum sharpness (all maxima are of minimum width). I did this by hand (slider) not with a function but a simple autofocus routine, maximizing the sum of the absolute differences between neighboring values would be OK to use.

In the last plot I show the 1st derivative (redundant now, but average slope of the horizontal regions should center around zero as a double-check).

The width of the crystal layers (step height you were asking for) is now given by the distances of the maxima in the histogram.

I'll leave the actual determination of each step as an exercise (e.g. threshold the histogram, then calculate the center of gravity for each peak, then get differences of neighboring peak centers of gravity).

Improve this answer
Sign up to request clarification or add additional context in comments.

7 Comments

Ya, that looks pretty good - thank you :) I have been trying to implement the autofocus for the histogram, but im kinda stuck - any suggestions?
also, how can I threshold the histogram?
Thresholding means you remove all values below a certain level, which you could do with a loop, or more elegantly with a numpy construct like : aa[aa<3.5e-9] =0. I'll take a look at the autofocus later, but try something yourself, like make a sum of the absolute values squared of the differences of all consecutive values in the histogram. This metric should peak, when the histogram is very sharp.
I have been trying to do the autofocus for all day, as the slider seems finicky - I could use your help for sure. thank you kindly :)
I don't think I want to get rid of values under a certain value as you mention. I want to get rid of the bins with a low count rate - or?
|
1

The 2nd step here with automatic leveling of the histogram. I define a function focus that tells me how sharp the histogram is then I use that function to find the angle that gives me the sharpest histogram. I then pick that angle, threshold the histogram and find the center of gravity of each histogram peak. The differences between those locations are the steps. If the numbers were real in meters the steps would be about 4 Angstrom.

It turns out that the same idea can be used in 2D to level AFM, or white light interferometer data and very precisely determine step heights not just of crystals, but of nanoscale coatings or etch depths. 

from __future__ import division
from ipywidgets import *
import numpy as np
import matplotlib.pyplot as p
%matplotlib inline

def rotatedata(x,y,a):
    cosa=np.cos(a)
    sina=np.sin(a)
    x = x*cosa-y*sina
    y = x*sina+y*cosa
    return x,y


data=np.loadtxt('plot3.txt')
data=data.T
x,y=data[0],data[1]

def rotateAndCheck(a2):
    x2,y2=rotatedata(x,y, a2) 
    vals,edges=np.histogram(y2,bins=230)
    focus=np.sqrt(np.sum((np.diff(vals))**2))
    return focus

focus=[]
amin,amax,astep=-0.01,0.01,0.0001
for i in np.arange(amin,amax,astep):
    focus.append(rotateAndCheck(i))


fig=p.figure(num=None, figsize=(18, 16), dpi= 80, facecolor='w', edgecolor='k')


p.subplot(311)    

p.plot(focus,'.-')
nm=np.argmax(focus)
angle=amin+astep*nm


p.subplot(312)
x2,y2=rotatedata(x,y, angle) 
vals,edges,_=p.hist(y2,bins=230)


#now threshold
p.subplot(313)
vals[vals<3]=0
#print len(edges),len(vals)
deltaedge=edges[1]-edges[0]
#print deltaedge
#p.bar(edges[:-1],vals,0.05e-10)
p.bar(np.arange(len(vals)),vals,0.05e-10)
p.show()

# now you go through the histogram from left to right, identify each group and compute the center of gravity for each group
# this could get trickier if the bin size is not well chosen.

from scipy.ndimage.measurements import center_of_mass
levels=[]

for i in range(1,len(edges)-2):
    if  vals[i-1]==0 and vals[i]>0:
        istart=i
        #print 'istart: ',istart
    if  vals[i]>0 and vals[i+1]==0:
        istop=i
        #print 'istop', istop
        sum=np.sum(vals[istart:istop+1])
        c= center_of_mass(vals[istart:istop+1])[0] 
        level= edges[istart]+c*deltaedge
        levels.append(level)

        #print i, sum

print 'levels: ',levels        
print
print 'steps: ' ,np.diff(levels)

Output: enter image description here

Improve this answer

2 Comments

Do you know why it always omits the last level? There is clearly 10 groups in your histo but only 9 level array enteries -- Im having the same problem in other examples as well
It leaves out the first one. I think it's because there is no zero bin to the left. That could be fixed.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.