如何以高效的方式使用 4D-STEM 数据集应用虚拟光圈？
How to apply virtual apperture with 4D-STEM dataset in EFFICIENT way?
我想将任意定义的位掩码应用为虚拟孔径，并以高效的方式将其应用于 4D-STEM 数据集。
我使用 SliceN 函数并逐像素应用掩码，这对于大型数据集来说非常慢。如何优化它以 运行 更快？
Image 4DSTEM := GetFrontImage() // dimention [ScanX, ScanY, Dx, Dy]
Image mask: = iradius // just an arbitrary mask (aperture)
Image out // dimention [ScanX, ScanY]
for (number i=0; i<ScanX; i++)
{ for (number j=0; j<ScanY; j++)
{
Diff2D = 4DSTEM.SliceN(4,2,i,j,0,0,2,Dx,1,3,Dy,1)
out.setpixel(i,j, sum(diff2D*mask))
}
}
out.showimage()
对于 [100,100,512,512] 数据集，需要几分钟才能完成。当我不得不多次重复该操作时，与矩阵操作相比，这是一种较慢的方式。但我不知道如何以有效的方式做到这一点。
谢谢！
您在这里遇到了脚本语言的限制。不幸的是，使用 sliceN 已经是您可以达到的最佳状态。速度优化中的其他一切都需要并行化的编译代码。 （即您可以编写 C++ 代码并使用 SDK 编译您自己的插件。）
但是，与您的示例相比还有一些改进空间。
首先，你上面的例子没有 运行 :c) 但很快就解决了。
第 1 点：
尽量避免数字类型转换。 DM 脚本只知道 number 但在内部正确的数字类型（整数、浮点数、signed/unsigned、字节大小）之间存在差异。脚本语言使用 real-4-byte 作为默认值，除非明确说明不同。有些方法默认会 return real-4-byte 。因此，如果数据和掩码都使用真实的 4 字节数据，处理速度会最快。
在我的测试中，运行 uint16 数据加 uint8 掩码 和 *real4 数据加 real4 掩码之间的时间差很明显！将近30%的时差
第 2 点：
请勿复制您的切片图像！ Dif2D.
使用 := 而不是 =
SliceN命令return是一个直接寻址所需内存的表达式。您可以直接在任何其他表达式中使用它（就像我在下面所做的那样），或者您可以使用 := 为其分配一个图像变量并为其命名。
速度提升不大，但每次循环迭代减少一次复制操作。
第 3 点：
你的额外知识：现在对于任意掩码你无能为力，但大多数情况下掩码在大范围内都是零值的，并且可以定义一个更小的包含所有非零点的 ROI。如果是这种情况，您可以将数学运算限制在该区域。
即不要将整个 DP 与相同大小的掩码相乘，只需使用较小的掩码并使用 DP 的相应子部分。
这实际上会有很大的不同，但这取决于您的面具。
当然要先“找到”这个ROI。在我下面的脚本中，我有一个辅助方法来做到这一点，利用相对较快的 max() 命令和图像旋转作为加速技巧。
第 4 点：
...将摆脱双循环并将其替换为图像表达式。不幸的是，DigitalMicrograph 目前 (GMS 3.3) 不支持 4D 或 5D 数据。
The script below executed on a [53 x 52 x 512 x 512] STEM DI (of real-4 byte data) gave me the following timings:
Original: 12.80910 sec
Test 1 : 10.77700 sec
Test 2 : 1.83017 sec
// Helper class for timing
class CTimer{
number s
string n
~CTimer(object self){result("\n"+n+": "+ (GetHighResTickCount()-s)/GetHighResTicksPerSecond()+" sec");}
object Start(object self, string n_) { n=n_; s=GetHighResTickCount(); return self;}
}
// Helper method to find best non-zero containing ROI
void GetNonZeroArea( image src, number &t, number &l, number &b, number &r )
{
image work = !!src // Make a binary image which is 0 only where src==0
number d
max(work,d,t) // get "first" non-zero pixel coordinate, this is y = dist from TOP
rotateRight(work) // rotate image right
max(work,d,l) // get "first" non-zero pixel coordinate, this is y = dist from LEFT
rotateRight(work) // rotate image right
max(work,d,b) // get "first" non-zero pixel coordinate, this is y = dist from BOTTOM
b = work.ImageGetDimensionSize(1) - b // Opposite side!
rotateRight(work) // rotate image right
max(work,d,r) // get "first" non-zero pixel coordinate
r = work.ImageGetDimensionSize(1) - r // Opposite side!
}
// The original proposed script (plus fixes to make it actually run)
image Original(image STEM4D, image mask)
{
Number ScanX = STEM4D.ImageGetDimensionSize(0)
Number ScanY = STEM4D.ImageGetDimensionSize(1)
Number Dx = STEM4D.ImageGetDimensionSize(2)
Number Dy = STEM4D.ImageGetDimensionSize(3)
Image out := RealImage("Test1",4,ScanX,ScanY)
for (number i=0; i<ScanX; i++)
{ for (number j=0; j<ScanY; j++)
{
image Diff2D = STEM4D.SliceN(4,2,i,j,0,0,2,Dx,1,3,Dy,1)
out.setpixel(i,j, sum(Diff2D*mask))
}
}
return out
}
// Remove copying the slice, just reference it
image Test1(image STEM4D, image mask)
{
Number ScanX = STEM4D.ImageGetDimensionSize(0)
Number ScanY = STEM4D.ImageGetDimensionSize(1)
Number Dx = STEM4D.ImageGetDimensionSize(2)
Number Dy = STEM4D.ImageGetDimensionSize(3)
Image out := RealImage("Test1",4,ScanX,ScanY)
for (number i=0; i<ScanX; i++)
{ for (number j=0; j<ScanY; j++)
{
image Diff2D := STEM4D.SliceN(4,2,i,j,0,0,2,Dx,1,3,Dy,1)
out.setpixel(i,j, sum(Diff2D*mask))
}
}
return out
}
// Limit mask size to what is needed!
image Test2(image STEM4D, image mask )
{
Number ScanX = STEM4D.ImageGetDimensionSize(0)
Number ScanY = STEM4D.ImageGetDimensionSize(1)
Number Dx = STEM4D.ImageGetDimensionSize(2)
Number Dy = STEM4D.ImageGetDimensionSize(3)
Image out := RealImage("Test1",4,ScanX,ScanY)
Number t,l,b,r
GetNonZeroArea(mask,t,l,b,r)
Number w = r - l
Number h = b - t
image subMask := mask.slice2(l,t,0, 0,w,1, 1,h,1 )
for (number i=0; i<ScanX; i++)
for (number j=0; j<ScanY; j++)
out.setpixel(i,j, sum(STEM4D.SliceN(4,2,i,j,l,t,2,w,1,3,h,1)*subMask))
return out
}
Image src := GetFrontImage() // dimention [ScanX, ScanY, Dx, Dy]
Number ScanX = src.ImageGetDimensionSize(0)
Number ScanY = src.ImageGetDimensionSize(1)
Number Dx = src.ImageGetDimensionSize(2)
Number Dy = src.ImageGetDimensionSize(3)
Number r = 50 // mask radius
Image maskImg := RealImage("Mask",4,Dx,Dy)
maskImg = iradius < r ? 1 : 0 // just an aperture mask
image resultImg
{
object timer = Alloc(CTimer).Start("Original")
resultImg := Original(src,maskImg)
}
resultImg.SetName("Oringal")
resultImg.ShowImage()
{
object timer = Alloc(CTimer).Start("Test 1")
Test1(src,maskImg).ShowImage()
}
resultImg.SetName("Test 1")
resultImg.ShowImage()
{
object timer = Alloc(CTimer).Start("Test 2")
Test2(src,maskImg).ShowImage()
}
resultImg.SetName("Test 2")
resultImg.ShowImage()
编译代码对比：
现在，应该补充一点，上面的脚本仍然相当慢。因为它是迭代和使用脚本语言。 DigitalMicrograph 完全编译的 c++ 代码 快 很多。因此，如果您拥有为您提供 SI 菜单的许可软件包，那么您需要使用 SI/Map/Signal 命令。对于我上面提到的示例 STEM DI，这几乎是瞬时的。我的 展示了如何通过脚本利用此功能。
正如我在 中提到的，使用经过编译的并行化代码时，真正的速度取胜。毕竟，DigitalMicrograph 在可用的 SI "signal" 地图功能中做到了这一点。此功能在免费版本中不可用，但如果您有 Spectrum-Imaging 采集，您很可能也有适当的许可证。
下面的答案通过使用命令 ChooseMenuItem() 访问 UI 并应用更多技巧来利用此功能。该脚本有点冗长，但它的各个部分还展示了一些其他值得了解的好技巧：
TestSignalIntegrationInSI is the main script demoing how things can work.CreatePickerByScript shows how one can create picker-spectra on SIs. This is used to open a 'Picker Diffraction Pattern' image from the STEM DI.AddTestMasksToDP_ROIs programmatically adds ROIs to the diffraction pattern to be used as mask AddTestMasksToDP_Threshold programmatically adds an intensity-threshold mask to be used as mask.AddTestMasksToDP_DPMasks programmatically adds the various types of diffraction-masks to be used as maskGetIntegratedSignalViaSIMenu is the central step of the script. With a picker-DP and required 'masks' on it front-most, the menu command is called to perform the signal-extraction (as fast as possible.) Then the displayed result-image is returned.GetNewestImage is just a utility method showing how on can access the latest memory-created image.
这是脚本：
image GetNewestImage()
{
// New images get the next higher imageID.
// This can be used to identify the "latest" created image.
if ( 0 == CountImages() ) Throw( "No image in memory!" )
// We create a temp. image to get the uppder limit
number lastID = RealImage("Dummy",4,1).ImageGetID()
// Then we search for the next lower existing one
image lastImg
for( number ID = lastID - 1; ID>0; ID-- )
{
lastImg := FindImageByID(ID)
if ( lastImg.ImageIsValid() ) break
}
return lastImg
}
image CreatePickerByScript( image SI, number t, number l, number b, number r )
{
if ( SI.ImageGetNumDimensions()<3 ) Throw( "Sorry, LineScans are not supprorted here." )
// Adding a non-volatile ROI of specific RoiNAME acts as if using
// the picker-tool. The ID string must be unique!
ROI pickerROI = NewROI()
pickerROI.RoiSetVolatile( 0 )
string uniqueID = GetDate(0)+"@"+GetTime(1)+";"+round(random()*1000)
pickerROI.RoiSetName( "SICursor(##"+uniqueID+"##)" )
SI.ImageGetImageDisplay(0).ImageDisplayAddROI( pickerROI )
// This creates the picker image.
// So the child is now the "newest" image in memory
image child := GetNewestImage()
return child
}
void AddTestMasksToDP_ROIs( image DP )
{
// Add ROIs to the DP which are your masks (any numebr and type of ROI works)
imageDisplay DPdisp = DP.ImageGetImageDisplay(0)
number dpX = DP.ImageGetDimensionSize(0)
number dpY = DP.ImageGetDimensionSize(1)
// Only simple RECT ROIs are supported
ROI maskRoi1 = NewROI()
maskRoi1.ROISetRectangle( dpY*0.1, dpX*0.1, dpY*0.8, dpX*0.3 )
DPdisp.ImageDisplayAddROI(maskRoi1)
// Arbitrary multi-vertex (use for ovals etc.)
ROI maskRoi2 = NewROI()
maskRoi2.ROISetRectangle( dpY*0.7, dpX*0.1, dpY*0.9, dpX*0.9 )
DPdisp.ImageDisplayAddROI(maskRoi2)
}
void AddTestMasksToDP_Threshold( image DP )
{
// Add intensity treshhold mask (highest 95% intensity range)
imageDisplay DPdisp = DP.ImageGetImageDisplay(0)
DPdisp.RasterImageDisplaySetThresholdOn( 1 )
number low = max(DP) * 0.05
number high = max(DP)
DPdisp.RasterImageDisplaySetThresholdLimits( low, high )
}
void AddTestMasksToDP_DPMasks( image DP )
{
// Add Diffraction masks to the DP
imageDisplay DPdisp = DP.ImageGetImageDisplay(0)
// Spot masks (always symmetric pair)
Component spotMask = NewComponent(8,0,0,0,0) // 8 = Spotmask
spotMask.ComponentSetControlPoint(4, 0, 0,0) // 4 = TopLeft of one spot [Size only]
spotMask.ComponentSetControlPoint(7,10,10,0) // 7 = BottomRight of one spot [Size only]
spotMask.ComponentSetControlPoint(8,150,0,0) // 8 = Spot position [center]
DPdisp.ComponentAddChildAtEnd(spotMask)
// Bandpass mask (Only circles are correctly supported)
Component bandpassMask = NewComponent(15,0,0,0,0) // 15 = Bandpass (ring)
number r1 = 100
number r2 = 120
bandpassMask.ComponentSetControlPoint(7,r1,r1,0) // 7 = BottomRight of one ring [Size only]
bandpassMask.ComponentSetControlPoint(14,r2,r2,0) // 14 = BottomRight of one ring [Size only]
DPdisp.ComponentAddChildAtEnd(bandpassMask)
// Wege mask (symmetric)
Component wedgeMask = NewComponent(19,0,0,0,0) // 19 = wedgemask (ringsegment)
wedgeMask.ComponentSetControlPoint(9,10,20,0) // 9 = One wedge vector
wedgeMask.ComponentSetControlPoint(10,-20,40,0) // 10 = Other wedge vector
DPdisp.ComponentAddChildAtEnd(wedgeMask)
// Array mask (symmetric)
Component arrayMask = NewComponent(9,0,0,0,0) // 9 = arrayMask (ringsegment)
arrayMask.ComponentSetControlPoint(9,-70,-60,0) // 9 = One array vector
arrayMask.ComponentSetControlPoint(10,99,-99,0) // 10 = Other array vector
arrayMask.ComponentSetControlPoint(4, 0, 0,0) // 4 = TopLeft of one spot [Size only]
arrayMask.ComponentSetControlPoint(7,20,20,0) // 7 = BottomRight of one spot [Size only]
DPdisp.ComponentAddChildAtEnd(arrayMask)
}
image GetIntegratedSignalViaSIMenu( image pickerChild )
{
// Call the Menu to do the work
// The picker-spectrum or DP needs to be front-most
pickerChild.SelectImage()
ChooseMenuItem("SI","Map","Signal")
// The created signal map is NOT the newest image
// (some internal iamges are created for the mask)
// but it is the front-most displayed one.
image signalMap := GetFrontImage()
return signalMap
}
image GetMaskFromSignalMap( image signalMap, number DPx, number DPy )
{
// The actual mask is stored in the tags
string tagPath = "Processing:[0]:Parameters:Mask"
tagGroup tg = signalMap.ImageGetTagGroup()
if ( !tg.TagGroupDoesTagExist(tagPath) )
Throw( "Sorry, no mask tag found." )
image mask := RealImage("Mask",4,DPx, DPy )
if ( !tg.TagGroupGetTagAsArray(tagPath,mask) )
Throw( "Sorry, could not retrieve mask. Maybe wrong size?" )
return mask
}
void TestSignalIntegrationInSI()
{
image STEMDI := GetFrontImage()
image DP := STEMDI.CreatePickerByScript(0,0,1,1)
if ( TwoButtonDialog( "Add ROIs as mask?", "Yes", "No" ) )
AddTestMasksToDP_ROIs( DP )
else if ( TwoButtonDialog( "Add intensity treshold as mask?", "Yes", "No" ) )
AddTestMasksToDP_Threshold( DP )
else if ( TwoButtonDialog( "Add diffraction masks as mask?", "Yes", "No" ) )
AddTestMasksToDP_DPMasks( DP )
image signalMap := GetIntegratedSignalViaSIMenu( DP )
number dpX = DP.ImageGetDimensionSize(0)
number dpY = DP.ImageGetDimensionSize(1)
// We may want to close the DP again. No longer needed
//DP.DeleteImage()
// Verification: Get Mask image form SignalMap
image usedMask := GetMaskFromSignalMap( signalMap, dpX, dpY )
usedMask.SetName( "This mask was used." )
usedMask.ShowImage()
}
TestSignalIntegrationInSI()
下面的解决方案通过执行就地乘法然后投影来利用内部表达式循环。
令人失望的是，事实证明解决方案实际上 比使用 SliceN 命令的 for 循环慢 。
For the same test-data of size [53 x 52 x 512 x 512] the resulting timing is:
Data copy: 1.28073 sec
Inplace multiply: 30.1978 sec
Project 1/2: 1.1208 sec
Project 2/2: 0.0019557 sec
InPlace multiplication with projections (total): 32.9045 sec
InPlace multiplication with projections (total): 34.9853 sec
// Helper class for timing
class CTimer{
number s
string n
~CTimer(object self){result("\n"+n+": "+ (GetHighResTickCount()-s)/GetHighResTicksPerSecond()+" sec");}
object Start(object self, string n_) { n=n_; s=GetHighResTickCount(); return self;}
}
image MaskMultipliedSum( image STEM4D, image MASK2D, number copyFirst )
{
// Boring feasability checks...
if ( 4 != STEM4D.ImageGetNumDimensions() )
Throw( "Input data is not 4D." )
if ( 2 != MASK2D.ImageGetNumDimensions() )
Throw( "Input mask is not 2D." )
Number ScanX = STEM4D.ImageGetDimensionSize(0)
Number ScanY = STEM4D.ImageGetDimensionSize(1)
Number Dx = STEM4D.ImageGetDimensionSize(2)
Number Dy = STEM4D.ImageGetDimensionSize(3)
if ( Dx != MASK2D.ImageGetDimensionSize(0) )
Throw ("X dimension of mask does not match input data." )
if ( Dy != MASK2D.ImageGetDimensionSize(1) )
Throw ("Y dimension of mask does not match input data." )
// Do the maths!
image workCopy4D
if ( copyFirst )
{
object timer = Alloc(CTimer).Start("Data copy")
workCopy4D = STEM4D
}
else
workCopy4D := STEM4D
{
object timer = Alloc(CTimer).Start("Inplace multiply")
workCopy4D *= MASK2D[idimindex(2),idimindex(3)]
}
// Now we want to "sum up" over Dx and Dy
image p1,p2
{
object timer = Alloc(CTimer).Start("Project 1/2")
p1 := project( workCopy4D, 3 )
}
{
object timer = Alloc(CTimer).Start("Project 2/2")
p2 := project( p1, 2 )
}
return p2
}
image stack4D, mask2D
If ( GetTwoLabeledImagesWithPrompt("Please select 4D data and 2D mask", "Select input", "4D data", stack4D, "2D mask", mask2D ) )
{
number doCopy = TwoButtonDialog("Create workcopy?","Yes (takes time)","No (overwrites input data!)")
object timer = Alloc(CTimer).Start("InPlace multiplication with projections (total)")
MaskMultipliedSum(stack4D,mask2D,doCopy).ShowImage()
}
我想将任意定义的位掩码应用为虚拟孔径，并以高效的方式将其应用于 4D-STEM 数据集。
我使用 SliceN 函数并逐像素应用掩码，这对于大型数据集来说非常慢。如何优化它以 运行 更快？
Image 4DSTEM := GetFrontImage() // dimention [ScanX, ScanY, Dx, Dy]
Image mask: = iradius // just an arbitrary mask (aperture)
Image out // dimention [ScanX, ScanY]
for (number i=0; i<ScanX; i++)
{ for (number j=0; j<ScanY; j++)
{
Diff2D = 4DSTEM.SliceN(4,2,i,j,0,0,2,Dx,1,3,Dy,1)
out.setpixel(i,j, sum(diff2D*mask))
}
}
out.showimage()
对于 [100,100,512,512] 数据集，需要几分钟才能完成。当我不得不多次重复该操作时，与矩阵操作相比，这是一种较慢的方式。但我不知道如何以有效的方式做到这一点。 谢谢！
您在这里遇到了脚本语言的限制。不幸的是，使用 sliceN 已经是您可以达到的最佳状态。速度优化中的其他一切都需要并行化的编译代码。 （即您可以编写 C++ 代码并使用 SDK 编译您自己的插件。）
但是，与您的示例相比还有一些改进空间。
首先，你上面的例子没有 运行 :c) 但很快就解决了。
第 1 点：
尽量避免数字类型转换。 DM 脚本只知道 number 但在内部正确的数字类型（整数、浮点数、signed/unsigned、字节大小）之间存在差异。脚本语言使用 real-4-byte 作为默认值，除非明确说明不同。有些方法默认会 return real-4-byte 。因此，如果数据和掩码都使用真实的 4 字节数据，处理速度会最快。
在我的测试中，运行 uint16 数据加 uint8 掩码 和 *real4 数据加 real4 掩码之间的时间差很明显！将近30%的时差
第 2 点：
请勿复制您的切片图像！ Dif2D.
:= 而不是 =
SliceN命令return是一个直接寻址所需内存的表达式。您可以直接在任何其他表达式中使用它（就像我在下面所做的那样），或者您可以使用 := 为其分配一个图像变量并为其命名。
速度提升不大，但每次循环迭代减少一次复制操作。
第 3 点：
你的额外知识：现在对于任意掩码你无能为力，但大多数情况下掩码在大范围内都是零值的，并且可以定义一个更小的包含所有非零点的 ROI。如果是这种情况，您可以将数学运算限制在该区域。
即不要将整个 DP 与相同大小的掩码相乘，只需使用较小的掩码并使用 DP 的相应子部分。
这实际上会有很大的不同，但这取决于您的面具。
当然要先“找到”这个ROI。在我下面的脚本中，我有一个辅助方法来做到这一点，利用相对较快的 max() 命令和图像旋转作为加速技巧。
第 4 点：
...将摆脱双循环并将其替换为图像表达式。不幸的是，DigitalMicrograph 目前 (GMS 3.3) 不支持 4D 或 5D 数据。
The script below executed on a
[53 x 52 x 512 x 512]STEM DI (of real-4 byte data) gave me the following timings:Original: 12.80910 sec Test 1 : 10.77700 sec Test 2 : 1.83017 sec
// Helper class for timing
class CTimer{
number s
string n
~CTimer(object self){result("\n"+n+": "+ (GetHighResTickCount()-s)/GetHighResTicksPerSecond()+" sec");}
object Start(object self, string n_) { n=n_; s=GetHighResTickCount(); return self;}
}
// Helper method to find best non-zero containing ROI
void GetNonZeroArea( image src, number &t, number &l, number &b, number &r )
{
image work = !!src // Make a binary image which is 0 only where src==0
number d
max(work,d,t) // get "first" non-zero pixel coordinate, this is y = dist from TOP
rotateRight(work) // rotate image right
max(work,d,l) // get "first" non-zero pixel coordinate, this is y = dist from LEFT
rotateRight(work) // rotate image right
max(work,d,b) // get "first" non-zero pixel coordinate, this is y = dist from BOTTOM
b = work.ImageGetDimensionSize(1) - b // Opposite side!
rotateRight(work) // rotate image right
max(work,d,r) // get "first" non-zero pixel coordinate
r = work.ImageGetDimensionSize(1) - r // Opposite side!
}
// The original proposed script (plus fixes to make it actually run)
image Original(image STEM4D, image mask)
{
Number ScanX = STEM4D.ImageGetDimensionSize(0)
Number ScanY = STEM4D.ImageGetDimensionSize(1)
Number Dx = STEM4D.ImageGetDimensionSize(2)
Number Dy = STEM4D.ImageGetDimensionSize(3)
Image out := RealImage("Test1",4,ScanX,ScanY)
for (number i=0; i<ScanX; i++)
{ for (number j=0; j<ScanY; j++)
{
image Diff2D = STEM4D.SliceN(4,2,i,j,0,0,2,Dx,1,3,Dy,1)
out.setpixel(i,j, sum(Diff2D*mask))
}
}
return out
}
// Remove copying the slice, just reference it
image Test1(image STEM4D, image mask)
{
Number ScanX = STEM4D.ImageGetDimensionSize(0)
Number ScanY = STEM4D.ImageGetDimensionSize(1)
Number Dx = STEM4D.ImageGetDimensionSize(2)
Number Dy = STEM4D.ImageGetDimensionSize(3)
Image out := RealImage("Test1",4,ScanX,ScanY)
for (number i=0; i<ScanX; i++)
{ for (number j=0; j<ScanY; j++)
{
image Diff2D := STEM4D.SliceN(4,2,i,j,0,0,2,Dx,1,3,Dy,1)
out.setpixel(i,j, sum(Diff2D*mask))
}
}
return out
}
// Limit mask size to what is needed!
image Test2(image STEM4D, image mask )
{
Number ScanX = STEM4D.ImageGetDimensionSize(0)
Number ScanY = STEM4D.ImageGetDimensionSize(1)
Number Dx = STEM4D.ImageGetDimensionSize(2)
Number Dy = STEM4D.ImageGetDimensionSize(3)
Image out := RealImage("Test1",4,ScanX,ScanY)
Number t,l,b,r
GetNonZeroArea(mask,t,l,b,r)
Number w = r - l
Number h = b - t
image subMask := mask.slice2(l,t,0, 0,w,1, 1,h,1 )
for (number i=0; i<ScanX; i++)
for (number j=0; j<ScanY; j++)
out.setpixel(i,j, sum(STEM4D.SliceN(4,2,i,j,l,t,2,w,1,3,h,1)*subMask))
return out
}
Image src := GetFrontImage() // dimention [ScanX, ScanY, Dx, Dy]
Number ScanX = src.ImageGetDimensionSize(0)
Number ScanY = src.ImageGetDimensionSize(1)
Number Dx = src.ImageGetDimensionSize(2)
Number Dy = src.ImageGetDimensionSize(3)
Number r = 50 // mask radius
Image maskImg := RealImage("Mask",4,Dx,Dy)
maskImg = iradius < r ? 1 : 0 // just an aperture mask
image resultImg
{
object timer = Alloc(CTimer).Start("Original")
resultImg := Original(src,maskImg)
}
resultImg.SetName("Oringal")
resultImg.ShowImage()
{
object timer = Alloc(CTimer).Start("Test 1")
Test1(src,maskImg).ShowImage()
}
resultImg.SetName("Test 1")
resultImg.ShowImage()
{
object timer = Alloc(CTimer).Start("Test 2")
Test2(src,maskImg).ShowImage()
}
resultImg.SetName("Test 2")
resultImg.ShowImage()
编译代码对比：
现在，应该补充一点，上面的脚本仍然相当慢。因为它是迭代和使用脚本语言。 DigitalMicrograph 完全编译的 c++ 代码 快 很多。因此，如果您拥有为您提供 SI 菜单的许可软件包，那么您需要使用 SI/Map/Signal 命令。对于我上面提到的示例 STEM DI，这几乎是瞬时的。我的
正如我在
下面的答案通过使用命令 ChooseMenuItem() 访问 UI 并应用更多技巧来利用此功能。该脚本有点冗长，但它的各个部分还展示了一些其他值得了解的好技巧：
TestSignalIntegrationInSIis the main script demoing how things can work.CreatePickerByScriptshows how one can create picker-spectra on SIs. This is used to open a 'Picker Diffraction Pattern' image from the STEM DI.AddTestMasksToDP_ROIsprogrammatically adds ROIs to the diffraction pattern to be used as maskAddTestMasksToDP_Thresholdprogrammatically adds an intensity-threshold mask to be used as mask.AddTestMasksToDP_DPMasksprogrammatically adds the various types of diffraction-masks to be used as maskGetIntegratedSignalViaSIMenuis the central step of the script. With a picker-DP and required 'masks' on it front-most, the menu command is called to perform the signal-extraction (as fast as possible.) Then the displayed result-image is returned.GetNewestImageis just a utility method showing how on can access the latest memory-created image.
这是脚本：
image GetNewestImage()
{
// New images get the next higher imageID.
// This can be used to identify the "latest" created image.
if ( 0 == CountImages() ) Throw( "No image in memory!" )
// We create a temp. image to get the uppder limit
number lastID = RealImage("Dummy",4,1).ImageGetID()
// Then we search for the next lower existing one
image lastImg
for( number ID = lastID - 1; ID>0; ID-- )
{
lastImg := FindImageByID(ID)
if ( lastImg.ImageIsValid() ) break
}
return lastImg
}
image CreatePickerByScript( image SI, number t, number l, number b, number r )
{
if ( SI.ImageGetNumDimensions()<3 ) Throw( "Sorry, LineScans are not supprorted here." )
// Adding a non-volatile ROI of specific RoiNAME acts as if using
// the picker-tool. The ID string must be unique!
ROI pickerROI = NewROI()
pickerROI.RoiSetVolatile( 0 )
string uniqueID = GetDate(0)+"@"+GetTime(1)+";"+round(random()*1000)
pickerROI.RoiSetName( "SICursor(##"+uniqueID+"##)" )
SI.ImageGetImageDisplay(0).ImageDisplayAddROI( pickerROI )
// This creates the picker image.
// So the child is now the "newest" image in memory
image child := GetNewestImage()
return child
}
void AddTestMasksToDP_ROIs( image DP )
{
// Add ROIs to the DP which are your masks (any numebr and type of ROI works)
imageDisplay DPdisp = DP.ImageGetImageDisplay(0)
number dpX = DP.ImageGetDimensionSize(0)
number dpY = DP.ImageGetDimensionSize(1)
// Only simple RECT ROIs are supported
ROI maskRoi1 = NewROI()
maskRoi1.ROISetRectangle( dpY*0.1, dpX*0.1, dpY*0.8, dpX*0.3 )
DPdisp.ImageDisplayAddROI(maskRoi1)
// Arbitrary multi-vertex (use for ovals etc.)
ROI maskRoi2 = NewROI()
maskRoi2.ROISetRectangle( dpY*0.7, dpX*0.1, dpY*0.9, dpX*0.9 )
DPdisp.ImageDisplayAddROI(maskRoi2)
}
void AddTestMasksToDP_Threshold( image DP )
{
// Add intensity treshhold mask (highest 95% intensity range)
imageDisplay DPdisp = DP.ImageGetImageDisplay(0)
DPdisp.RasterImageDisplaySetThresholdOn( 1 )
number low = max(DP) * 0.05
number high = max(DP)
DPdisp.RasterImageDisplaySetThresholdLimits( low, high )
}
void AddTestMasksToDP_DPMasks( image DP )
{
// Add Diffraction masks to the DP
imageDisplay DPdisp = DP.ImageGetImageDisplay(0)
// Spot masks (always symmetric pair)
Component spotMask = NewComponent(8,0,0,0,0) // 8 = Spotmask
spotMask.ComponentSetControlPoint(4, 0, 0,0) // 4 = TopLeft of one spot [Size only]
spotMask.ComponentSetControlPoint(7,10,10,0) // 7 = BottomRight of one spot [Size only]
spotMask.ComponentSetControlPoint(8,150,0,0) // 8 = Spot position [center]
DPdisp.ComponentAddChildAtEnd(spotMask)
// Bandpass mask (Only circles are correctly supported)
Component bandpassMask = NewComponent(15,0,0,0,0) // 15 = Bandpass (ring)
number r1 = 100
number r2 = 120
bandpassMask.ComponentSetControlPoint(7,r1,r1,0) // 7 = BottomRight of one ring [Size only]
bandpassMask.ComponentSetControlPoint(14,r2,r2,0) // 14 = BottomRight of one ring [Size only]
DPdisp.ComponentAddChildAtEnd(bandpassMask)
// Wege mask (symmetric)
Component wedgeMask = NewComponent(19,0,0,0,0) // 19 = wedgemask (ringsegment)
wedgeMask.ComponentSetControlPoint(9,10,20,0) // 9 = One wedge vector
wedgeMask.ComponentSetControlPoint(10,-20,40,0) // 10 = Other wedge vector
DPdisp.ComponentAddChildAtEnd(wedgeMask)
// Array mask (symmetric)
Component arrayMask = NewComponent(9,0,0,0,0) // 9 = arrayMask (ringsegment)
arrayMask.ComponentSetControlPoint(9,-70,-60,0) // 9 = One array vector
arrayMask.ComponentSetControlPoint(10,99,-99,0) // 10 = Other array vector
arrayMask.ComponentSetControlPoint(4, 0, 0,0) // 4 = TopLeft of one spot [Size only]
arrayMask.ComponentSetControlPoint(7,20,20,0) // 7 = BottomRight of one spot [Size only]
DPdisp.ComponentAddChildAtEnd(arrayMask)
}
image GetIntegratedSignalViaSIMenu( image pickerChild )
{
// Call the Menu to do the work
// The picker-spectrum or DP needs to be front-most
pickerChild.SelectImage()
ChooseMenuItem("SI","Map","Signal")
// The created signal map is NOT the newest image
// (some internal iamges are created for the mask)
// but it is the front-most displayed one.
image signalMap := GetFrontImage()
return signalMap
}
image GetMaskFromSignalMap( image signalMap, number DPx, number DPy )
{
// The actual mask is stored in the tags
string tagPath = "Processing:[0]:Parameters:Mask"
tagGroup tg = signalMap.ImageGetTagGroup()
if ( !tg.TagGroupDoesTagExist(tagPath) )
Throw( "Sorry, no mask tag found." )
image mask := RealImage("Mask",4,DPx, DPy )
if ( !tg.TagGroupGetTagAsArray(tagPath,mask) )
Throw( "Sorry, could not retrieve mask. Maybe wrong size?" )
return mask
}
void TestSignalIntegrationInSI()
{
image STEMDI := GetFrontImage()
image DP := STEMDI.CreatePickerByScript(0,0,1,1)
if ( TwoButtonDialog( "Add ROIs as mask?", "Yes", "No" ) )
AddTestMasksToDP_ROIs( DP )
else if ( TwoButtonDialog( "Add intensity treshold as mask?", "Yes", "No" ) )
AddTestMasksToDP_Threshold( DP )
else if ( TwoButtonDialog( "Add diffraction masks as mask?", "Yes", "No" ) )
AddTestMasksToDP_DPMasks( DP )
image signalMap := GetIntegratedSignalViaSIMenu( DP )
number dpX = DP.ImageGetDimensionSize(0)
number dpY = DP.ImageGetDimensionSize(1)
// We may want to close the DP again. No longer needed
//DP.DeleteImage()
// Verification: Get Mask image form SignalMap
image usedMask := GetMaskFromSignalMap( signalMap, dpX, dpY )
usedMask.SetName( "This mask was used." )
usedMask.ShowImage()
}
TestSignalIntegrationInSI()
下面的解决方案通过执行就地乘法然后投影来利用内部表达式循环。
令人失望的是，事实证明解决方案实际上 比使用 SliceN 命令的 for 循环慢 。
For the same test-data of size
[53 x 52 x 512 x 512]the resulting timing is:Data copy: 1.28073 sec Inplace multiply: 30.1978 sec Project 1/2: 1.1208 sec Project 2/2: 0.0019557 sec InPlace multiplication with projections (total): 32.9045 sec InPlace multiplication with projections (total): 34.9853 sec
// Helper class for timing
class CTimer{
number s
string n
~CTimer(object self){result("\n"+n+": "+ (GetHighResTickCount()-s)/GetHighResTicksPerSecond()+" sec");}
object Start(object self, string n_) { n=n_; s=GetHighResTickCount(); return self;}
}
image MaskMultipliedSum( image STEM4D, image MASK2D, number copyFirst )
{
// Boring feasability checks...
if ( 4 != STEM4D.ImageGetNumDimensions() )
Throw( "Input data is not 4D." )
if ( 2 != MASK2D.ImageGetNumDimensions() )
Throw( "Input mask is not 2D." )
Number ScanX = STEM4D.ImageGetDimensionSize(0)
Number ScanY = STEM4D.ImageGetDimensionSize(1)
Number Dx = STEM4D.ImageGetDimensionSize(2)
Number Dy = STEM4D.ImageGetDimensionSize(3)
if ( Dx != MASK2D.ImageGetDimensionSize(0) )
Throw ("X dimension of mask does not match input data." )
if ( Dy != MASK2D.ImageGetDimensionSize(1) )
Throw ("Y dimension of mask does not match input data." )
// Do the maths!
image workCopy4D
if ( copyFirst )
{
object timer = Alloc(CTimer).Start("Data copy")
workCopy4D = STEM4D
}
else
workCopy4D := STEM4D
{
object timer = Alloc(CTimer).Start("Inplace multiply")
workCopy4D *= MASK2D[idimindex(2),idimindex(3)]
}
// Now we want to "sum up" over Dx and Dy
image p1,p2
{
object timer = Alloc(CTimer).Start("Project 1/2")
p1 := project( workCopy4D, 3 )
}
{
object timer = Alloc(CTimer).Start("Project 2/2")
p2 := project( p1, 2 )
}
return p2
}
image stack4D, mask2D
If ( GetTwoLabeledImagesWithPrompt("Please select 4D data and 2D mask", "Select input", "4D data", stack4D, "2D mask", mask2D ) )
{
number doCopy = TwoButtonDialog("Create workcopy?","Yes (takes time)","No (overwrites input data!)")
object timer = Alloc(CTimer).Start("InPlace multiplication with projections (total)")
MaskMultipliedSum(stack4D,mask2D,doCopy).ShowImage()
}