VI_OnDiffraction.md

1. キャリブレーション

基本的なキャリブレーション方法

正しい解析を行うためには、回折図形のキャリブレーション(カメラ長)が正しくなくてはいけません。STEM像やTEM像と異なり回折図形には中心があるのでscale(1画素当たりの単位)だけでなくorigin(中心座標)も正しく行う必要があります。

手入力による設定
Ctrl+Dでimage tagを開き、Calibrationでx軸とy軸の原点[pixel]とscale[nm-1/pixel]を設定します。

ROIを使う方法
画像中の参考となる長さの部分にlineのROIを描き、メニューのAnalysis > Calibrate... で単位と数値を入力します。

キャリブレーションされたかどうかの確認方法
マウスを動かしてimage status(左上)から、X: Y:を確認すると、キャリブレーションした座標が単位付きで表示されます。あるいは、Ctrl+Dでimage tagを開いたのち、Captionsにチェックを入れると、データの外側に軸とscaleが表示されます。

スクリプトを使ったキャリブレーション方法

ScriptでROI(Region Of Interest)を使い設定する方法
下記のscriptを使うと、ROIで回折図形の中心にpoint(+)をおいたり、あるいは同じ面間隔の回折スポットを通るようにoval(〇)を描き、原点を決めることができます。下記のscriptは丸ごとコピーペーストして使えます。if ～ else ifの行はエレガントではありませんが、ROIの種類毎に後々コピペすることも考え、繰り返し(repeat myself)しています。
ImageDisplayCountROIS( )を使ってROIを数えているので、for loopとImageDisplayGetROI( )で全て見ても良いのですが、自分で使う分には数を確認するだけで十分と考えて簡素にしています。

スクリプト例1 ROIを使ったキャリブレーション

// Origin calibration using oval/rectangle/line/point ROI
image 			img :=getfrontimage()
imageDisplay	imgDsp = img.ImageGetImageDisplay( 0 )
ImageDisplaySetCaptionOn( imgDsp, 1 )   // Captionをつけて原点を分かりやすくします。
// ROIが１つであることを確認します。
      number n_ROI = imgDsp.ImageDisplayCountROIS() 
        if(n_ROI ==0) {OkDialog("Put oval/point/rectangle/line ROI\nCalibration is aborted."); exit(0)
    			}
        if(n_ROI !=1) {OkDialog("Put only one ROI.\nCalibration is aborted."); exit(0)
    			}
// Variables
  number centerX, centerY	
  number top,left, bottom, right		// oval, rectangle
  number startX, startY, endX, endY	// line
  number pointX, pointY				      // point
  string unit
  number originalX, originalY 
  number scaleX, scaleY
// 現在のcalibrationの情報を得ます。
  ImageGetDimensionCalibration(img, 0, originalX, scaleX, unit, 1 )
  ImageGetDimensionCalibration(img, 1, originalY, scaleY, unit, 1 )
// ROIを読み出します
  result("\n\nROI coordinates:")
	ROI crrROI = ImgDsp.ImageDisplayGetROI( 0 )
if (ROIisPoint(crrROI)){ 	        ROIGetPoint(crrROI, pointX, pointY) ; result("\tPoint")
	     centerX = pointX			; centerY = pointY
	     result("\n\t(pointX,pointY):("+pointX+" ,"+pointY+")")
	     }
  else if (ROIisOval(crrROI)){		ROIGetOval(crrROI, top, left, bottom, right) ; result("\tOval")	
	     centerX = (right  + left)/2; centerY = (bottom + top )/2
	     result("\n\t(top,left,bottom,right):("+top+" ,"+left+" ,"+bottom+" ,"+right+")")
	     }
  else if (ROIisRectangle(crrROI)){ ROIGetRectangle(crrROI,top,left,bottom,right);result("\tRectangle")
	     centerX = (right  + left)/2; centerY = (bottom + top )/2
	     result("\n\t(top,left,bottom,right):("+top+" ,"+left+" ,"+bottom+" ,"+right+")")
	     }
  else if (ROIisLine(crrROI)){ 	    ROIGetline(crrROI, startX, startY, endX, endY) ; result("\tLine")
  	centerX = (startX + endX)/2; centerY = (startY + endY)/2
	     result("\n\t(startX,Y,endX,Y):(" + startX +" ," + startY +" ," +endX+" ,"+ endY + ")")
	     }
  else {
    OkDialog("Loop/Curve ROI cannot be used.\nPut oval/point/rectangle/line ROI.\nCalibration is aborted.")
    exit(0)
    }
result("\n\t(centerX, centerY)     :(" + centerX +" ," + centerY +" )")
if(centerX==0 || centerY==0){
    if(!OkCancelDialog("Origin is set to (" + centerX + ", "+ centerY + " ). ")) exit(0)
    }
