ij.measure.Calibration Java Examples

/**
 * Identifies labels within image and computes an instance of the generic
 * type T for each region in input label image.
 * 
 * @param image
 *            a label image (8, 16 or 32 bits)
 * @param calib
 *            the spatial calibration of the image
 * @return a map between the region label and the result of analysis for
 *         each region
 */
public Map<Integer, T> analyzeRegions(ImageProcessor image, Calibration calib)
{
	// extract particle labels
	fireStatusChanged(this, "Find Labels");
	int[] labels = LabelImages.findAllLabels(image);
	int nLabels = labels.length;
	
	// compute analysis result for each label
	fireStatusChanged(this, "Analyze regions");
	T[] results = analyzeRegions(image, labels, calib);

	// encapsulate into map
	fireStatusChanged(this, "Convert to map");
	Map<Integer, T> map = new TreeMap<Integer, T>();
	for (int i = 0; i < nLabels; i++)
	{
		map.put(labels[i], results[i]);
	}

       return map;
}
 

/**
 * Test method for {@link inra.ijpb.measure.region2d.Convexity#analyzeRegions(ij.process.ImageProcessor, int[], ij.measure.Calibration)}.
 */
@Test
public final void testAnalyzeRegionsImageProcessorIntArrayCalibration()
{
	ImageProcessor image = new ByteProcessor(8, 8);
	image.set(1, 1, 255);
	image.set(5, 1, 255);
	image.set(1, 5, 255);
	image.set(5, 5, 255);
	int[] labels = new int[] {255};
	Calibration calib = new Calibration();
	
	Convexity algo = new Convexity();
	Convexity.Result[] results = algo.analyzeRegions(image, labels, calib);
	assertEquals(results.length, 1);
	
	Convexity.Result res1 = results[0];
	assertEquals(4, res1.area, .01);
	assertEquals(25, res1.convexArea, .01);
	assertEquals(4.0/25.0, res1.convexity, .01);
}
 

/**
 * Calculate the neighborhood size based on the calibration of the image.
 * @param calibration image calibration
 */
private void calculateNeighborhoodOffsets(Calibration calibration) 
{
	double max = calibration.pixelDepth;
	if(calibration.pixelHeight > max)
		max = calibration.pixelHeight;
	if(calibration.pixelWidth > max)
		max = calibration.pixelWidth;

	this.x_offset = ((int) Math.round(max / calibration.pixelWidth) > 1) ? 
						(int) Math.round(max / calibration.pixelWidth) : 1;
	this.y_offset = ((int) Math.round(max / calibration.pixelHeight) > 1) ? 
						(int) Math.round(max / calibration.pixelHeight) : 1;
	this.z_offset = ((int) Math.round(max / calibration.pixelDepth) > 1) ? 
						(int) Math.round(max / calibration.pixelDepth) : 1;

	if(debug)
	{
		IJ.log("x_offset = " + this.x_offset);
		IJ.log("y_offset = " + this.y_offset);
		IJ.log("z_offset = " + this.z_offset);
	}

}
 

/**
 * Computes the Look-up table used to measure mean breadth within 3D images.
 * 
 * @param calib
 *            the spatial calibration of image
 * @param nDirs
 *            the number of directions (3 or 13)
 * @param conn2d
 *            the connectivity to use on square faces of plane sections (4 or 8)
 * @return a look-up table with 256 entries
 */
public static final double[] meanBreadthLut(Calibration calib, int nDirs, int conn2d)
{
    if (nDirs == 3)
    {
        return meanBreadthLutD3(calib, conn2d);
    }
    else if (nDirs == 13)
    {
        return meanBreadthLutD13(calib, conn2d);
    }
    else
    {
        throw new IllegalArgumentException("Number of directions should be either 3 or 13, not " + nDirs);
    }
}
 

