working on Stack alignment

render-ws-java-client/src/main/java/org/janelia/render/client/multisem/RecapKensAlignment.java

@@ -1,29 +1,39 @@
 package org.janelia.render.client.multisem;
 
+import java.awt.Rectangle;
 import java.io.File;
 import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.HashMap;
 
-import loci.common.DebugTools;
+import mpicbg.models.AffineModel2D;
 import mpicbg.models.InvertibleBoundable;
+import mpicbg.models.RigidModel2D;
 import mpicbg.models.TranslationModel2D;
 import mpicbg.stitching.ImageCollectionElement;
 import mpicbg.stitching.fusion.Fusion;
+import mpicbg.trakem2.transform.CoordinateTransform;
+import mpicbg.trakem2.transform.CoordinateTransformList;
 
 public class RecapKensAlignment
 {
 	public static String basePath = "/nrs/hess/from_mdas/Hayworth/Wafer53/";
 	public static String stitchedSlabs = "/scan_corrected_equalized_target_dir/stitched_stacks";
+	public static String rigidSlabs = "/scan_corrected_equalized_target_dir/stitched_stacks_substrate_replaced";
 
 	public static class TransformedImage
 	{
 		public String fileName;
-		public int slab;
-		public int z;
+		public int slab, z, width, height;
 		public ArrayList< InvertibleBoundable > models = new ArrayList<>();
 	}
 
+	public static class TransformedZLayer
+	{
+		public ArrayList< TransformedImage > transformedImages;
+		public int width, height;
+	}
+
 	/**
 	 * This already ignores slices that were dropped at a later point
 	 */
@@ -45,7 +55,7 @@ public static ArrayList< Integer > numSlices( final int slab )
 
 	public static void reconstruct( final int slab )
 	{
-		final HashMap< Integer, ArrayList< TransformedImage > > models = new HashMap<>();
+		final HashMap< Integer, TransformedZLayer > models = new HashMap<>();
 
 		// find out numSlices and indices
 		final ArrayList< Integer > slices = numSlices( slab );
@@ -59,7 +69,7 @@ public static void reconstruct( final int slab )
 
 
 		// Stitching
-		System.out.println( "STITCHING" );
+		System.out.println( "\nSTITCHING" );
 
 		for ( int zIndex = 0; zIndex < numSlices; ++zIndex )
 		{
@@ -89,14 +99,17 @@ public static void reconstruct( final int slab )
 				imgSizes[ i ][ 0 ] = 1996;
 				imgSizes[ i ][ 1 ] = 1748;
 				stitchingModels.add( t );
-				++i;
 
 				final TransformedImage tI = new TransformedImage();
 				tI.fileName = e.getFile().getAbsolutePath();
 				tI.slab = slab;
 				tI.z = z;
+				tI.width = imgSizes[ i ][ 0 ];
+				tI.height = imgSizes[ i ][ 1 ];
 				tI.models.add( t );
 
+				++i;
+
 				modelsZ.add( tI );
 				//System.out.println( e.getFile().getAbsolutePath() );
 				//System.out.println( Arrays.toString( e.getOffset() ) + ", " + e.getDimensionality() );
@@ -116,9 +129,93 @@ public static void reconstruct( final int slab )
 			bbStitching.set( -offset[ 0 ], -offset[ 1 ] );
 
 			modelsZ.forEach( m -> m.models.add( bbStitching ) );
-			models.put( z, modelsZ );
+
+			final TransformedZLayer tzl = new TransformedZLayer();
+			tzl.transformedImages = modelsZ;
+			tzl.width = size[ 0 ];
+			tzl.height = size[ 1 ];
+
+			models.put( z, tzl );
 		}
 
+
+		// Rigid registration
+		System.out.println( "\nRIGID REGISTRATION" );
+
+		// this directory contains all XML's of FIJI Register_Virtual_Stack_MT
+		final File dir = new File( basePath, String.format( "%s/%03d_Enhanced_Transforms", rigidSlabs, slab ) );
+
+		System.out.println( dir.getAbsolutePath() );
+		System.out.println( Arrays.toString( dir.list( (d,fn) -> fn.toLowerCase().endsWith(".xml" ) ) ) );
+
+		// Common bounds to create common frame for all images
+		final Rectangle commonBounds = new Rectangle(0, 0, models.get( slices.get( 0 ) ).width, models.get( slices.get( 0 ) ).height );
+
+		for ( int zIndex = 0; zIndex < numSlices; ++zIndex )
+		{
+			final int z = slices.get( zIndex );
+			final TransformedZLayer tzl = models.get( z );
+
+			final String fn = new File( dir, String.format( "StichedImage_scan_%03d.xml", z ) ).getAbsolutePath();
+			System.out.println( "Processing z=" + z + " (" + fn + ")" );
+
+			final CoordinateTransformList<CoordinateTransform> transforms = RecapKensAlignmentTools.readCoordinateTransform( fn );
+
+			// the second transformation (translation, only exists from the 2nd z layer on) is the bounding offset of the previous z layer
+			// I think the idea is that the offset is ignored for each individual z-plane, and rather the next one is moved accordingly
+			// I guess Albert did not care that a bit of the data is potentially cut off(?)
+			final int numTransforms = transforms.getList( new ArrayList<>() ).size();
+			//System.out.println( "Number of transforms: " + numTransforms );
+
+			// update the transformation list of all images in this z-plane
+			for ( int t = 0; t < numTransforms; ++t )
+			{
+				final CoordinateTransform m = transforms.get( t );
+
+				if ( mpicbg.trakem2.transform.RigidModel2D.class.isInstance( m ) )
+				{
+					//System.out.println( "rigid: " + m );
+					models.get( z ).transformedImages.forEach( ti -> ti.models.add( (RigidModel2D)m ) );
+				}
+				else if ( mpicbg.trakem2.transform.TranslationModel2D.class.isInstance( m ) )
+				{
+					//System.out.println( "translation: " + m );
+					models.get( z ).transformedImages.forEach( ti -> ti.models.add( (TranslationModel2D)m ) );
+				}
+				else if ( mpicbg.trakem2.transform.AffineModel2D.class.isInstance( m ) )
+				{
+					//System.out.println( "affine: " + m );
+					models.get( z ).transformedImages.forEach( ti -> ti.models.add( (AffineModel2D)m ) );
+				}
+				else
+				{
+					throw new RuntimeException( "Don't know how to process model: " + m.getClass().getName() );
+				}
+			}
+
+			final Rectangle boundingbox = RecapKensAlignmentTools.getBoundingBox( tzl.width, tzl.height, transforms );
+			System.out.println( "Bounding box = " + boundingbox);
+
+			// Update common bounds
+			int min_x = commonBounds.x;
+			int min_y = commonBounds.y;
+			int max_x = commonBounds.x + commonBounds.width;
+			int max_y = commonBounds.y + commonBounds.height;
+
+			if(boundingbox.x < commonBounds.x)
+				min_x = boundingbox.x;
+			if(boundingbox.y < commonBounds.y)
+				min_y = boundingbox.y;
+			if(boundingbox.x + boundingbox.width > max_x)
+				max_x = boundingbox.x + boundingbox.width;
+			if(boundingbox.y + boundingbox.height > max_y)
+				max_y = boundingbox.y + boundingbox.height;
+
+			commonBounds.x = min_x;
+			commonBounds.y = min_y;
+			commonBounds.width = max_x - min_x;
+			commonBounds.height = max_y - min_y;
+		}
 	}
 
 	public static void main( String[] args )

render-ws-java-client/src/main/java/org/janelia/render/client/multisem/RecapKensAlignmentTools.java

@@ -1,17 +1,27 @@
 package org.janelia.render.client.multisem;
 
+import java.awt.Rectangle;
 import java.io.BufferedReader;
 import java.io.File;
+import java.io.FileNotFoundException;
+import java.io.FileReader;
 import java.io.IOException;
 import java.util.ArrayList;
 import java.util.HashMap;
 import java.util.Map;
+import java.util.concurrent.Callable;
 
 import ij.ImagePlus;
+import ij.io.FileSaver;
+import ij.process.ImageProcessor;
 import mpicbg.models.TranslationModel2D;
 import mpicbg.models.TranslationModel3D;
 import mpicbg.stitching.ImageCollectionElement;
 import mpicbg.stitching.TextFileAccess;
+import mpicbg.trakem2.transform.CoordinateTransform;
+import mpicbg.trakem2.transform.CoordinateTransformList;
+import mpicbg.trakem2.transform.TransformMesh;
+import mpicbg.trakem2.transform.TransformMeshMapping;
 import stitching.utils.Log;
 
 public class RecapKensAlignmentTools
@@ -170,4 +180,91 @@ public static ArrayList< ImageCollectionElement > getLayoutFromFile( final Strin
 		return elements;
 	}
 
+	// from: register_virtual_stack_slices/src/main/java/register_virtual_stack/Transform_Virtual_Stack_MT.java
+	/**
+	 * Read coordinate transform from file (generated in Register_Virtual_Stack)
+	 *
+	 * @param filename  complete file name (including path)
+	 * @return true if the coordinate transform was properly read, false otherwise.
+	 */
+	public static CoordinateTransformList<CoordinateTransform> readCoordinateTransform( String filename )
+	{
+		final CoordinateTransformList<CoordinateTransform> ctl = new CoordinateTransformList<CoordinateTransform>();
+		try 
+		{
+			final FileReader fr = new FileReader(filename);
+			final BufferedReader br = new BufferedReader(fr);
+			String line = null;
+			while ((line = br.readLine()) != null) 
+			{
+				int index = -1;
+				if( (index = line.indexOf("class=")) != -1)
+				{
+					// skip "class"
+					index+= 5;
+					// read coordinate transform class name
+					final int index2 = line.indexOf("\"", index+2); 
+					final String ct_class = line.substring(index+2, index2);
+					final CoordinateTransform ct = (CoordinateTransform) Class.forName(ct_class).newInstance();
+					// read coordinate transform info
+					final int index3 = line.indexOf("=", index2+1);
+					final int index4 = line.indexOf("\"", index3+2); 
+					final String data = line.substring(index3+2, index4);
+					ct.init(data);
+					ctl.add(ct);
+				}
+			}
+			br.close();
+
+		} catch (FileNotFoundException e) {
+			System.err.println("File not found exception" + e);
+
+		} catch (IOException e) {
+			System.err.println("IOException exception" + e);
+
+		} catch (NumberFormatException e) {
+			System.err.println("Number format exception" + e);
+
+		} catch (InstantiationException e) {
+			System.err.println("Instantiation exception" + e);
+
+		} catch (IllegalAccessException e) {
+			System.err.println("Illegal access exception" + e);
+
+		} catch (ClassNotFoundException e) {
+			System.err.println("Class not found exception" + e);
+
+		}
+		return ctl;
+	}
+
+	// adapted from register_virtual_stack_slices/src/main/java/register_virtual_stack/Register_Virtual_Stack_MT.applyTransformAndSave()
+	// to only return the bounding box
+	public static Rectangle getBoundingBox(
+			final int width,
+			final int height,
+			final CoordinateTransform transform)
+	{
+		// Open next image
+		//final ImagePlus imp2 = readImage(source_dir + file_name);
+
+		// Calculate transform mesh
+		final TransformMesh mesh = new TransformMesh(transform, 32, width, height);
+		//TransformMeshMapping mapping = new TransformMeshMapping(mesh);
+
+		// Create interpolated deformed image with black background
+		//imp2.getProcessor().setValue(0);
+		//final ImageProcessor ip2 = interpolate ? mapping.createMappedImageInterpolated(imp2.getProcessor()) : mapping.createMappedImage(imp2.getProcessor()); 
+		//imp2.setProcessor(imp2.getTitle(), ip2);
+
+		//imp2.show();
+
+		// Accumulate bounding boxes, so in the end they can be reopened and re-saved with an enlarged canvas.
+		final Rectangle currentBounds = mesh.getBoundingBox();
+		return currentBounds;
+		//bounds[i] = currentBounds;
+
+		// Save target image
+		//return new FileSaver(imp2).saveAsTiff(makeTargetPath(target_dir, file_name));
+	}
 }