// Calibration
Result("\n written by Kimoto")
ImageSetDimensionCalibration(img, 0, centerX , scaleX, unit, 1 )	//ここでx軸
ImageSetDimensionCalibration(img, 1, centerY , scaleY, unit, 1 )  //ここでy軸を設定します。
result("\nOriginal  origin (x,y): " + originalX + " , " + OriginalY )  
result("\nCorrected origin (x,y): " + centerX + " , " + centerY )  
Result("\n\t Written by Koji Kimoto") //一応著作権を宣言しておきます。

ROIの座標はintegerですが、originは小数点以下で指定できます。
次のようなスクリプトで、originを少しずつ動かして確認することもできます。数字の4(左)、6(右)、8(上)、2(下)で、動くステップは+(2倍)、-(半分)、/で元のステップ値に戻し、スペースバーで決定します。例えば、ROIのBand Pass tool(ドーナッツのようなマーク)をおいておくと、原点を中心とする円を描きますのでそれを見ながらスポット位置の調整ができます。

スクリプト例2 originの微調整

  image img:=getfrontimage()
    string unit
    number originalX, originalY, scaleX, scaleY
    ImageGetDimensionCalibration(img, 0, originalX, scaleX, unit, 1 )
    ImageGetDimensionCalibration(img, 1, originalY, scaleY, unit, 1 )
    result("\nOriginal  origin (x,Y): " + originalX + " , " + originalY )
    number preX = originalX
    number preY = originalY
    number OriginalStrength =1
    number strength = OriginalStrength
    number key
    OpenAndSetProgressWindow("Aligh Origin","( "+ preX +" , "+preY + " )  step: " + strength ,"Space to abort")
    while(key!=32){ // Space bar (key=32)が押されるまで続けます。
      	key=GetKey()
        if(key ==56) preY -= strength			// up      8  
        if(key ==50) preY += strength			// down    2  
        if(key ==52) preX -= strength			// left    4  
        if(key ==54) preX += strength			// right   6  
        if(key ==47){					//  Back to original origin(x,y)  5
          preY = originalX
          preY = originalY
          }
        if(key ==43) strength *= 2                //  +
        if(key ==45) strength /= 2                //  -
        if(key ==42) strength = OriginalStrength  //  /
    OpenAndSetProgressWindow("Aligh Origin","( "+preX+" , "+preY+" )  step: "+strength,"Space to abort")  
      ImageSetDimensionCalibration(img, 0, preX , scaleX, unit, 1 )  
      ImageSetDimensionCalibration(img, 1, preY , scaleY, unit, 1 )  
      UpdateImage(img)  
      }
    result("\nCorrected origin (x,y): " + preX + " , " + preY )
    CloseProgressWindow()
    Result("\n\t Written by Koji Kimoto")

2. 回転平均の計算

回転平均の計算
回折強度を散乱角に対して1Dプロットするために、回転平均を取りたい場面が多くあります。画像データの中心でrotational averageを行うscriptは、DigitalMicrographのexample scriptとしてGatanから提供(Example Script folder参照)されています。元々のscriptはTEM用CCDカメラの実用化を実現したGatanのグルの一人であるDr. Ming Panが書いたものです。
ここではそのscriptを参考に、画像のoriginを中心に回転平均を取り1Dデータ表示するscriptを示します。通常、新たな画像を定義すると原点は画像の中心ではなく左上の(0,0)となっているので、その場合には画像の中心を原点とするif文を挿入しています。GetSize(),GetOrigin(),CreateFloatImage()などのコマンド群には歴史を感じますが、動きますのでそのままにしてあります。

スクリプト例3 Off-centered rotational averaging

image OffCentRotAve1D(image img)
{
number sizeX,sizeY
GetSize(img, sizeX,sizeY)
number oriX, oriY
GetOrigin(img, oriX, oriY)
	if(oriX==0 && oriY==0){
	oriX = sizeX/2; oriY = sizeY/2; result("\nOrigins X,Y are set to the center of the front image.")
	}
number newsizeX = 2 * min(sizeX-oriX, oriX)
number newsizeY = 2 * min(sizeY-oriY, oriY)
number newsizeXY = min(newsizex, newsizeY)
number halfminor = newsizeXY/2
number centerX   = newsizeX/2
number centerY   = newsizeY/2
number sampling = halfminor*8	// Why 8? Note aliasing of Nyquist frequency in case of low sampling
image band
Band := CreateFloatImage( "map", halfMinor, sampling )
number k = 2 * pi() / sampling
Band = warp(img, icol*sin(irow*k) + OriX, icol*cos(irow*k) + OriY ) // bug fix
image proj := CreateFloatImage( "line projection", halfMinor, 1 )
proj =0
proj[icol,0] += band
proj /= sampling
return proj
Result("\nWritten by Koji Kimoto")  
}