/** Make a mesh as a calibrated list of 3D triangles.*/
static private List<Point3f> makeTriangles(final Profile[] p, final double scale) {
	try {
		final VectorString2D[] sv = new VectorString2D[p.length];
		boolean closed = true; // dummy initialization
		final Calibration cal = p[0].getLayerSet().getCalibrationCopy();
		cal.pixelWidth *= scale;
		cal.pixelHeight *= scale;
		for (int i=0; i<p.length; i++) {
			if (0 == p[i].n_points) continue;
			if (0 == i) closed = p[i].closed;
			else if (p[i].closed != closed) {
				Utils.log2("All profiles should be either open or closed, not mixed.");
				return null;
			}
			sv[i] = p[i].getPerimeter2D(cal);
		}
		return SkinMaker.generateTriangles(sv, -1, -1, closed);
	} catch (final Exception e) {
		IJError.print(e);
	}
	return null;
}
 

/**
 * Test method for {@link inra.ijpb.measure.region3d.InterfaceSurfaceArea#process(ij.ImageStack, int, int, ij.measure.Calibration)}.
 */
@Test
public final void testProcess_SingleVoxel()
{
    ImageStack stack = ImageStack.create(5, 5, 5, 8);
    for (int z = 0; z < 5; z++)
    {
        for (int y = 0; y < 5; y++)
        {
            for (int x = 0; x < 5; x++)
            {
                stack.setVoxel(x, y, z, 1);
            }
        }
    }
    stack.setVoxel(2, 2, 2, 2);
   
    Calibration calib = new Calibration();
    
    InterfaceSurfaceArea algo = new InterfaceSurfaceArea();
    
    // labels 1-2  -> R = 1, S ~= 3.14 
    double S12 = algo.process(stack, 1, 2, calib);
    double expS12 = Math.PI * 4 * .5 * .5; 
    assertEquals(expS12, S12, expS12 * .05);
}
 

@Override
public ResultsTable measure(ResultsTable rt) {
	if (-1 == n_points) setupForDisplay(); //reload
	if (0 == n_points) return rt;
	if (null == rt) rt = Utils.createResultsTable("Polyline results", new String[]{"id", "length", "name-id"});
	// measure length
	double len = 0;
	final Calibration cal = layer_set.getCalibration();
	if (n_points > 1) {
		final VectorString3D vs = asVectorString3D();
		vs.calibrate(cal);
		len = vs.computeLength(); // no resampling
	}
	rt.incrementCounter();
	rt.addLabel("units", cal.getUnit());
	rt.addValue(0, this.id);
	rt.addValue(1, len);
	rt.addValue(2, getNameId());
	return rt;
}
 

/** Returns a listing of all balls contained here, one per row with index, x, y, z, and radius, all calibrated.
 * 'name-id' is a column that displays the title of this Ball object only when such title is purely a number.
 */
@Override
public ResultsTable measure(ResultsTable rt) {
	if (-1 == n_points) setupForDisplay(); //reload
	if (0 == n_points) return rt;
	if (null == rt) rt = Utils.createResultsTable("Ball results", new String[]{"id", "index", "x", "y", "z", "radius", "name-id"});
	final Object[] ob = getTransformedData();
	final double[][] p = (double[][])ob[0];
	final double[] p_width = (double[])ob[1];
	final Calibration cal = layer_set.getCalibration();
	for (int i=0; i<n_points; i++) {
		rt.incrementCounter();
		rt.addLabel("units", cal.getUnit());
		rt.addValue(0, this.id);
		rt.addValue(1, i+1);
		rt.addValue(2, p[0][i] * cal.pixelWidth);
		rt.addValue(3, p[1][i] * cal.pixelHeight);
		rt.addValue(4, layer_set.getLayer(p_layer[i]).getZ() * cal.pixelWidth);
		rt.addValue(5, p_width[i] * cal.pixelWidth);
		rt.addValue(6, getNameId());
	}
	return rt;
}
 

/**
 * Test method for {@link inra.ijpb.measure.IntrinsicVolumes3D#surfaceAreaDensity(ij.ImageStack, ij.measure.Calibration, int)}.
 */
