Semi-automatic stitching of filamentous structures in image stacks from serial-section electron tomography

We present a software-assisted workflow for the alignment and matching of filamentous structures across a 3D stack of serial images. This is achieved by combining automatic methods, visual validation, and interactive correction. After an initial alignment, the user can continuously improve the result by interactively correcting landmarks or matches of filaments. Supported by a visual quality assessment of regions that have been already inspected, this allows a trade-off between quality and manual labor. The software tool was developed to investigate cell division by quantitative 3D analysis of microtubules (MTs) in both mitotic and meiotic spindles. For this, each spindle is cut into a series of semi-thick physical sections, of which electron tomograms are acquired. The serial tomograms are then stitched and non-rigidly aligned to allow tracing and connecting of MTs across tomogram boundaries. In practice, automatic stitching alone provides only an incomplete solution, because large physical distortions and a low signal-to-noise ratio often cause experimental difficulties. To derive 3D models of spindles despite the problems related to sample preparation and subsequent data collection, semi-automatic validation and correction is required to remove stitching mistakes. However, due to the large number of MTs in spindles (up to 30k) and their resulting dense spatial arrangement, a naive inspection of each MT is too time consuming. Furthermore, an interactive visualization of the full image stack is hampered by the size of the data (up to 100 GB). Here, we present a specialized, interactive, semi-automatic solution that considers all requirements for large-scale stitching of filamentous structures in serial-section image stacks. The key to our solution is a careful design of the visualization and interaction tools for each processing step to guarantee real-time response, and an optimized workflow that efficiently guides the user through datasets. Author summary Electron tomography of biological samples is used for a 3D reconstruction of filamentous structures, such as microtubules (MTs) in mitotic and meiotic spindles. Large-scale electron tomography can be applied to increase the reconstructed volume for the visualization of full spindles. For this, each spindle is cut into a series of semi-thick physical sections, of which electron tomograms are acquired. The serial tomograms are then stitched and non-rigidly aligned to allow tracing and connecting of MTs across tomogram boundaries. Previously, we presented fully automatic approaches for this 3D reconstruction pipeline. However, large volumes often suffer from imperfections (i.e. physical distortions) caused during sectioning and imaging, making it difficult to apply fully automatic approaches for matching and stitching of numerous tomograms. Therefore, we developed an interactive, semi-automatic solution that considers all requirements for large-scale stitching of microtubules in serial-section image stacks. We achieved this by combining automatic methods, visual validation and interactive error correction, thus allowing the user to continuously improve the result by interactively correcting landmarks or matches of filaments. We present large-scale reconstructions of spindles in which the automatic workflow failed and where different steps of manual corrections were needed. Our approach is also applicable to other biological samples showing 3D distributions of MTs in a number of different cellular contexts.

174 The serial tomogram alignment is computed automatically with a linear 175 transformation. The networks that can be handled are sparse and small compared 176 to the microtubules data sets that are the focus here.
177 Our goal was to stitch microtubules across a whole stack of serial-section electron 178 tomograms in a semi-automatic way allowing for corrections during the stitching 179 process. A number of requirements were identified for such an approach. The user 180 should be able to influence all steps of the semi-automatic stitching process. In 181 addition, the pipeline should be flexible such that the user can change the alignment 182 or matching at any time. In the worst case, the user should be able to perform the 183 whole stitching manually. In the best case, the user should be able to apply the 184 stitching completely automatically. For cases of medium quality, minor manual 185 interactions should support the automatic algorithms to rapidly increase the quality 9 186 of the reconstruction. The whole pipeline should be supported by visualizations to 187 verify the current state and offer options to interactively correct errors. Finally, the 188 application should be able to interactively handle all relevant data sizes, which could 189 be up to 30 tomograms with 5 GB each.
190 In this paper, we present a software-assisted workflow to address these 191 requirements. Two related problems need to be solved: (1) the alignment problem, 192 where the serial tomograms are transformed such that they fit together; and (2) the 193 matching problem, where the filament structures are combinatorially connected 194 between neighboring tomograms. In summary, we make the following key 195 contributions:

196
 Design of a user-controlled workflow (Fig 2) (Fig 1B), 1 ,…, 240 we expect that each tomogram has been pre-processed, i.e distortions along the 241 thin tomogram side, which are usually the z-direction, have been removed in a 242 flattening pre-processing step, followed by a cropping step. Furthermore, the 243 filaments have been segmented and verified for each tomogram. We will denote the 244 set of filaments of tomogram by , where each filament is a piecewise ∈ 245 linear curve with and . We further assume that the stack is 246 always built in z-direction. We call a slice of the tomogram that lies in the x-y-plane 247 a z-slice w.r.t. a certain z-value.
248 In addition to the tomograms and their corresponding filament sets, the 'Serial 249 Section Stack' also stores the stitching state for each of the interfaces between -1 250 neighboring tomograms. A stitching state contains sets of landmarks on the 251 tomogram boundaries and sets of matchings representing corresponding filament 252 ends.

253
Serial section aligner 254 The 'Serial Section Aligner' module provides two modes. The alignment mode 255 allows the user to manually add, modify, and remove landmarks. The matching 256 mode, on the other hand, enables the automatic computation of a matching in 257 combination with manual editing and verification. In order to visualize the current 258 state of the matching and alignment and to allow for manual interaction, the module 259 uses a four-view visualization (Fig 3). In addition, the module has a graphical user 260 interface (GUI) to adjust the module parameters (Fig 2A). The first part of the GUI 261 is independent of the specific mode. It allows selecting the currently investigated   14 301 Finally, the bottom right view shows an abstract illustration of the whole stack and 302 the selected interface with the pair of tomograms. The user can change the interface 303 with a mouse click onto the illustration. Again, the same colors are used, yellow for 304 the bottom and blue for the top tomogram. Next to the stack, a vertical histogram is 305 shown. Each bin of the histogram shows the density of the unmatched filament ends 306 for 25% of a tomogram extent in z-direction. At the beginning, the histogram usually 307 shows higher densities at the interface boundaries, because the filaments are not 308 yet matched. During the stitching process these differences decrease as more and 309 more filament ends are matched. The plot helps the user to evaluate the overall 310 quality but also to identify the most critical interfaces that need further investigation. 313 In addition to manual landmarks, automatic landmarks (yellow) can be generated 314 from matched filaments. This approach is described in the following section. The 315 user can also remove and manipulate these landmarks. However, if the automatic 316 landmark placement is activated, these landmarks will be updated by changes in 317 the matching.
318 Matching mode. In this mode, the user defines the matching between the filaments 319 of the two tomograms. This is supported by four views and the possibility to 320 automatically compute the matching in combination with manual editing (Fig 4).  (Fig 2A), the ends of matched filaments 357 will be used to generate automatic landmarks for the alignment of the tomograms.

Recommended workflow
413 For applying the 'Serial Section Aligner', we recommend the following overall 414 workflow (Fig 2B). After creating the 'Serial Section Stack' object and detecting the 415 correct tomogram order, the section interfaces will be processed one after another.  (Fig 5B,D). Finally, the same procedure is followed as described 440 above for a direct successful automatic matching.

