如何在 SI 图像中移动列

编程问答 2022-05-17

问题描述

为了在SI图像中进行漂移校正，如下图所示：

drift corrected SI

我写代码：

number max_shift=5
image src := GetFrontimage()
number sx,sy,sz
src.Get3DSize(sx,sz)
result("sx: "+sx+"\n")
result("sy: "+sy+"\n")
result("sz: "+sz+"\n")

// assume a random shift in x
image shift := IntegerImage("xcorrd",4,1,sy)
shift = max_shift*Random() 

// make a coordinate table
image col := IntegerImage("col",sx,sy)
image row := IntegerImage("row",sy)
image plane := IntegerImage("plane",sy)
col = icol
row = irow
plane = iplane

// to expand the shift as the same size with source image
image ones := IntegerImage("ones",sy)
ones = 1

// create a random column shift of the source SI image
for(number i=0; i<sy; i++) {
    col[i,i+1,sx] = col[i,sx]+shift.GetPixel(0,i)*ones[i,sx]
};

// drift corrected
image im := RealImage("test si",sx+max_shift,sz)
im=0
im[col,row,plane] = src[icol,irow,iplane]

im.ImageGetTagGroup().TagGroupcopyTagsFrom(src.ImageGetTagGroup())
im.ImagecopyCalibrationFrom(src)
im.SetName(src.GetName()+"-drift corrected")
im.showimage()

图像可以被校正，但是光谱不能转移到校正后的 SI，如图所示：

the result

我只是想知道我的脚本有什么问题。

提前致谢。

解决方法

im[col,row,plane] = src[icol,irow,iplane]

内在变量 icol、irow、iplane 将通过行中唯一的固定大小 图像表达式求值。在您的情况下，col、row 和 plane（大小相同）

但是，它们都是 2D 的，因此内部发生的是您迭代 X 和 Y，然后写入值：

im[ col(x,y),row(x,plane(x,y) ] = src[x,y,0] // iterated over all x/y

正如 Don I 在评论中提到的，您可能想要遍历 z 维度。

或者，您可以在脚本中制作所有大小为 (sx,sy,sz) 的图像。 这对表达式有效，但效率极低。

一般来说，这里最好的解决方案是根本不使用 icol、irow、iplane，而是使用 Slice 命令。见this answer：

我可能会为如下所示的 SI 编写按行 x-shift 的代码： 该脚本利用了一个事实，即可以在 x 方向上移动整个“块”（X x 1 x Z），迭代 y。

number sx = 256
number sy = 256
number sz = 100
image testSI := realImage("SI",4,sx,sz)
testSI = sin(itheta/(idepth-iplane)*idepth) + (iplane % (icol+1))/idepth
testSI.ShowImage()

image xdrift := RealImage("XDrift per line",sy)
xdrift = trunc(random()*5 + 20*sin(icol/iwidth*3*PI()))
xdrift.ShowImage()

// Apply linewise Drift to SI,assuming xDrift holds this data
xDrift -= min(xDrift)   // ensure only positive shifts
image outSI := realImage("SI shifted",sx+max(xDrift),sz)
outSI.ShowImage()

for( number y=0; y<sy; y++ ){
    number yShift = sum(xDrift[y,0])
    outSI.slice2( yShift,1,2,sz,1 ) = testSI.slice2(0,1)
}

下面的脚本执行“逐平面”迭代，但对逐平面移动没有限制。 事实上，这里每个像素都有一个指定的 XY 位移。

请注意，如果您想使用值的双线性插值（以及有效范围之外的 0 截断），您可以使用 warp(source,xexpr,yexpr ) 代替二维寻址 source[ xexpr,yexpr ]。

number sx = 256
number sy = 256
number sz = 100
image testSI := realImage("SI",sz)
testSI = sin(itheta/(idepth-iplane)*idepth) + (iplane % (icol+1))/idepth
testSI.ShowImage()

image xdrift := RealImage("pixelwise XDrift",sy)
xdrift = irow%10*random() + 20*cos(irow/iheight*5*PI())
xdrift.ShowImage()

image ydrift := RealImage("pixelwise yDrift",sy)
ydrift = 10*abs(cos(icol/iwidth* (irow/iheight) * 10 * PI())) + irow/iheight * 10 
ydrift.ShowImage()

// Apply pixelwise Drift to SI
xDrift -= min(xDrift)   // ensure only positive shifts
yDrift -= min(yDrift)   // ensure only positive shifts
number mx = max(xDrift)
number my = max(yDrift)
image outSI := realImage("SI shifted",sx+mx,sy+my,sz)
outSI.ShowImage()
for( number z=0;z<sz;z++){
    image outPlane := outSI.Slice2( 0,z,1) 
    image srcPlane := testSI.Slice2( 0,1)
    outPlane = srcPlane[ icol + xdrift[icol,irow] - mx,irow + ydrift[icol,irow] - my ]
    // outPlane = srcPlane.warp( icol + xdrift[icol,irow] - my)
}
dm-script