@Test
public void testSurfaceAreaDensity()
{
    String fileName = getClass().getResource("/files/microstructure3D_10x10x10.tif").getFile();
    ImagePlus imagePlus = IJ.openImage(fileName);
    
    // basic check up
    assertNotNull(imagePlus);
    assertTrue(imagePlus.getStackSize() > 0);

    // get image stack and calibration
    ImageStack image = imagePlus.getStack();
    Calibration calib = imagePlus.getCalibration();

    // compare with a value computed with Matlab (2018.09.11)
    double surfaceDensity = IntrinsicVolumes3D.surfaceAreaDensity(image, calib, 13);
    assertEquals(0.56, surfaceDensity, .001);
}
 

/**
 * Test method for {@link inra.ijpb.measure.IntrinsicVolumes2D#perimeter(ij.process.ImageProcessor, Calibration, int)}.
 */
@Test
public final void testPerimeter_smallSquare_D4()
{
	// create a binary image containing a square
	ImageProcessor image = new ByteProcessor(8, 8);
	for (int y = 2; y < 6; y++)
	{
		for (int x = 2; x < 6; x++)
		{
			image.set(x, y, 255);
		}
	}
	
	double perim = IntrinsicVolumes2D.perimeter(image, new Calibration(), 4);
	assertEquals(14.0582, perim, .01);
}
 

/**
 * Test method for {@link inra.ijpb.measure.region2d.LargestInscribedCircle#analyzeRegions(ij.ImagePlus)}.
 */
@Test
public final void testAnalyzeRegionsImagePlus_Rectangles()
{
	ImageProcessor image = new ByteProcessor(14, 9);
	image.set(1, 1, 1);
	fillRect(image, 3, 4, 1, 2, 2);		// radius 2
	fillRect(image, 1, 1+3, 4, 4+3, 3); // radius 4
	fillRect(image, 6, 6+6, 1, 1+6, 4); // radius 7
	
	LargestInscribedCircle op = new LargestInscribedCircle();
	ResultsTable table = op.createTable(op.analyzeRegions(image, new Calibration()));
	
	assertEquals(1, table.getValue("InscrCircle.Radius", 0), .1);
	assertEquals(1, table.getValue("InscrCircle.Radius", 1), .1);
	assertEquals(2, table.getValue("InscrCircle.Radius", 2), .1);
	assertEquals(4, table.getValue("InscrCircle.Radius", 3), .1);
}
 

/**
 * Test method for {@link inra.ijpb.measure.IntrinsicVolumes2D#perimeter(ij.process.ImageProcessor, Calibration, int)}.
 */
@Test
public final void testPerimeter_disk_D4()
{
	double radius = 16.0;

	// create a binary image containing a square
	ImageProcessor image = new ByteProcessor(40, 40);
	for (int y = 0; y < 40; y++)
	{
		double y2 = (y - 20.2);
		for (int x = 0; x < 40; x++)
		{
			double x2 = (x - 20.3);
			image.set(x, y, Math.hypot(x2, y2) < radius ? 255 : 0);
		}
	}
	
	// compute perimeter with default (1,1) calibration
	Calibration calib = new Calibration();
	double perim = IntrinsicVolumes2D.perimeter(image, calib, 4);
	
	// check to expected value with a tolerance of 5 percents
	double exp = 2 * Math.PI * radius;
	assertEquals(exp, perim, exp * 0.05);
}
 

/**
 * Test method for {@link inra.ijpb.measure.IntrinsicVolumes3D#surfaceArea(ij.ImageStack, ij.measure.Calibration, int)}.
 */
@Test
public final void testSurfaceArea_SmallCube_D3()
{
    ImageStack image = ImageStack.create(4, 4, 4, 8);
    image.setVoxel(1, 1, 1, 255);
    image.setVoxel(2, 1, 1, 255);
    image.setVoxel(1, 2, 1, 255);
    image.setVoxel(2, 2, 1, 255);
    image.setVoxel(1, 1, 2, 255);
    image.setVoxel(2, 1, 2, 255);
    image.setVoxel(1, 2, 2, 255);
    image.setVoxel(2, 2, 2, 255);
    Calibration calib = new Calibration();
    
    double surface = IntrinsicVolumes3D.surfaceArea(image, calib, 3);
    
    double exp = 16.0;
    assertEquals(exp, surface, 16.0*0.2);
}
 