// main
image img:=getfrontimage()
image img1D
img1D := OffCentRotAve1D(img)
showimage(img1D)

3. [1/nm]から[mrad]への変換

通常、画像取得装置では、面間隔[nm]への変換が簡単になるように、回折図形は[1/nm]でキャリブレーションされています。一方、収束角などの実験条件は多くの場合には角度[mrad]で指定されています。しばしば実験条件を[mrad]で確認したいためこのscriptを作りました。
加速電圧(V)によって波長(λ)が変わるため、散乱角2θと面間隔の逆数(1/d)との関係も変わります(2θ＝λ(V)/d)。そのため変換には加速電圧情報が必要です。通常の実験データでは加速電圧はimage tag（"Microscope Info:Voltage"）に保存されています。次のscriptでは、[1/nm]を[mrad]に変換し、ファイル名に_mradをつけます(生データを間違って上書きしないようにするため)。加速電圧のtag情報が無い場合には、GetNumber()で入力を求めます。

スクリプト例4 Conversion from 1/nm to mrad

number kVLamNM(number KV)
{
number LAM		//	wavelength in nm
number RelVol   // relativistic V 
number Voltage	//  V 
Voltage = kV * 1000 
RelVol = voltage*(1+0.000000978502*voltage)
LAM = 1.2263/sqrt(RelVol)	// in nm
return LAM
}

image img				  
number scaleX, scaleY	// original scale
number newScaleX, newScaleY	// converted scale
string unitstring	// unitstring, e.g., "1/nm" and "mrad"
number V		
number kV		
number lamNM	// wavelength in nm

img := getfrontimage()
  getscale(img, scaleX, scaleY)
if(!getnumbernote(img,"Microscope Info:Voltage" ,V)) getnumber("Input acc. voltage [V]", V, V)
  kV = V / 1000
  lamNM = kVLamNM(KV)

GetUnitString(img, unitstring)
  if(unitstring!="1/nm") 
  {
  OkDialog("Scale unit should be (1/nm)")
  exit(0)
  }

NewScaleX = 1000 * ScaleX * lamNM	
SetName(img, getname(img)+"_mrad")
SetUnitString(img, "mrad")
SetScale(img,NewScaleX, newScaleX) 
Result("\nWritten by Koji Kimoto")

4. 4D-STEMのデータの座標の順番

4D-STEMはプローブを走査しながらdiffraction patternを取得するもので、位置x,yと回折図形の座標u,vの4つの座標のデータが取得されるので4D-STEMと呼ばれています。先日とある会議で「4D-STEMって誰が言い出したの？」とDr. Colin Ophusに尋ねたところ、「私です」と言っていましたので、命名はDr. Ophusだと思います。木本は2011年のUltramicroscopyの論文で、spatially-resolved diffractometryと名前を付けたのですが普及しませんでした。
　4D-STEM dataは古いバージョンと新しいバージョンで座標の順番が変わっています。
　当初は特定バージョンのバグかと思ったのですが、その後、石塚和夫先生(HREM Research Inc.)と話をしたところ、GMS3のある時点でx,y,u,vがu,v,x,yに変わったと教えていただきました。その話を、早稲田大学の平田先生にしたところ「そうなんだよ、最初は分からなかったよ」と平田先生もおっしゃっていました。ユーザー同士の情報交換は重要だと思います。
　xyとuv座標の入れ替えは、以前はsliceN(4,4,...)とMeta Data:Data Order Swappedから自作scriptで行っていましたが、GMS3.4以降ではVolume > Reorder for fast spatial/diffraction dimension access というメニューがあり、簡単に変換できます。以前はscriptでしかできなかったことが、いつの間にかメニューから標準的な機能としてできるようになっている事はしばしばあります。例えば10年前、Volume > Rotate は無かったので自作scriptで処理していました。当時はFourier transform (FT)も2Dしかできなかったので、2D-FTと1D-FT(complex)とのloopを回していたのですが、まもなくHREM Research Inc.が提供するIPU(フリーライブラリー)で高速にできるようになり、その後のDigitalMicrographでは標準機能となりました。
　DM上のpythonでデータ解析をする際には、imageをNumPy arrayにしています。NumPyは.reshape()や.transpose()が便利で使いやすい(計算も速い)のですが、NumPyの座標は3次元の場合z, y, xの順番なのでDMデータの変換時に混乱します。Gatanのhelp fileにその記述があるのに気が付いたのは、ひとしきり試行錯誤した後でした。

VI. Diffraction解析に使う小さなツール群

1. キャリブレーション

基本的なキャリブレーション方法

スクリプトを使ったキャリブレーション方法

2. 回転平均の計算

3. [1/nm]から[mrad]への変換

4. 4D-STEMのデータの座標の順番