Skip to content

Roi axis stack #1219

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
sergey-yaroslavtsev wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Roi axis stack #1219

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
fc53f5a
extra elif for flatten data
sergey-yaroslavtsev May 11, 2026
8acc5eb
Scale are scalar not full array
sergey-yaroslavtsev May 11, 2026
9f592b3
add protectin, cleaning, add suggested dimensions to dialog
sergey-yaroslavtsev May 12, 2026
b2e7126
some renaming and comment correction
sergey-yaroslavtsev May 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion src/PyMca5/PyMcaGui/pymca/QStackWidget.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -246,7 +246,8 @@ def setStack(self, *var, **kw):
if (1 in self._stack.data.shape) and\
isinstance(self._stack.data, numpy.ndarray):
oldshape = self._stack.data.shape
dialog = ImageShapeDialog(self, shape=oldshape[0:2])
suggested = self._stack.info.get("SuggestedImageShape", oldshape[0:2])
dialog = ImageShapeDialog(self, shape=suggested)
dialog.setModal(True)
ret = dialog.exec()
if ret:
Expand Down
83 changes: 77 additions & 6 deletions src/PyMca5/PyMcaIO/HDF5Stack1D.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -330,7 +330,7 @@ def loadFileList(self, filelist, selection, scanlist=None):
if self.__dtype0 is None:
if (bytefactor == 8) and (neededMegaBytes < (2*physicalMemory)):
# try reading as float32
print("Forcing the use of float32 data")
_logger.info("Forcing the use of float32 data")
self.__dtype = numpy.float32
else:
raise MemoryError("Force dynamic loading")
Expand Down Expand Up @@ -428,7 +428,7 @@ def loadFileList(self, filelist, selection, scanlist=None):
else:
del mcaObjectPaths["elapsed_time"]
else:
# we have to have as many elpased times as MCA spectra
# we have to have as many elapsed times as MCA spectra
_time = numpy.zeros((self.data.shape[0] * self.data.shape[1]),
numpy.float32)
if "calibration" in mcaObjectPaths:
Expand Down Expand Up @@ -800,13 +800,47 @@ def loadFileList(self, filelist, selection, scanlist=None):
# try to get scales
scaleList = []
if xSelectionList is not None:
# McaIndex not always eqaul mcaIndex
Comment thread
woutdenolf marked this conversation as resolved.
Outdated
stackMcaIndex = self.info.get("McaIndex", mcaIndex)
if len(xDatasetList) == 1:
xDataset = xDatasetList[0]
if xDataset.size == shape[self.info['McaIndex']]:
# assuming providing channels
self.x = [xDataset.reshape(-1)]
else:
# user could not set coordiantes for one axis only
_logger.warning("Ignoring channels selection %s" % xSelectionList)

# for dataset of N*M self.data.shape = (1, N, M)
# if positions are flatten then both X, Y, coordinates should be of size N.
elif len(xDatasetList) in (2, 3) and len(xDatasetList[0]) == len(xDatasetList[1]) == self.data.shape[1]:
# assuming providing spatial coordinates (and maybe channels)
origin_x, origin_y = self.shortenScales(xDatasetList[0], xDatasetList[1])
Comment thread
woutdenolf marked this conversation as resolved.
Outdated
if origin_x.size > 1:
delta_x = numpy.mean(origin_x[1:] - origin_x[:-1], dtype=numpy.float32)
else:
delta_x = 1.0
xScale = [origin_x[0], delta_x]

if origin_y.size > 1:
delta_y = numpy.mean(origin_y[1:] - origin_y[:-1], dtype=numpy.float32)
else:
delta_y = 1.0
yScale = [origin_y[0], delta_y]

scaleList = [xScale, yScale]
self.info["SuggestedImageShape"] = (len(origin_y), len(origin_x))

if len(xDatasetList) == 2:
pass
elif len(xDatasetList) == 3 and xDatasetList[2].size == shape[self.info['McaIndex']]:
# assuming providing spatial coordinates and channels
self.x = [xDatasetList[2].reshape(-1)]
Comment thread
woutdenolf marked this conversation as resolved.
else:
_logger.warning("Coordinates were applied, but channels could not be")
_logger.warning("Ignoring channels selection %s" % xSelectionList)


elif len(xDatasetList) == len(self.data.shape):
# assuming providing spatial coordinates and channels
goodScale = 0
Expand All @@ -823,7 +857,7 @@ def loadFileList(self, filelist, selection, scanlist=None):
for i in range(len(self.data.shape)):
dataset = xDatasetList[i].reshape(-1)
datasize = self.data.shape[i]
if i == mcaIndex:
if i == stackMcaIndex:
self.x = [dataset]
else:
origin = dataset[0]
Expand All @@ -842,8 +876,8 @@ def loadFileList(self, filelist, selection, scanlist=None):
_logger.warning("Ignoring dimension selections %s" % xSelectionList)
elif len(xDatasetList) == (len(self.data.shape) - 1):
scaleList = []
for i in range(len(self.data.shape)):
if i == mcaIndex:
for i in range(len(xDatasetList)):
if i == stackMcaIndex:
continue
dataset = xDatasetList[i].reshape(-1)
datasize = self.data.shape[i]
Expand Down Expand Up @@ -901,7 +935,7 @@ def loadFileList(self, filelist, selection, scanlist=None):
if len(dims) == len(self.data.shape):
scaleList = []
for i in range(len(self.data.shape)):
if i == mcaIndex:
if i == self.info['McaIndex']:
Copy link
Copy Markdown
Collaborator

