|
39 | 39 | add_micron_coord_sys, |
40 | 40 | ) |
41 | 41 |
|
42 | | -# def _is_multiscale(element): |
43 | | -# return ( |
44 | | -# hasattr(element, 'keys') |
45 | | -# and callable(element.keys) |
46 | | -# and not isinstance(element, GeoDataFrame) |
47 | | -# ) |
48 | | - |
49 | | -# def get_microns_scales(sdata, element_name): |
50 | | -# el = sdata[element_name] |
51 | | -# if _is_multiscale(el): |
52 | | -# img = sd.get_pyramid_levels(el, n=0) |
53 | | -# img = img.image if hasattr(img, 'image') else img |
54 | | -# else: |
55 | | -# img = el.image if hasattr(el, 'image') else el |
56 | | - |
57 | | -# # Get transforms from the actual image DataArray, not the DataTree |
58 | | -# el_transforms = get_transformation(img, get_all=True) |
59 | | -# microns_tf = el_transforms.get('microns', None) |
60 | | -# if microns_tf is None: |
61 | | -# ps = sdata['table'].uns['section_metadata']['pixel_size'] |
62 | | -# microns_tf = Scale([ps, ps], axes=('x', 'y')) |
63 | | -# set_transformation(img, microns_tf, to_coordinate_system='microns') |
64 | | -# if len(microns_tf.scale) >= 2: |
65 | | -# x_y_axes = ('x', 'y') |
66 | | -# x_y_tf = [microns_tf.axes.index(axis) for axis in x_y_axes if axis in microns_tf.axes] |
67 | | -# microns_tf = Scale([microns_tf.scale[i] for i in x_y_tf], x_y_axes) |
68 | | -# return microns_tf |
69 | | - |
70 | 42 | def get_zstack_zarr(alignment_folder, paths, zstack_fov_size=None): |
71 | 43 | zstack_folder = alignment_folder / 'zstacks' |
72 | 44 | zstack_folder.mkdir(parents=True, exist_ok=True) |
@@ -322,265 +294,6 @@ def make_element_3d( |
322 | 294 |
|
323 | 295 | return sdata |
324 | 296 |
|
325 | | - |
326 | | - |
327 | | -# ####### Functions to handle landmarks ######### |
328 | | -# def get_biwarp_params(bigwarp_json_path): |
329 | | -# import json |
330 | | -# with open(bigwarp_json_path, 'r') as f: |
331 | | -# bw_json = json.load(f) |
332 | | -# for source, source_data in bw_json['Sources'].items(): |
333 | | -# if source_data['isMoving']: |
334 | | -# if 'confocal' in source_data['uri'].lower(): |
335 | | -# continue |
336 | | -# moving_image_path = Path(source_data['uri']) |
337 | | -# if moving_image_path.name=='?': |
338 | | -# moving_image_path = moving_image_path.parent |
339 | | -# else: |
340 | | -# reference_image_path = Path(source_data['uri']) |
341 | | -# if reference_image_path.name=='?': |
342 | | -# reference_image_path = reference_image_path.parent |
343 | | -# if 'zstack' in moving_image_path.name.lower() or 'z_stack' in moving_image_path.name.lower(): |
344 | | -# moving_image = 'czstack' |
345 | | -# reference_image = 'xenium_section' |
346 | | -# transform_type = bw_json['Transform']['type'] |
347 | | -# bigwarp_params = {'moving_image': moving_image, |
348 | | -# 'reference_image': reference_image, |
349 | | -# 'moving_image_path': moving_image_path, |
350 | | -# 'reference_image_path': reference_image_path, |
351 | | -# 'transform_type': transform_type} |
352 | | -# return bigwarp_params |
353 | | - |
354 | | -# def invert_xenium_y_landmarks(landmarks, landmarked_image_path): |
355 | | -# with tifffile.TiffFile(landmarked_image_path) as tif: |
356 | | -# landmarked_image_shape = tif.pages[0].shape |
357 | | -# full_y_size = landmarked_image_shape[0] |
358 | | -# landmarks['xenium_y'] = full_y_size - landmarks['xenium_y'] |
359 | | -# return landmarks |
360 | | - |
361 | | -# def remove_landmark_buffer(landmarks, czstack_buffer=None, xenium_buffer=None): |
362 | | -# if czstack_buffer is not None: |
363 | | -# landmarks['czstack_y'] = landmarks['czstack_y'] - czstack_buffer.get('y', 0) |
364 | | -# landmarks['czstack_x'] = landmarks['czstack_x'] - czstack_buffer.get('x', 0) |
365 | | -# landmarks['czstack_z'] = landmarks['czstack_z'] - czstack_buffer.get('z', 0) |
366 | | -# if xenium_buffer is not None: |
367 | | -# landmarks['xenium_y'] = landmarks['xenium_y'] - xenium_buffer.get('y', 0) |
368 | | -# landmarks['xenium_x'] = landmarks['xenium_x'] - xenium_buffer.get('x', 0) |
369 | | -# landmarks['xenium_z'] = landmarks['xenium_z'] - xenium_buffer.get('z', 0) |
370 | | -# return landmarks |
371 | | - |
372 | | -# def get_section_landmarks(landmarks_path, dims_order=['x','y','z'], bigwarp_project_path=None, moving_img=None): |
373 | | -# if bigwarp_project_path is not None: |
374 | | -# bigwarp_params = get_biwarp_params(bigwarp_project_path) |
375 | | -# elif moving_img is not None: |
376 | | -# bigwarp_params = {'moving_image': moving_img} |
377 | | -# else: |
378 | | -# print("No BigWarp project path or moving image specified - assuming czstack was moving image") |
379 | | -# bigwarp_params = {'moving_image': 'czstack'} |
380 | | -# print(f"Loading landmarks from: {landmarks_path}") |
381 | | -# landmarks = pd.read_csv(landmarks_path, header=None) |
382 | | -# if bigwarp_params['moving_image']=='czstack': |
383 | | -# # Flatten the lists using + operator to concatenate them |
384 | | -# landmarks.columns = ['landmark_name', 'active'] + [f'czstack_{dim}' for dim in dims_order] + [f'xenium_{dim}' for dim in dims_order] |
385 | | -# else: |
386 | | -# landmarks.columns = ['landmark_name', 'active'] + [f'xenium_{dim}' for dim in dims_order] + [f'czstack_{dim}' for dim in dims_order] |
387 | | - |
388 | | -# return landmarks, bigwarp_params |
389 | | - |
390 | | -# def get_landmarked_image_props(landmarked_image_path, sdata, landmarks, |
391 | | -# section_n, invert_y=False): |
392 | | -# with tifffile.TiffFile(landmarked_image_path) as tif: |
393 | | -# landmarked_image_shape = tif.pages[0].shape |
394 | | -# print(f"Landmarked image shape: {landmarked_image_shape}") |
395 | | - |
396 | | -# bbox_xmin = 0 |
397 | | -# bbox_ymin = 0 |
398 | | - |
399 | | -# cell_labels = sd.get_pyramid_levels(sdata['cell_labels'], n=0) |
400 | | -# section_shape_y, section_shape_x = cell_labels.data.shape |
401 | | - |
402 | | -# bbox_dict = sdata['table'].uns.get('sections_bboxes', None) |
403 | | - |
404 | | -# if bbox_dict is not None and str(section_n) in bbox_dict: |
405 | | -# # Paired: invert within image space FIRST, then scale to full slide + offset |
406 | | -# if invert_y: |
407 | | -# landmarks['xenium_y'] = landmarked_image_shape[0] - landmarks['xenium_y'] |
408 | | - |
409 | | -# section_bbox = bbox_dict[str(section_n)] |
410 | | -# bbox_xmin = section_bbox['x_min'] |
411 | | -# bbox_ymin = section_bbox['y_min'] |
412 | | -# full_slide_shape_y = np.max([bbox['y_max'] for bbox in bbox_dict.values()]) |
413 | | -# full_slide_shape_x = np.max([bbox['x_max'] for bbox in bbox_dict.values()]) |
414 | | -# scale_factor_y = full_slide_shape_y / landmarked_image_shape[0] |
415 | | -# scale_factor_x = full_slide_shape_x / landmarked_image_shape[1] |
416 | | -# print(f"Paired section — full slide shape: {[full_slide_shape_y, full_slide_shape_x]}") |
417 | | -# print(f"Image downsampled by: (yx) [{scale_factor_y:.4f}, {scale_factor_x:.4f}]") |
418 | | -# print(f"Section bbox offset: {section_bbox}") |
419 | | -# else: |
420 | | -# # Standalone: scale to full section first, then invert within section space |
421 | | -# scale_factor_y = section_shape_y / landmarked_image_shape[0] |
422 | | -# scale_factor_x = section_shape_x / landmarked_image_shape[1] |
423 | | -# print(f"Standalone section — section array shape: {[section_shape_y, section_shape_x]}") |
424 | | -# print(f"Image downsampled by: (yx) [{scale_factor_y:.4f}, {scale_factor_x:.4f}]") |
425 | | - |
426 | | -# landmarks['xenium_x'] = landmarks['xenium_x'] * scale_factor_x |
427 | | -# landmarks['xenium_y'] = landmarks['xenium_y'] * scale_factor_y |
428 | | -# landmarks['xenium_x'] = landmarks['xenium_x'] - bbox_xmin |
429 | | -# landmarks['xenium_y'] = landmarks['xenium_y'] - bbox_ymin |
430 | | - |
431 | | -# if bbox_dict is None or str(section_n) not in bbox_dict: |
432 | | -# # Standalone: invert after scaling, using full section height |
433 | | -# if invert_y: |
434 | | -# landmarks['xenium_y'] = section_shape_y - landmarks['xenium_y'] |
435 | | - |
436 | | -# return landmarks |
437 | | - |
438 | | -# def format_landmarks(sdata, |
439 | | -# landmarks_path, |
440 | | -# section_n, |
441 | | -# bigwarp_project_paths=None, |
442 | | -# moving_img=None, |
443 | | -# czstack_buffer=None, |
444 | | -# invert_lm_y=False, |
445 | | -# fix_cropped_landmarks=False, |
446 | | -# landmarked_image_path=None, |
447 | | -# dims_order=['x','y','z']): |
448 | | - |
449 | | -# landmarks, bigwarp_params = get_section_landmarks( |
450 | | -# landmarks_path=landmarks_path, |
451 | | -# bigwarp_project_path=bigwarp_project_paths, |
452 | | -# moving_img=moving_img, |
453 | | -# dims_order=dims_order) |
454 | | - |
455 | | -# if czstack_buffer is not None: |
456 | | -# landmarks = remove_landmark_buffer(landmarks, |
457 | | -# czstack_buffer=czstack_buffer) |
458 | | - |
459 | | -# if fix_cropped_landmarks and landmarked_image_path is not None: |
460 | | -# # invert_y is handled inside get_landmarked_image_props at the right point |
461 | | -# landmarks = get_landmarked_image_props( |
462 | | -# landmarked_image_path, |
463 | | -# sdata, |
464 | | -# landmarks, |
465 | | -# section_n, |
466 | | -# invert_y=invert_lm_y) |
467 | | -# elif invert_lm_y and landmarked_image_path is not None: |
468 | | -# # sections where fix_cropped_landmarks=False |
469 | | -# landmarks = invert_xenium_y_landmarks(landmarks, landmarked_image_path) |
470 | | - |
471 | | -# # Adjust landmarks resolutions |
472 | | -# full_scale_pixel_size = sdata['table'].uns['section_metadata']['pixel_size'] |
473 | | -# landmarks = landmarks.rename(columns={'xenium_x': 'x', 'xenium_y': 'y', 'xenium_z': 'z'}) |
474 | | -# landmarks = PointsModel.parse(landmarks) |
475 | | -# set_transformation(landmarks, Identity(), to_coordinate_system='global') |
476 | | -# set_transformation( |
477 | | -# landmarks, |
478 | | -# Scale([full_scale_pixel_size, full_scale_pixel_size], axes=('x', 'y')), |
479 | | -# to_coordinate_system='microns' |
480 | | -# ) |
481 | | -# return landmarks |
482 | | - |
483 | | -# ####### Functions to find transforms ########## |
484 | | -# def get_affine_from_landmarks_flat(moving_coords, ref_coords): |
485 | | -# moving_2d = moving_coords[:, :2] |
486 | | -# ref_2d = ref_coords[:, :2] |
487 | | -# n = moving_2d.shape[0] |
488 | | -# A_xy = np.hstack([moving_2d, np.ones((n, 1))]) |
489 | | -# result, _, _, _ = np.linalg.lstsq(A_xy, ref_2d, rcond=None) |
490 | | - |
491 | | -# # Z offset: mean difference between ref and moving Z coords |
492 | | -# # mat[2,2] = 1.0 keeps matrix invertible (required by Napari) |
493 | | -# z_offset = float(np.mean(ref_coords[:, 2] - moving_coords[:, 2])) |
494 | | - |
495 | | -# mat = np.eye(4) |
496 | | -# mat[0, 0] = result[0, 0]; mat[0, 1] = result[1, 0]; mat[0, 3] = result[2, 0] |
497 | | -# mat[1, 0] = result[0, 1]; mat[1, 1] = result[1, 1]; mat[1, 3] = result[2, 1] |
498 | | -# mat[2, 2] = 1.0 # MUST be 1.0 - keeps matrix invertible for Napari rendering |
499 | | -# mat[2, 3] = z_offset |
500 | | - |
501 | | -# return Affine(mat, input_axes=('x', 'y', 'z'), output_axes=('x', 'y', 'z')) |
502 | | - |
503 | | -# def tilt_affines(moving_pts, fixed_pts, flat_affine): |
504 | | -# # moving_pts = Xenium (2D: x, y, 0) |
505 | | -# # fixed_pts = CZStack (3D: x, y, z) |
506 | | - |
507 | | -# # Fit Z_stack = a*X_xe + b*Y_xe + c |
508 | | -# A = np.column_stack([moving_pts[:, 0], moving_pts[:, 1], np.ones(len(moving_pts))]) |
509 | | -# coeffs, _, _, _ = lstsq(A, fixed_pts[:, 2]) |
510 | | -# a, b, c = coeffs |
511 | | - |
512 | | -# # Build the 3D matrix using the 2D XY results |
513 | | -# tilt_mat = np.array(flat_affine.matrix).copy() |
514 | | -# tilt_mat[2, 0] = a # Effect of Xenium X on Stack Z |
515 | | -# tilt_mat[2, 1] = b # Effect of Xenium Y on Stack Z |
516 | | -# tilt_mat[2, 2] = 1.0 # Keeps matrix invertible |
517 | | -# tilt_mat[2, 3] = c # The base Z-slice offset |
518 | | - |
519 | | -# return Affine(tilt_mat, input_axes=('x', 'y', 'z'), output_axes=('x', 'y', 'z')) |
520 | | - |
521 | | -# def _is_multiscale(element): |
522 | | -# return ( |
523 | | -# hasattr(element, 'keys') |
524 | | -# and callable(element.keys) |
525 | | -# and not isinstance(element, GeoDataFrame) |
526 | | -# ) |
527 | | - |
528 | | -# def get_microns_scales(sdata, element_name): |
529 | | -# el = sdata[element_name] |
530 | | -# if _is_multiscale(el): |
531 | | -# img = sd.get_pyramid_levels(el, n=0) |
532 | | -# img = img.image if hasattr(img, 'image') else img |
533 | | -# else: |
534 | | -# img = el.image if hasattr(el, 'image') else el |
535 | | - |
536 | | -# # Get transforms from the actual image DataArray, not the DataTree |
537 | | -# el_transforms = get_transformation(img, get_all=True) |
538 | | -# microns_tf = el_transforms.get('microns', None) |
539 | | -# if microns_tf is None: |
540 | | -# ps = sdata['table'].uns['section_metadata']['pixel_size'] |
541 | | -# microns_tf = Scale([ps, ps], axes=('x', 'y')) |
542 | | -# set_transformation(img, microns_tf, to_coordinate_system='microns') |
543 | | -# if len(microns_tf.scale) >= 2: |
544 | | -# x_y_axes = ('x', 'y') |
545 | | -# x_y_tf = [microns_tf.axes.index(axis) for axis in x_y_axes if axis in microns_tf.axes] |
546 | | -# microns_tf = Scale([microns_tf.scale[i] for i in x_y_tf], x_y_axes) |
547 | | -# return microns_tf |
548 | | - |
549 | | -# def add_affine_to_element(element, affine_tf, coord_sys_name, |
550 | | -# microns_tf=None, microns_tf_position='before', |
551 | | -# overwrite_existing=False): |
552 | | -# """ |
553 | | -# microns_tf_position: 'before' puts microns_tf before affine_tf (section px → µm → czstack px) |
554 | | -# 'after' puts microns_tf after affine_tf (section px → czstack px → czstack µm) |
555 | | -# """ |
556 | | -# def _build_tf(tf_to_fullres): |
557 | | -# if microns_tf is not None: |
558 | | -# if microns_tf_position == 'after': |
559 | | -# return Sequence([tf_to_fullres, affine_tf, microns_tf]) |
560 | | -# else: |
561 | | -# return Sequence([tf_to_fullres, microns_tf, affine_tf]) |
562 | | -# return Sequence([tf_to_fullres, affine_tf]) |
563 | | - |
564 | | -# if _is_multiscale(element): |
565 | | -# for n_l, level in enumerate(element.keys()): |
566 | | -# element_obj = sd.get_pyramid_levels(element, n=n_l) |
567 | | -# tf_to_fullres = get_transformation(element_obj, to_coordinate_system='global') |
568 | | -# if isinstance(tf_to_fullres, Identity): |
569 | | -# tf_to_fullres = Scale([1.0, 1.0, 1.0], axes=('x', 'y', 'z')) |
570 | | -# if coord_sys_name in element_obj.attrs['transform'] and not overwrite_existing: |
571 | | -# continue |
572 | | -# set_transformation(element_obj, _build_tf(tf_to_fullres), |
573 | | -# to_coordinate_system=coord_sys_name) |
574 | | -# else: |
575 | | -# element_obj = element |
576 | | -# tf_to_fullres = get_transformation(element_obj, to_coordinate_system='global') |
577 | | -# if isinstance(tf_to_fullres, Identity): |
578 | | -# tf_to_fullres = Scale([1.0, 1.0, 1.0], axes=('x', 'y', 'z')) |
579 | | -# if coord_sys_name in element_obj.attrs['transform'] and not overwrite_existing: |
580 | | -# return |
581 | | -# set_transformation(element_obj, _build_tf(tf_to_fullres), |
582 | | -# to_coordinate_system=coord_sys_name) |
583 | | - |
584 | 297 | ####### QC alignment ######### |
585 | 298 | def get_section_z_stats(affines_dict, landmarks_dict, czstack_shape_yx=(512, 512)): |
586 | 299 | y_max, x_max = czstack_shape_yx |
|
0 commit comments