3

I would like to multiply the intensity of a EELS datacube with the thickness map. I tried with simple math command but I only obtain the result for the first slice. I think that the calculation is as simple as doing Result(x,y,E) = SI(x,y,E) * Thickness(x,y) and to loop on x and y, but I don't jnow how to do this with DM script. Thank you Regards

BmyGuest
  • 6,331
  • 1
  • 21
  • 35
Eric Leroy
  • 55
  • 5
  • In your eqauation you actually need to loop over all three variables (x,y and E), as each value of E represents one "image slice" of the stack. The first variant in the answer below does exactly that. If it would only be for one single value of E, you would just have a pixelwise multiplication of two image, which in DM-script is as easy as `A *= B`, that is what the second variant is doing (iterating over E). – BmyGuest Aug 01 '18 at 09:37

2 Answers2

1

There are different ways to do that. The following script creates the two test images for below: 

Image SIStack := RealImage( "Stack", 4, 32,32,64)
SIStack = iplane + 1 + random()
Image Map2D := RealImage( "Map", 4, 32, 32 )
Map2D = iradius * sin( irow/iwidth * pi() * 4 )

SIStack.ShowImage()
Map2D.ShowImage()

Variant 1

To multiply the two together, you can indeed use an expression very similar to what you have written. The intrinsic variables icol , irow and iplane (or more generally iindex(n)) are replaced by the respective pixel indices of X,Y and Z (starting with 0 for first index, mind you!). You can use them as the "running variables" in your expression. So the following will do what you want:

Image ResultStack := RealImage( "RStack", 4, 32,32,64)
ResultStack = SIStack[icol,irow,iplane] * Map2D[icol,irow]
ResultStack.ShowImage()

The above is the 1:1 translation of your equation and I typed it out to make that clear. Note that in this type of expression, at least one part of it needs to be of "static" size, i.e. not being indexed. This will define over which range the icol/irow/iplane indices will run. In the example above, the lefthand-side ResultStackis providing this, so icol goes from 0 to size-x of ResultStack minus 1.

As an asside: This even works if Map2D is smaller than the running variables! Any index "out of bounds" will be considered as the last available index, essentiallaly extrapolationg the "border" value.

Now the script above could be written more compactly (and without the addtional 'ResultsStack') as:

SIStack *= Map2D[icol,irow]

Variant 2

You can also use slicing and iterate over the z dimension:

for( number i = 0; i<64; i++ )
    SIStack.Slice2(0,0,i, 0,32,1 ,1,32,1 ) *= Map2D

For more details on either Variant, you might want to check out this answer

BmyGuest
  • 6,331
  • 1
  • 21
  • 35
0

Thank you. Here is the script I wrote. The only problem, is that the created datacube is not recognised as an EELS datacube.

image thick, si, res, out
number width, height, xsize, ysize, zsize,i
//select images
string title = "Image Selection Dialog"
string prompt = "Please select two images."
string label_si = "SI :"
string label_thick = "Thickness map:"
if ( !GetTwoLabeledImagesWithPrompt( prompt, title, label_si, si, label_thick, thick ) )
 Throw( "User pressed cancel." )
res=si
ImageCopyCalibrationFrom(res, si)
// size of the different images
thick.GetSize(width, height)
si.Get3DSize(xsize, ysize, zsize)
Result("width="+width+" pixels"+", height="+height+" pixels\n")
Result("SI xsize:"+xsize+" ysize:"+ysize+" Nb Energy channels:"+zsize+"\n")
for(number i=0; i<zsize; i++)
{
  res.Slice2(0,0,i,0,xsize,1,1,ysize,1) *= thick  
}
res.ShowImage()
Eric Leroy
  • 55
  • 5
  • This is because `image = image` just copies the values over and the SI info is stored in the meta data. It is often better to create a proper (deep) clone by using `image img := sourceImg.imageClone()` You then also don't need to copy over any calibrations. So just use `res := si.imageClone()` in your script above. – BmyGuest Aug 02 '18 at 18:21