/**
 * Test method for {@link inra.ijpb.measure.IntrinsicVolumes2D#perimeter(ij.process.ImageProcessor, Calibration, int)}.
 */
@Test
public final void testPerimeters_smallSquare_D2()
{
	// create a binary image containing a square
	ImageProcessor image = new ByteProcessor(8, 8);
	for (int y = 2; y < 6; y++)
	{
		for (int x = 2; x < 6; x++)
		{
			image.set(x, y, 255);
		}
	}
	
	// compute perimeter with default (1,1) calibration
	Calibration calib = new Calibration();
	int[] labels = new int[] {255};
	double[] perims = IntrinsicVolumes2D.perimeters(image, labels, calib, 2);
	
	assertEquals(12.5664, perims[0], .01);
}
 

/**
 * Test method for {@link inra.ijpb.measure.IntrinsicVolumes3D#meanBreadth(ij.ImageStack, ij.measure.Calibration, int)}.
 */
@Test
public final void testMeanBreadth_SmallCube_D3()
{
    ImageStack image = ImageStack.create(4, 4, 4, 8);
    image.setVoxel(1, 1, 1, 255);
    image.setVoxel(2, 1, 1, 255);
    image.setVoxel(1, 2, 1, 255);
    image.setVoxel(2, 2, 1, 255);
    image.setVoxel(1, 1, 2, 255);
    image.setVoxel(2, 1, 2, 255);
    image.setVoxel(1, 2, 2, 255);
    image.setVoxel(2, 2, 2, 255);
    Calibration calib = new Calibration();
    
    double meanBreadth = IntrinsicVolumes3D.meanBreadth(image, calib, 3, 8);
    
    double exp = 2.0;
    assertEquals(exp, meanBreadth, exp * 0.2);
}
 

@Test
public final void testGeodesicDiameter_FiveTouchingRects_Borgefors()
{
	ByteProcessor labelImage = new ByteProcessor(17, 11);
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 9; j++)
		{
			labelImage.set(i +  1, j + 1, 1); 
			labelImage.set(j +  4, i + 1, 2); 
			labelImage.set(j +  4, i + 4, 3); 
			labelImage.set(j +  4, i + 7, 4); 
			labelImage.set(i + 13, j + 1, 5); 
		}
	}
	
	GeodesicDiameter algo = new GeodesicDiameter(ChamferWeights.BORGEFORS);
	int[] labels = new int[]{1, 2, 3, 4, 5};
	GeodesicDiameter.Result[] geodDiams = algo.analyzeRegions(labelImage, labels, new Calibration());
	
	for (int i = 0; i < 5; i++)
	{
		assertEquals((26.0/3.0)+1.41, geodDiams[i].diameter, .1);
	}
}
 

/**
 * Test method for {@link inra.ijpb.measure.region3d.EquivalentEllipsoid#equivalentEllipsoids(ij.ImageStack, int[], ij.measure.Calibration)}.
 */
@Test
public final void testEquivalentEllipsoid3()
{
    String fileName = getClass().getResource("/files/ellipsoid_A45_B35_C25_T30_P20_R10.tif").getFile();
    ImagePlus imagePlus = IJ.openImage(fileName);
    assertNotNull(imagePlus);
    ImageStack image = imagePlus.getStack();

    Calibration calib = imagePlus.getCalibration();
    Ellipsoid elli = EquivalentEllipsoid.equivalentEllipsoids(image, new int[] {255}, calib)[0];
    
    assertEquals(45, elli.radius1(), .1);
    assertEquals(35, elli.radius2(), .1);
    assertEquals(25, elli.radius3(), .1);
    assertEquals(30, elli.phi(), 1.0);
    assertEquals(20, elli.theta(), 1.0);
    assertEquals(10, elli.psi(), 1.0);
}
 

