Previous post about spot finding only covered the signal filtering, but a key part of the work is also finding the connected components to identify the complete spot, the total intensity and the centre of mass. This requires folding in an additional analysis step after identifying the signal pixels, to figure out which ones belong to the same spot. In the DIALS implementation of spot finding we use a 4-connected rule i.e. a connection graph which looks like:
010
111
010
So that pixels which are edge-adjacent are described as being in the same spot. As a stream analysis for a given pixel we need only consider the pixel above and to the left of the current one, since the other two will be considered later. If either of these are signal, then the current pixel should be joined to an existing spot. If neither are, we create a new spot. If both are signal and they belong to different spots we need to combine those two spots then add the contribution from the current pixel.
Code from previous post extended to include connected component labelling:
/* Spot structure:
*
* sum of intensity of pixels, intensity * x offset, * y offset, .. etc.
* bounding box, then pointer to parent spot if this gets merged with
* another blob (whereapon n will be reset to 0) - n is the number of
* pixels contributing to the spot.
*
* Sized to be 6 words / 24 bytes per spot.
*/
typedef struct spot {
uint32_t i_sum;
uint32_t ix_sum;
uint32_t iy_sum;
uint16_t x0, x1, y0, y1;
uint16_t parent, n;
} spot;
#define MAX_SPOTS 1024
spot *spots = NULL;
uint32_t nspots = 0;
This is just the structure we will use to store spots - important to note here that not every record will contain a spot on output, since the combination of spots may eliminate counts from a spot but we don’t want to delete the label.
Need to allocate room for the spots we find - we can’t do this dynamically with such a small memory:
spots = (spot *)m_malloc(sizeof(spot) * MAX_SPOTS);
Add a way of returning the number of spots found as an indicator of success:
uint32_t signal_filter_reset(void) {
uint32_t result = 0;
for (int j = 1; j <= nspots; j++) {
if (spots[j].n > 2) result ++;
}
nspots = 0;
return result;
}
Working out how to extract the actual spot positions to µPython would be necessary if we were actually going to use this but happily this is just a proof of concept.
Additional connected component code - this is where we do the heavy lifting. First NULL’ing the 0th spot since we cannot usefully label that one without a lot of extra code complexity.
spot no_spot;
no_spot.i_sum = 0;
no_spot.ix_sum = 0;
no_spot.iy_sum = 0;
no_spot.n = 0;
no_spot.parent = 0;
no_spot.x0 = 0;
no_spot.x1 = 0;
no_spot.y0 = 0;
no_spot.y1 = 0;
spots[0] = no_spot;
And then the actual implementation:
// perform connected component analysis - replaces im[i * width + j]
// pixel value with spot id or 0
if (signal) {
uint16_t above = i > 0 ? im[(i - 1) * width + j] : 0;
uint16_t left = j > 0 ? im[i * width + (j - 1)] : 0;
// FIXME everything to do with x0, y0 etc. will be
// wrong as this is a rolling window - need to keep
// track of a true i here
while (spots[above].parent > 0) {
above = spots[above].parent;
}
while (spots[left].parent > 0) {
left = spots[left].parent;
}
if (above == 0 && left == 0) {
// create a new spot record
nspots++;
spots[nspots].i_sum = p;
spots[nspots].ix_sum = p * j;
spots[nspots].iy_sum = p * i;
spots[nspots].n = 1;
spots[nspots].x0 = j;
spots[nspots].x1 = j;
spots[nspots].y0 = i;
spots[nspots].y1 = i;
spots[nspots].parent = 0;
im[i * width + j] = nspots;
} else if ((above == left) || (above == 0 && left > 0) ||
(above > 0 && left == 0)) {
// merge with correct one
uint16_t keep = MAX(above, left);
spots[keep].i_sum += p;
spots[keep].ix_sum += p * j;
spots[keep].iy_sum += p * i;
spots[keep].n++;
spots[keep].x0 = MIN(j, spots[keep].x0);
spots[keep].x1 = MAX(j, spots[keep].x1);
spots[keep].y0 = MIN(i, spots[keep].y0);
spots[keep].y1 = MAX(i, spots[keep].y1);
im[i * width + j] = keep;
} else if (above != left) {
// deal with the collision and add this spot while we are here
uint16_t keep = MIN(above, left);
uint16_t reject = MAX(above, left);
spot r = spots[reject];
spots[keep].i_sum += r.i_sum;
spots[keep].ix_sum += r.ix_sum;
spots[keep].iy_sum += r.iy_sum;
spots[keep].n += r.n;
spots[keep].x0 = MIN(r.x0, spots[keep].x0);
spots[keep].x1 = MAX(r.x1, spots[keep].x1);
spots[keep].y0 = MIN(r.y0, spots[keep].y0);
spots[keep].y1 = MAX(r.y1, spots[keep].y1);
spots[reject] = no_spot;
spots[reject].parent = keep;
spots[keep].i_sum += p;
spots[keep].ix_sum += p * j;
spots[keep].iy_sum += p * i;
spots[keep].n++;
spots[keep].x0 = MIN(j, spots[keep].x0);
spots[keep].x1 = MAX(j, spots[keep].x1);
spots[keep].y0 = MIN(i, spots[keep].y0);
spots[keep].y1 = MAX(i, spots[keep].y1);
im[i * width + j] = keep;
}
} else {
im[i * width + j] = 0;
}
This is obviously tagged onto the end of the loop in the previous post about identifying signal pixels: we now overwrite the image table with the spot ID (or zero) to keep track of the pixels which belong to this spot.
/sd/frame_00000.raw => 275 pixels / 50 spots in 131.34s
/sd/frame_00001.raw => 693 pixels / 119 spots in 131.55s
/sd/frame_00002.raw => 734 pixels / 124 spots in 131.50s
/sd/frame_00003.raw => 853 pixels / 146 spots in 138.43s
/sd/frame_00004.raw => 754 pixels / 125 spots in 131.43s
/sd/frame_00005.raw => 757 pixels / 133 spots in 131.40s
/sd/frame_00006.raw => 731 pixels / 118 spots in 131.36s
/sd/frame_00007.raw => 766 pixels / 126 spots in 131.31s
/sd/frame_00008.raw => 645 pixels / 111 spots in 138.23s
/sd/frame_00009.raw => 704 pixels / 110 spots in 131.23s
Takes a couple of minutes per image, which won’t win any prizes, but it works and gives the right answers (I also found a bug)
for ((j = 1; j< 11; j++)); do dials.import ../Insulin_6_2.nxs image_range=${j},${j}; dials.find_spots imported.expt algorithm=dispersion ; done| grep "spots found"
=>
50 spots found on 1 image
119 spots found on 1 image
124 spots found on 1 image
146 spots found on 1 image
125 spots found on 1 image
133 spots found on 1 image
118 spots found on 1 image
126 spots found on 1 image
111 spots found on 1 image
110 spots found on 1 image
Qapla’