448
Technical aspects 449 Stack data structure. As described above (see section Serial Section Stack), the 450 serial section data is handled in a stack data structure. If is the number of sections, 451 then the stack stores tomograms with their corresponding filament data sets .
452 The data sets are stored in the file system. The stack data structure loads 453 tomograms and filaments on demand and hands them over to the 'Serial Section 454 Aligner' module. During alignment and matching, only two tomograms are loaded 455 into RAM at the same time.
456 For each of the section interfaces between consecutive tomograms, a set of -1 457 landmarks and a matching is stored. As briefly described above, landmarks always 458 occur as pairs of 2D points in x-y-space, one point for each of the neighboring 459 tomograms. The points mark corresponding locations on the two sides of an 460 interface and define an x-y-transformation that warps one tomogram into the space 461 of the neighboring one such that the tomograms and filaments are aligned. Because 462 the tomograms are thin with a flat boundary in z-direction, after the pre-processing 463 step, warping is only defined as an x-y-transformation, ignoring z. In addition to the 464 pair of points, each landmark has a type attribute that stores the source of the 465 landmark. The source can be manual or automatic. As described earlier, these types 466 will be colored in red and light yellow, respectively. 512 For all points within a grid cell, bilinear interpolation of the warped cell is applied to 513 compute the warped position (Fig 6). This accelerates the rendering of the warped 514 slice, so that it can be interactively updated during landmark or slice changes.
515 Additionally, it creates a kind of level of detail: the correctness of the warping 516 increases the closer the user zooms into the slice. Since the warping affects only 517 the x-y-space of the tomogram, it is sufficient to use a 1D grid for side views.

518
We can further increase the rendering performance for the warped slices using an 519 approximation of the moving least-squares approach. In its original implementation, 520 all landmarks are used to transform a point. However, since the influence of a 521 landmark decreases with increasing distance, it is sufficient to detect the closest k 522 landmarks to transform a point. To do this in an efficient way, all landmarks are 523 stored in a grid data structure.
524 Landmark rendering. Each landmark is visualized in the alignment mode as a 525 circle on the z-slice (Fig 3A,B). The circles have a constant size in screen space 526 with a scaling parameter. This seems to be superior to object-space radii, because 535 This allows the user to see the image feature that is related to the landmark. 547 Similarly, for each cylinder, two 3D vertices are created for the start and end points, 548 respectively. Then in the geometry shader, a quad is created that encloses the 549 cylinder after projection. Finally, ray casting is used in the fragment shader to find 550 the intersection of the view ray with the cylinder.
551 To quickly update the filament visualization of the warped tomogram in the overlay 552 view, all spheres and cylinders that intersect the slice are detected in a first step.
553 Then, only these spheres and cylinders are warped and visualized in a second step.
554 Both, filaments and landmarks, are rendered first in a frame buffer. During the 555 rendering on the slice, a shader is used to apply anti-aliasing to the boundaries. To 556 do so, for each fragment at the boundary, the tangent line is computed based on an 557 analysis of the neighboring fragments in the four main directions. The line is used 558 to compute a blending value for the fragment. Furthermore, for close-up views on 559 the filaments or landmarks, the line is used to create a pixel-wide anti-aliased black 560 silhouette. This procedure drastically increases the visual quality and helps to 561 distinguish between slice data and filaments or landmarks.
562 Filament end state map. In matching mode, a map depicts which filament ends are 563 already confirmed, which ends are automatically matched but not confirmed, and 564 which ends are unmatched and not confirmed (Fig 7). This map is generated and 565 updated in the following way.
25 566 First, a discrete 2D scalar field is created that has the same boundary in the x-and 567 y-directions as the bottom tomogram. Then for each sample position of the scalar 568 field, the closest distance to all filament ends in the matching region of the bottom 569 section is computed. The identifier of the closest filament end is stored at this 570 position in the scalar field. Hence, the scalar field is a discrete representation of the 571 Voronoi diagram of the filament ends in the matching region. Note that this map 572 needs to be updated only when the user changes the matching region.
573 Then the z-slice of the bottom section is rendered as described above. However, 574 during the detection of the intensity value in the tomogram for a pixel, additionally 575 the closest filament end is determined from the pre-computed Voronoi diagram.
576 Then the state of the end is used to modulate the brightness of the intensity value,  586 In the following paragraphs, three challenging datasets are described in detail. For 587 the first two, experimental problems had to be tackled, such as damaged sections 588 or poor image quality. The third example represents the largest spindle that has 589 been reconstructed so far. We used it as benchmark for our system.