/**
 * Display the result of bounding box extraction as overlay on a given
 * image.
 * 
 * @param target
 *            the ImagePlus used to display result
 * @param results
 *            the associative map between region label and bounding box
 */
private void showResultsAsOverlay(Map<Integer, Box2D> results, ImagePlus target)	
{
	// get spatial calibration of target image
	Calibration calib = target.getCalibration();
	
	// create overlay
	Overlay overlay = new Overlay();
	Roi roi;
	
	// add each box to the overlay
	for (int label : results.keySet()) 
	{
		// Coordinates of inscribed circle, in pixel coordinates
		Box2D box = results.get(label);
		box = uncalibrate(box, calib);
		roi = createRoi(box);
		
		// draw inscribed circle
		roi.setStrokeColor(Color.BLUE);
		overlay.add(roi);
	}
	
	target.setOverlay(overlay);
}
 

/**
 * Test method for {@link inra.ijpb.measure.region2d.OrientedBoundingBox2D#analyzeRegions(ij.process.ImageProcessor, int[], ij.measure.Calibration)}.
 */
@Test
public final void testAnalyzeRegions_circles()
{
	String fileName = getClass().getResource("/files/circles.tif").getFile();
	ImagePlus imagePlus = IJ.openImage(fileName);
	assertNotNull(imagePlus);
	
	ImageProcessor image = imagePlus.getProcessor();
	int[] labels = new int[] {255};
	Calibration calib = new Calibration();

	OrientedBoundingBox2D algo = new OrientedBoundingBox2D();
	OrientedBox2D[] boxes = algo.analyzeRegions(image, labels, calib);
	
	assertEquals(1, boxes.length);
	OrientedBox2D box = boxes[0];
	
	// Length of oriented box
	assertEquals(272.23, box.length(), .05);
	// width of oriented box
	assertEquals(108.86, box.width(), .05);
}
 

private static final ArrayList<Point2D> calibrate(ArrayList<Point2D> points, Calibration calib)
{
	if (!calib.scaled())
	{
		return points;
	}
	
	ArrayList<Point2D> res = new ArrayList<Point2D>(points.size());
	for (Point2D point : points)
	{
		double x = point.getX() * calib.pixelWidth + calib.xOrigin;
		double y = point.getY() * calib.pixelHeight + calib.yOrigin;
		res.add(new Point2D.Double(x, y));
	}
	return res;
}
 

/**
 * Display the result of oriented bounding box extraction as overlay on a
 * given image.
 * 
 * @param target
 *            the ImagePlus used to display result
 * @param results
 *            the associative map between region label and bounding box
 */
private void showResultsAsOverlay(Map<Integer, OrientedBox2D> results, ImagePlus target)	
{
	// get spatial calibration of target image
	Calibration calib = target.getCalibration();
	
	// create overlay
	Overlay overlay = new Overlay();
	Roi roi;
	
	// add each box to the overlay
	for (int label : results.keySet()) 
	{
		// Coordinates of inscribed circle, in pixel coordinates
		OrientedBox2D box = results.get(label);
		roi = createUncalibratedRoi(box, calib);
		
		// draw inscribed circle
		roi.setStrokeColor(Color.GREEN);
		overlay.add(roi);
	}
	
	target.setOverlay(overlay);
}
 

/**
 * Test method for {@link inra.ijpb.measure.region2d.Convexity#analyzeRegions(ij.process.ImageProcessor, int[], ij.measure.Calibration)}.
 */