@woutdenolf woutdenolf May 12, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

mcaIndex, self.info['McaIndex'] and stackMcaIndex.

I'm going to assume you know what you're doing because I have no idea. Same for all the other places where you changed it.

Copy link
Copy Markdown
Collaborator Author

@sergey-yaroslavtsev sergey-yaroslavtsev May 13, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

mcaIndex is initial index of channels. In most cases it goes to self.info["McaIndex"] to be stored but in some specific cases self.info["McaIndex"] stores not the mcaIndex but exact value.
line 459 in HDF5Stack1D.py for example:

        if (not DONE) and (not considerAsImages):
            _logger.info("Data in memory as spectra")
            self.info["McaIndex"] = 2
            n = 0

So i believe it is what should be used. But in most cases they are the same.

Copy link
Copy Markdown
Collaborator

@woutdenolf woutdenolf May 14, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How is 2 the exact value? Looks like an index to me, although it seems very random to suddenly choose 2.

Copy link
Copy Markdown
Collaborator Author

@sergey-yaroslavtsev sergey-yaroslavtsev May 20, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, probably not the best wording on my side. I meant "specific value/index" instead of referring to mcaIndex.

Copy link
Copy Markdown
Member

@vasole vasole May 20, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So i believe it is what should be used. But in most cases they are the same.

Stop believing.

info['McaIndex'] can be -1 to indicate any arbitrarily shaped dataset containing spectra at its last dimension.

Because of that McaIndex is recalculated and used based on the shape of the dataset. It will be 2 if your dataset is a three-dimensional array of spectra.

Again, think generic and not particular.

Copy link
Copy Markdown
Collaborator Author

@sergey-yaroslavtsev sergey-yaroslavtsev May 20, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I’m not sure how that contradicts what I said. So was it actually not a bug? If that’s the case, could you please explain why?

if initially mcaIndex was 1 then considerAsImages is False and DONE could be False.
Therefore, self.info["McaIndex"] = 2.
Thus, using mcaIndex (=1) instead self.info["McaIndex"] (=2) can lead to unwanted behavior.

continue
dataset = dims[i]
origin = dataset[0]
Expand All @@ -921,6 +955,43 @@ def loadFileList(self, filelist, selection, scanlist=None):
self.info["xScale"] = xScale
self.info["yScale"] = yScale

def shortenScales(self, A, B):
A = numpy.asarray(A)
B = numpy.asarray(B)

if len(A) == 0 or len(B) == 0:
raise ValueError("Arrays must not be empty")

# Detect which array is related to the slow motor
if len(B) == 1 or (len(A)>1 and abs(A[1] - A[0]) > abs(B[1] - B[0])):
fast = A
slow = B
fast_is_A = True
else:
fast = B
slow = A
fast_is_A = False

# Find repetition length n in the slow array
diff_idx = numpy.where(slow != slow[0])[0]
Copy link
Copy Markdown
Collaborator

@woutdenolf woutdenolf May 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To be clear this assumes the slow motor readout is not noisy.

Copy link
Copy Markdown
Collaborator Author

@sergey-yaroslavtsev sergey-yaroslavtsev May 13, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed i thought it is most regular case, but it appears to be 50-50.
I need to fix it.

Copy link
Copy Markdown
Collaborator

@woutdenolf woutdenolf May 14, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes but the route of analyzing noise motors of regular scans to find the regular grid has been tried many times and is inherently flaky.

Your next implementation might involve some kind of noise threshold which of course is impossible to know.

Then you might realize there is only one property you can rely on: the values of the fast axis are "piecewise monotonic". Est is using that to split black-and-forth energy scans, but that's only a 1D problem. You need to first find the fast axis. Silx is using this for silx view.

And after all this work you realize you need to handle CTRL-C etc.

In my opinion all this is a dead end and we should not rely on regular vs. irregular. ewoksfluo uses an optimal grid finding approach for any list of nD coordinates.

# Protection if one motor did not move at all
if len(diff_idx) == 0:
n = len(slow)
else:
n = diff_idx[0]

short_fast = fast[:n]
short_slow = slow[::n]
Comment on lines +982 to +983
Copy link
Copy Markdown
Collaborator

@woutdenolf woutdenolf May 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We sample the coordinates and take this as the final regular grid coordinates. To make it more robust we could take the median of all possible samplings. Not sure we care. All this is very approximate anyway.

Copy link
Copy Markdown
Collaborator Author

@sergey-yaroslavtsev sergey-yaroslavtsev May 13, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not exactly. We find the grid from Z-like type. Then to reuse existing logic we take the first element and find mean step of the grid - this is xScale and yScale. So full grid is not going further.

Copy link
Copy Markdown
Collaborator

@woutdenolf woutdenolf May 14, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What I mean is you are using fast[:n] and slow[::n] which is a single sample. You are not using all the values in fast and slow.

You could do something like

ny = len(fast) // n

fast_scale = np.array([np.median(fast[j::n]) for j in range(n)])

slow_scale = np.array([np.median(slow[i * n : (i + 1) * n]) for i in range(ny)])

But as I said, trying to recognize zig-zag and snake coordinates is very flaky anyway. So this is not the biggest problem of the approach.


if fast_is_A:
X_sh = short_fast
Y_sh = short_slow
else:
X_sh = short_slow
Y_sh = short_fast

return X_sh, Y_sh

def getDimensions(self, nFiles, nScans, shape, index=None):
#somebody may want to overwrite this
"""
Expand Down
Loading