@Test
public final void testAnalyzeRegionsImageProcessorIntArrayCalibration_FourLabels()
{
    // create a 8x8 image with four regions corresponding to the corners of a 3x3 square
    ImageProcessor image = new ByteProcessor(8, 8);
    for (int y = 0; y < 4; y+=2)
    {
        for (int x = 0; x < 4; x+=2)
        {
            image.set(x + 1, y + 1, 2);
            image.set(x + 4, y + 1, 3);
            image.set(x + 1, y + 4, 4);
            image.set(x + 4, y + 4, 7);
        }
    }
    
    int[] labels = new int[] {2, 3, 4, 7};
    Calibration calib = new Calibration();
    
    Convexity algo = new Convexity();
    Convexity.Result[] results = algo.analyzeRegions(image, labels, calib);
    assertEquals(results.length, 4);
    
    for (int i = 0; i < 4; i++)
    {
        Convexity.Result res = results[0];
        assertEquals(4, res.area, .01);
        assertEquals(9, res.convexArea, .01);
        assertEquals(4.0/9.0, res.convexity, .01);
    }
}
 

public void drawDiameters(ImagePlus target, Map<Integer, PointPair2D> geodDiams)
{
	Overlay overlay = new Overlay();
	Calibration calib = target.getCalibration();
	
	for (PointPair2D result : geodDiams.values())
	{
		Roi roi = createDiametersRoi(result, calib);
	    roi.setStrokeColor(Color.BLUE);
	    overlay.add(roi);
	}

	target.setOverlay(overlay);
}
 

/**
 * Test method for {@link inra.ijpb.measure.IntrinsicVolumes2D#eulerNumberDensity(ij.process.ImageProcessor, int)}.
 */
@Test
public final void testPerimeterDensity_OhserMecklich_D4()
{
    ImageProcessor image = createOhserMuecklichImage();
    Calibration calib = new Calibration();
    
    double density = IntrinsicVolumes2D.perimeterDensity(image, calib, 4);
    
    assertEquals(0.5, density, .05);
}
 

public ResultsTable process(ImagePlus inputImage, int nDirs) 
  {
      // Check validity of parameters
      if (inputImage==null) 
          return null;

      if (debug) 
      {
      	System.out.println("Compute Crofton perimeter on image '" 
      			+ inputImage.getTitle());
      }
      
      ImageProcessor proc = inputImage.getProcessor();
      
      // Extract spatial calibration
      Calibration cal = inputImage.getCalibration();
      double[] resol = new double[]{1, 1};
      if (cal.scaled()) 
      {
      	resol[0] = cal.pixelWidth;
      	resol[1] = cal.pixelHeight;
      }

      ResultsTable results = GeometricMeasures2D.analyzeRegions(proc, resol);

// return the created array
return results;
  }
 

private final static double[] interceptsContributionsD13(Calibration calib)
{
    // distances between a voxel and its neighbors.
    // di refer to orthogonal neighbors
    // dij refer to neighbors on the same plane 
    // dijk refer to the opposite voxel in a tile
    double d1 = calib.pixelWidth;
    double d2 = calib.pixelHeight;
    double d3 = calib.pixelDepth;
    double vol = d1 * d2 * d3;
    
    double d12 = Math.hypot(d1, d2);
    double d13 = Math.hypot(d1, d3);
    double d23 = Math.hypot(d2, d3);
    double d123= Math.hypot(d12, d3);

    // direction weights corresponding to area of Voronoi partition on the
    // unit sphere, when germs are the 26 directions on the unit cube
    // Sum of (c1+c2+c3 + c4*2+c5*2+c6*2 + c7*4) equals 1.
    double[] weights = IntrinsicVolumes3DUtils.directionWeights3d13(calib);
    
    // contribution to number of intersections for each direction
    double[] kei = new double[7];
    
    kei[0] = 4 * weights[0] * vol / d1   / 8.0;
    kei[1] = 4 * weights[1] * vol / d2   / 8.0;
    kei[2] = 4 * weights[2] * vol / d3   / 8.0;
    kei[3] = 4 * weights[3] * vol / d12  / 4.0;
    kei[4] = 4 * weights[4] * vol / d13  / 4.0;
    kei[5] = 4 * weights[5] * vol / d23  / 4.0;
    kei[6] = 4 * weights[6] * vol / d123 / 2.0; 

    return kei;     
}
 

public void setOutputSize(int tgtW, int tgtH) {
	renderer.setTargetSize(tgtW, tgtH);
	Calibration cal = outputImage.getCalibration();
	CombinedTransform trans = renderer.getRenderingState().getFwdTransform();
	trans.adjustOutputCalibration(cal);
	renderer.getScalebar().setDefaultLength(trans.getOutputSpacing()[0] / trans.getScale());
	push(tgtW, tgtH);
}
 

/**
 * Test method for {@link inra.ijpb.measure.IntrinsicVolumes3D#surfaceAreas(ij.ImageStack, int[], ij.measure.Calibration, int)}.
 */
@Test
public final void testSurfaceAreas_SingleBall_D3_HalfResolY()
{
    ImageStack image = createBallImage_HalfResolY();
    int[] labels = new int[] {255};
    Calibration calib = new Calibration();
    calib.pixelHeight = 2.0;
    
    double[] surfaces = IntrinsicVolumes3D.surfaceAreas(image, labels, calib, 3);
    
    double exp = 5026.0;
    assertEquals(1, surfaces.length);
    assertEquals(exp, surfaces[0], exp * 0.02);
}
 

static public void zeroPad(final ImagePlus imp) {
	final ImageStack stack = imp.getStack();
	final int w = stack.getWidth();
	final int h = stack.getHeight();
	final int d = stack.getSize();
	final int type = imp.getType();
	// create new stack
	final ImageStack st = new ImageStack(w+2, h+2);

	// retrieve 1st processor
	ImageProcessor old = stack.getProcessor(1);

	// enlarge it and add it as a first slide.
	ImageProcessor ne = createProcessor(w+2, h+2, type);
	st.addSlice("", ne);

	// now do the same for all slices in the old stack
	for(int z = 0; z < d; z++) {
		old = stack.getProcessor(z+1);
		ne = createProcessor(w+2, h+2, type);
		ne.insert(old, 1, 1);
		st.addSlice(Integer.toString(z+1), ne);
	}

	// now add an empty new slice
	ne = createProcessor(w+2, h+2, type);
	st.addSlice(Integer.toString(d+1), ne);

	imp.setStack(null, st);

	// update the origin
	final Calibration cal = imp.getCalibration();
	cal.xOrigin -= cal.pixelWidth;
	cal.yOrigin -= cal.pixelHeight;
	cal.zOrigin -= cal.pixelDepth;
	imp.setCalibration(cal);
}
 

/**
 * Test method for {@link inra.ijpb.measure.IntrinsicVolumes3D#surfaceAreas(ij.ImageStack, int[], ij.measure.Calibration, int)}.
 */
@Test
public final void testSurfaceAreas_TouchingLabels_D3()
{
    ImageStack image = ImageStack.create(9, 9, 9, 8);
    for (int i = 0; i < 3; i++)
    {
        for (int j = 0; j < 3; j++)
        {
            for (int k = 0; k < 3; k++)
            {
                image.setVoxel(    i,     j,     k,  1);
                image.setVoxel(i + 3,     j,     k,  2);
                image.setVoxel(    i, j + 3,     k,  3);
                image.setVoxel(i + 3, j + 3,     k,  4);
                image.setVoxel(    i,     j, k + 3,  5);
                image.setVoxel(i + 3,     j, k + 3,  6);
                image.setVoxel(    i, j + 3, k + 3,  7);
                image.setVoxel(i + 3, j + 3, k + 3,  8);
            }
        }
    }
    int[] labels = new int[] {1, 2, 3, 4, 5, 6, 7, 8};
    Calibration calib = new Calibration();
    
    double[] surfaces = IntrinsicVolumes3D.surfaceAreas(image, labels, calib, 3);
    double exp = 36.0;
    assertEquals(8, surfaces.length);
    for (int i = 0; i < 8; i++)
    {
        assertEquals(exp, surfaces[i], exp * 0.2);
    }
}