HHTracker.java

Go to the documentation of this file.

package theba.trackers;

import java.awt.Point;
import java.awt.Rectangle;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Iterator;
import java.util.LinkedList;

import javax.swing.JOptionPane;
import javax.swing.JToggleButton;

import theba.core.BackTrack;
import theba.core.FiberSlice;
import theba.core.ImageFunctions;
import theba.core.Lumen;
import theba.core.LumenCandidate;
import theba.core.RegionMask2D;
import theba.core.Seed;
import theba.core.Stack;
import theba.core.Tracker;
import theba.core.gui.ThebaGUI;
import theba.core.io.SliceWriter;
import theba.core.math.DistancePath;
import theba.core.math.EDM;
import theba.core.math.Point3D;
import watershed.Watershed_Algorithm;
import theba.descriptors.AvgYoungCurvature;
import theba.descriptors.CircularityDescriptor;
import theba.descriptors.ConvexityDescriptor;
import theba.descriptors.EccentricityDescriptor;
import theba.descriptors.InternalRegionDescriptor;
import theba.descriptors.SolidityDescriptor;
import theba.descriptors.YoungBendingEnergy;


public class HHTracker extends Tracker {
        private LinkedList fibers = new LinkedList();

        private JToggleButton seedButton;

        private LinkedList<Seed> seedQueue;

        private LinkedList<LumenCandidate> seeds = new LinkedList<LumenCandidate>();

        boolean seedStillValid = true;

        private short currentId;

        private short[] map;

        private short[] map2;

        // configurable options
        private boolean automerge = false;

        private int minimumFiberSize = 50;

        private boolean removeShortFibers = false;

        @SuppressWarnings("unchecked")
        @Override
        public void track() {
                int width = control.getStack().getWidth();
                int height = control.getStack().getHeight();
                int depth = control.getStack().getDepth();
                automerge = (1 == control.getPreferences().getInt("autoMerge", 0));
                Point[][] skeleton = new Point[control.MAX_FIBERS][depth];
                Collections.sort(seeds);
                Point3D s = null;
                Iterator<LumenCandidate> iter = seeds.iterator();
                int count = 0;
                keepTracking = true;
                while (iter.hasNext() && keepTracking && !control.isStopped()) {
                        while (paused)
                                Thread.yield();
                        s = iter.next();
                        seedQueue = new LinkedList<Seed>();
                        count++;
                        System.out.println("\nTracking seed " + count + " of "
                                        + seeds.size() + " starting at " + s);
                        short[] data = control.getSlice(s.z);
                        if (data[s.x + s.y * width] != 0) {
                                continue;
                        }
                        currentId = (control.getNewFiberId());
                        ImageFunctions.delete3D(control.getStack(), s.x, s.y, s.z);
                        seedQueue.addLast(new Seed(s.x, s.y, s.z, 1, null, null, 0));
                        log("***************************");
                        log("Tracking seed with origin " + s);
                        StringBuffer results = new StringBuffer();
                        // results.append("Tracked lumens =" + (control.MAX_FIBERS -
                        // idList.size() - 2) + "\n");
                        results.append("Tracking seed  " + count + " of " + seeds.size()
                                        + "\n from " + s);
                        control.showResults(results, "Lumen tracking");
                        short[] lumenMask = new short[data.length];
                        int slice = s.z;
                        ImageFunctions.floodFill2D(s.x, s.y, width, height,
                                        data, lumenMask, currentId);
                        short[] lcopy = lumenMask.clone();
                        seedStillValid = true;
                        int length = 0;
                        // Trace forward from seed
                        while (slice < depth && seedStillValid && keepTracking) {
                                control.setProgress(slice);
                                removeCrack(slice, lumenMask);
                                removeStuff(slice, lumenMask);
                                if (automerge)
                                        skeleton[currentId][slice] = ImageFunctions
                                        .getAveragePoint(lumenMask, width, height);
                                if (slice % 2 == 0)
                                        control.showImage(control.getSlice(slice));
                                slice++;
                        }
                        seedStillValid = true;
                        slice = s.z - 1;
                        lumenMask = lcopy;
                        // Trace backward from seed
                        while (slice >= 0 && seedStillValid && keepTracking
                                        && !control.isStopped()) {
                                control.setProgress(slice);
                                removeCrack(slice, lumenMask);
                                removeStuff(slice, lumenMask);
                                if (slice % 2 == 0)
                                        control.showImage(control.getSlice(slice));
                                if (automerge)
                                        skeleton[currentId][slice] = ImageFunctions
                                        .getAveragePoint(lumenMask, width, height);
                                slice--;
                        }
                        control.showImage(control.getSlice(slice + 1));
                        if (automerge)
                                checkMerges(skeleton);
                        // We must recalculate length, since the total length may have grown
                        // after
                        // merges have been applied!
                        length = 0;
                        for (int z = 0; z < depth; z++) {
                                if (skeleton[currentId][z] != null)
                                        length++;
                        }
                        if (removeShortFibers) {
                                if (length < minimumFiberSize) {
                                        System.out.println("Fibre rejected, reclaiming space!");
                                        for (int z = 0; z < depth; z++) {
                                                short[] pixels = control.getSlice(z);
                                                for (int i = 0; i < pixels.length; i++) {
                                                        if (pixels[i] == currentId)
                                                                pixels[i] = control.INVALID;
                                                }
                                        }
                                        control.releaseFiberId(currentId);
                                        skeleton[currentId] = new Point[depth];
                                }
                        }
                }
                control.setProgressComplete();
                control.updateImage();
                seeds.clear(); // Clear seeds marks
                control.showResults("Tracking lumens complete");
        }

        private final class FindAllCandidatesAction implements Runnable {
                public void run() {
                        int depth = control.getStack().getDepth();
                        String s = JOptionPane.showInputDialog("How many slices to skip?",
                        "10");
                        try {
                                final int num = Integer.parseInt(s);
                                for (int z = 0; z < depth; z += num) {
                                        control.setProgress(z);
                                        findCandidates(z);
                                }
                        } catch (NumberFormatException ee) {
                        }
                }
        }

        private final class FindCandidatesAction implements Runnable {
                public void run() {
                        findCandidates(control.currentSlice());
                }
        }

        public HHTracker(ThebaGUI f) {
                super(f);
        }

        @Override
        public void setup() {
                seedButton = new JToggleButton("Mark seed");
                control.addToolbarButton(seedButton);
                Runnable lumenAction = new Runnable() {
                        public void run() {
                                track();
                        }
                };
                Runnable wallAction2 = new Runnable() {
                        public void run() {
                                traceWalls3d();
                        }
                };
                FindCandidatesAction findCandidates = new FindCandidatesAction();
                FindAllCandidatesAction findAllCandidates = new FindAllCandidatesAction();
                Runnable clearAction = new Runnable() {
                        public void run() {
                                seeds.clear();
                        }
                };
                control.addMenuItem("Track lumens", lumenAction);
                control.addMenuItem("Track walls (2D)", wallAction2);
                control.addMenuSeparator();
                control.addMenuItem("Find seeds in this slice", findCandidates);
                control.addMenuItem("Find seeds for every n'th slice",
                                findAllCandidates);
                control.addMenuSeparator();
                control.addMenuItem("Clear seeds", clearAction);
        }

        static final short max(short x, short y) {
                if ((x) > (y))
                        return x;
                return y;
        }

        public static void dilate3d_toWhite(Stack voxels, long[] counts) {
                int width = voxels.getWidth();
                int height = voxels.getHeight();
                int depth = voxels.getDepth();
                short[] prevdata = new short[width * height];
                short[] data = new short[width * height];
                short[] nextdata = new short[width * height];
                ThebaGUI control = ThebaGUI.getInstance();
                System.arraycopy(voxels.getSlice(1), 0, nextdata, 0, nextdata.length);
                System.arraycopy(voxels.getSlice(0), 0, data, 0, data.length);
                for (int z = 0; z < depth; z++) {
                        short[] origdata = voxels.getSlice(z);
                        control.setProgress(z);
                        if (control.isStopped())
                                return;
                        int index = 0;
                        for (int y = 0; y < height; y++) {
                                for (int x = 0; x < width; x++) {
                                        short max = data[index];
                                        if (max == 255) {
                                                /*
                                                 * Consider only propagating pixels that are not marked
                                                 * by a negative value in the counts table, this
                                                 * indicates that they are beyond their detected wall
                                                 * size and should not be propagated further
                                                 */
                                                if (x > 0)
                                                        max = max(max, data[index - 1]);
                                                if (x < width - 1)
                                                        max = max(max, data[index + 1]);
                                                if (y > 0)
                                                        max = max(max, data[index - width]);
                                                if (y < height - 1)
                                                        max = max(max, data[index + width]);
                                                if (z > 0)
                                                        max = max(max, prevdata[index]);
                                                if (z < depth - 1)
                                                        max = max(max, nextdata[index]);
                                                int id = max & (~(1 << 12));
                                                if (max > data[index] && counts[id] != -1) {
                                                        origdata[index] = (short) (max | (1 << 12));
                                                        counts[id]++;
                                                }
                                        }
                                        index++;
                                }
                        }
                        if (z > 0)
                                System.arraycopy(data, 0, prevdata, 0, prevdata.length);
                        System.arraycopy(nextdata, 0, data, 0, data.length);
                        if (z < depth - 2)
                                System.arraycopy(voxels.getSlice(z + 2), 0, nextdata, 0,
                                                nextdata.length);
                }
        }

        private void traceWalls3d() {
                Stack stack = control.getStack();
                int width = stack.getWidth();
                int height = stack.getHeight();
                int depth = stack.getDepth();
                long[] startvals = new long[control.MAX_FIBERS];
                long[] lumenvals = new long[control.MAX_FIBERS];
                for (int i = 0; i < 15; i++) {
                        // reset all non-negative values
                        for (int k = 0; k < lumenvals.length; k++) {
                                if (lumenvals[k] > 0)
                                        lumenvals[k] = 0;
                        }
                        dilate3d_toWhite(control.getStack(), lumenvals);
                        for (int k = 0; k < lumenvals.length; k++) {
                                if (lumenvals[k] > 0) {
                                        if (i == 0)
                                                startvals[k] = lumenvals[k];
                                        if (lumenvals[k] < startvals[k] / 4) {
                                                lumenvals[k] = -1;
                                        }
                                }
                        }
                        control.updateImage();
                }
                ThebaGUI control = ThebaGUI.getInstance();
                for (int z = 0; z < depth; z++) {
                        short[] data = control.getSlice(z);
                        control.setProgress(z);
                        if (control.isStopped())
                                return;
                        for (int x = 0; x < width; x++) {
                                for (int y = 0; y < height; y++) {
                                        int index = x + y * width;
                                        if (data[index] == 255)
                                                continue;
                                        if ((data[index] & (1 << 12)) != 0) {
                                                data[index] = (short) (data[index] & ~(1 << 12));
                                        } else {
                                                data[index] = 0;
                                        }
                                }
                        }
                }
                control.updateImage();
        }

        /*
         * Returns a factor describing the percentage of border-pixels that are
         * adjacent to a white border-area A too low percentage is a strong
         * indication that this region is not a valid fiber
         */
        public float checkReg(int xp, int yp, final short[] input) {
                short[] mark = new short[input.length];
                int width = control.getStack().getWidth();
                int height = control.getStack().getHeight();
                //int depth = control.getStack().getDepth();
                LinkedList<Point> queue = new LinkedList<Point>();
                queue.addFirst(new Point(xp, yp));
                float whiteCount = 0;
                float borderCount = 0;
                short col = input[xp + yp * width];
                while (queue.isEmpty() == false) {
                        Point p = queue.removeLast();
                        int x = p.x;
                        int y = p.y;
                        int index = x + y * width;
                        if (x + 1 < width) {
                                if (input[index + 1] == 255)
                                        whiteCount++;
                                else if (input[index + 1] == 0)
                                        borderCount++;
                                if (mark[index + 1] == 0 && input[index + 1] == col) {
                                        mark[index + 1] = 1;
                                        queue.addFirst(new Point(x + 1, y));
                                }
                        }
                        if (y + 1 < height) {
                                if (input[index + width] == 255)
                                        whiteCount++;
                                else if (input[index + width] == 0)
                                        borderCount++;
                                if (mark[index + width] == 0 && input[index + width] == col) {
                                        mark[index + width] = 1;
                                        queue.addFirst(new Point(x, y + 1));
                                }
                        }
                        if (y - 1 >= 0) {
                                if (input[index - width] == 255)
                                        whiteCount++;
                                else if (input[index - width] == 0)
                                        borderCount++;
                                if (mark[index - width] == 0 && input[index - width] == col) {
                                        mark[index - width] = 1;
                                        queue.addFirst(new Point(x, y - 1));
                                }
                        }
                        if (x - 1 >= 0) {
                                if (input[index - 1] == 255)
                                        whiteCount++;
                                else if (input[index - 1] == 0)
                                        borderCount++;
                                if (mark[index - 1] == 0 && input[index - 1] == col) {
                                        mark[index - 1] = 1;
                                        queue.addFirst(new Point(x - 1, y));
                                }
                        }
                }
                borderCount += whiteCount;
                return whiteCount / borderCount;
        }

        public void findCandidates(int slice) {
                int width = control.getStack().getWidth();
                int height = control.getStack().getHeight();
                //int depth = control.getStack().getDepth();
                short[] pixels = control.getSlice(slice);
                short[] copy = pixels.clone();
                for (int y = 0; y < height; y++) {
                        for (int x = 0; x < width; x++) {
                                int index = x + y * width;
                                if (copy[index] == 0) {
                                        short[] mask = new short[pixels.length];
                                        short color = (short) 0xf0;
                                        int lumenArea = ImageFunctions.floodFill2D(x, y, width,
                                                        height, copy, mask, color);
                                        for (int i = 0; i < pixels.length; i++)
                                                if (mask[i] != 0)
                                                        copy[i] = (short) 1;
                                        // Decide if this area is a ok lumen.
                                        boolean okLumen = true;
                                        okLumen = okLumen && lumenArea < 8000 && lumenArea > 50;
                                        RegionMask2D maskreg = new RegionMask2D(mask, width, height);
                                        // get bounds, use bounds
                                        if (okLumen) {
                                                double circularity = (new CircularityDescriptor())
                                                .measure(mask, width, height);
                                                okLumen = okLumen && circularity < 3.8
                                                && circularity > 0.7;
                                        }
                                        if (okLumen) {
                                                double curvature = (new AvgYoungCurvature()).measure(
                                                                mask, width, height);
                                                okLumen = okLumen && curvature > 0 && curvature < 1.6;
                                        }
                                        if (okLumen) {
                                                double bending_energy = (new YoungBendingEnergy())
                                                .measure(mask, width, height);
                                                okLumen = okLumen && bending_energy > 0
                                                && bending_energy < 8.0; // 4.0
                                        }
                                        if (okLumen) {
                                                double convexity = (Double) (new ConvexityDescriptor())
                                                .measure(maskreg);
                                                okLumen = okLumen && convexity > 0.8 && convexity < 1.0;
                                        }
                                        if (okLumen) {
                                                double solidity = (new SolidityDescriptor()).measure(
                                                                mask, width, height);
                                                okLumen = okLumen && solidity < 1.3 && solidity > 0.7;
                                        }
                                        if (okLumen) {
                                                double eccentricity = (new EccentricityDescriptor())
                                                .measure(mask, width, height);
                                                okLumen = okLumen && eccentricity < 15
                                                && eccentricity > 1; // 4
                                        }
                                        // if the area is very small and the eccentricity is large,
                                        // then we have a small thin region found in
                                        // ring-like fiber
                                        if (!okLumen && lumenArea < 200 && lumenArea > 80) {
                                                double eccentricity = (new EccentricityDescriptor())
                                                .measure(mask, width, height);
                                                okLumen = eccentricity > 15;
                                        }
                                        if (okLumen) {
                                                int internalRegionCount = (Integer) (new InternalRegionDescriptor())
                                                .measure(maskreg);
                                                okLumen = okLumen && internalRegionCount < 2;
                                        }
                                        if (okLumen) {
                                                okLumen = okLumen
                                                && !ImageFunctions.touchBorder(mask, width,
                                                                height);
                                        }
                                        // Ok, let's assume this is a lumen now
                                        if (okLumen) {
                                                LumenCandidate newSeed = new LumenCandidate(x, y,
                                                                slice, (lumenArea));
                                                seeds.add(newSeed);
                                                for (int i = 0; i < pixels.length; i++)
                                                        if (mask[i] != 0)
                                                                copy[i] = 0xcc;
                                        }
                                }
                        }
                }
                control.showImage(copy);
        }

        public LinkedList getFibers() {
                return fibers;
        }

        public boolean isReserved(short val) {
                if (val == 255 || val == control.INVALID)
                        return true;
                return false;
        }

        @Override
        public void mouseClicked(Point3D p) {
                if (seedButton.isSelected())
                        selectSeed(new Point(p.x, p.y));
        }

        public boolean removeCrack(int slice, short[] lumen) {
                int width = control.getStack().getWidth();
                int height = control.getStack().getHeight();
                int depth = control.getStack().getDepth();
                short[] data = control.getSlice(slice);
                short[] invdata = ImageFunctions.invertMask(data, false);
                ImageFunctions.invertMask(lumen, false);
                Rectangle bounds = ImageFunctions.getBounds(lumen, width, height, 12);
                // End search if fibre is outside volume (bounds==null indicates no
                // pixels found)
                if (bounds == null || bounds.x < 0 || bounds.y < 0
                                || bounds.width > width || bounds.height > height) {
                        seedStillValid = false;
                        slice = depth;
                        return true;
                }
                if (map == null) {
                        map = new short[data.length];
                        map2 = new short[data.length];
                }
                // Create a distancemap for use by backtracking, where lumen= distance 1
                for (int i = 0; i < map.length; i++) {
                        if (lumen[i] != 0)
                                map[i] = 1;
                        else
                                map[i] = 0;
                }
                map = ImageFunctions.dilate(lumen, map, width, height, bounds,
                                (short) 2);
                int maxDilations = control.getPreferences().getInt("maxDilations", 12);
                for (int i = 3; i < maxDilations; i++) {
                        System.arraycopy(map, 0, map2, 0, map.length);
                        map = ImageFunctions.dilateMasked(map2, map, invdata, width,
                                        height, bounds, (short) i);
                }
                map = ImageFunctions.mask(map, invdata, true);
                BackTrack.backTrack(map, 9, width, height, bounds);
                // copy new lumen back to source
                for (int i = 0; i < map.length; i++) {
                        lumen[i] = map[i];
                        if (lumen[i] != 0)
                                data[i] = currentId;
                }
                return false;
        }

        @SuppressWarnings("unchecked")
        public void removeStuff(int slice, short[] lumen) {
                int width = control.getStack().getWidth();
                int height = control.getStack().getHeight();
                //int depth = control.getStack().getDepth();
                int regId = 1;
                short[] mark = new short[width * height];
                short[] slicedata = control.getSlice(slice);
                ArrayList<LumenCandidate> liste = new ArrayList<LumenCandidate>();
                for (int x = 0; x < width; x++) {
                        for (int y = 0; y < height; y++) {
                                if (mark[x + y * width] == 0 && lumen[x + y * width] != 0) {
                                        regId++;
                                        int size = ImageFunctions.floodFill2D(x, y, width, height,
                                                        lumen, mark, (short) regId);
                                        LumenCandidate thisCandidate = new LumenCandidate(x, y, 0,
                                                        size);
                                        float borderRatio = checkReg(x, y, slicedata);
                                        if (borderRatio < 0.7) {
                                                // no pixels have a white border, so we remove it!
                                                log("Removed isolated region");
                                                ImageFunctions.floodFill2D(x, y, width, height,
                                                                slicedata, (short) 0);
                                                ImageFunctions.floodFill2D(x, y, width, height, lumen,
                                                                (short) 0);
                                        } else {
                                                liste.add(thisCandidate);
                                        }
                                }
                        }
                }
                int numberOfRegionsToKeep = 2;
                if (liste.size() > 1) {
                        Collections.sort(liste);
                        boolean useMR_rules = true;
                        if (useMR_rules) {
                                LumenCandidate A = liste.get(0);
                                LumenCandidate B = liste.get(1);
                                // bruk maritrunar-regler
                                if (A.getSize() > 5 * B.getSize()) {
                                        // System.out.println("Deleted small lumen-part.");
                                        numberOfRegionsToKeep--;
                                }
                        }
                        for (int i = numberOfRegionsToKeep; i < liste.size(); i++) {
                                LumenCandidate l = liste.remove(i);
                                ImageFunctions.floodFill2D(l.x, l.y, width, height, slicedata,
                                                (short) 0);
                                ImageFunctions.floodFill2D(l.x, l.y, width, height, lumen,
                                                (short) 0);
                        }
                }
        }

        public void selectSeed(Point e) {
                int width = control.getStack().getWidth();
                int height = control.getStack().getHeight();
                //int depth = control.getStack().getDepth();
                short[] data = control.getCurrentPixels();
                short[] filledRegion = new short[data.length];
                if (data[(int) e.getX() + (int) e.getY() * width] == 255)
                        return;
                int count = ImageFunctions.floodFill2D(data, filledRegion, (int) e
                                .getX(), (int) e.getY(), width, height, (short) 200);
                Point p = e.getLocation();
                short[] displayedImage = data.clone();
                if (p != null) {
                        // copy marked pixels to output
                        for (int i = 0; i < data.length; i++) {
                                if (filledRegion[i] != 0) {
                                        displayedImage[i] = filledRegion[i];
                                }
                        }
                        control.showImage(displayedImage);
                        if (count > 10000) {
                                int ans = JOptionPane
                                .showConfirmDialog(
                                                null,
                                                "This regions has "
                                                + count
                                                + " pixels \nAre you sure you want to add this region?",
                                                "Confirm seed", JOptionPane.YES_NO_OPTION);
                                if (ans == JOptionPane.NO_OPTION) {
                                        control.updateImage();
                                        return;
                                }
                        }
                        seeds.add(new LumenCandidate(p.x, p.y, control.currentSlice(),
                                        count));
                }
        }

        public void setFibers(LinkedList fibers) {
                this.fibers = fibers;
        }

        @SuppressWarnings("unchecked")
        public void traceWalls() {
                int width = control.getStack().getWidth();
                int height = control.getStack().getHeight();
                int depth = control.getStack().getDepth();
                control.clearChannels();
                SliceWriter lumenlog = new SliceWriter("lumens.raw", width, height,
                                depth, 5);
                control.addChannel(lumenlog, "Lumens");
                SliceWriter walls = new SliceWriter("walls.raw", width, height, depth,
                                5);
                control.addChannel(walls, "Walls");
                SliceWriter skeletonLog = new SliceWriter("skeleton.raw", width,
                                height, depth, 5);
                control.addChannel(skeletonLog, "Skeleton");
                keepTracking = true;
                for (int z = 0; z < depth && keepTracking; z++) {
                        long time = System.currentTimeMillis();
                        short[] pixels = control.getSlice(z).clone();
                        short[] map = ImageFunctions.distance2d(pixels, width, height);
                        short[] mark = new short[width * height];
                        short[] sout = new short[pixels.length];
                        ArrayList<Lumen> lumens = new ArrayList<Lumen>();
                        for (int x = 0; x < width; x++) {
                                for (int y = 0; y < height; y++) {
                                        int index = x + y * width;
                                        if (mark[index] == 0 && pixels[index] != 0
                                                        && !isReserved(pixels[index])) {
                                                short[] newLumenMask = new short[pixels.length];
                                                short fid = pixels[index];
                                                int count = 0;
                                                for (int i = 0; i < newLumenMask.length; i++) {
                                                        if (pixels[i] == fid) {
                                                                newLumenMask[i] = fid;
                                                                mark[i] = (short) 1;
                                                                count++;
                                                        }
                                                }
                                                Rectangle bounds = ImageFunctions.getBounds(
                                                                newLumenMask, width, height, 12);
                                                int regionID = 1;
                                                short[] mark2 = new short[width * height];
                                                ArrayList<FiberSlice> liste = new ArrayList<FiberSlice>();
                                                for (int xx = 0; xx < width; xx++) {
                                                        for (int yy = 0; yy < height; yy++) {
                                                                if (mark2[xx + yy * width] == 0
                                                                                && newLumenMask[xx + yy * width] == fid) {
                                                                        regionID++;
                                                                        ImageFunctions.floodFill2D(xx, yy, width,
                                                                                        height, newLumenMask, mark2,
                                                                                        (short) regionID);
                                                                        short[] lumenPart = new short[newLumenMask.length];
                                                                        int size = ImageFunctions.floodFill2D(xx,
                                                                                        yy, width, height, newLumenMask,
                                                                                        lumenPart, (short) regionID);
                                                                        FiberSlice storedLumenPart = new FiberSlice(
                                                                                        lumenPart, size, z, fid);
                                                                        liste.add(storedLumenPart);
                                                                }
                                                        }
                                                }
                                                // Find path between the two largest lumen parts
                                                if (liste.size() > 1) {
                                                        Collections.sort(liste);
                                                        // Get the two largest lumen parts
                                                        FiberSlice lumen1 = liste.remove(0);
                                                        FiberSlice lumen2 = liste.remove(0);
                                                        // find border of lumen 2
                                                        short[] lumenBorder = ImageFunctions.dilate(
                                                                        lumen2.lumendata.clone(), width, height);
                                                        lumenBorder = ImageFunctions.mask(lumenBorder,
                                                                        ImageFunctions.invertMask(lumen2.lumendata,
                                                                                        true), true);
                                                        // find geodesic distance from lumen 1
                                                        short[] geodist = new short[lumen1.lumendata.length];
                                                        geodist = ImageFunctions.geodist(lumen1.lumendata,
                                                                        pixels, width, height);
                                                        for (int i = 0; i < geodist.length; i++) {
                                                                if (geodist[i] == 0)
                                                                        geodist[i] = 255;
                                                        }
                                                        // find path that connects lumen1 with lumen2 and
                                                        // goes through non-void part in the image
                                                        DistancePath bestPath = new DistancePath();
                                                        short[] path = bestPath.getPath(geodist,
                                                                        lumenBorder, map, width, height);
                                                        for (int i = 0; i < path.length; i++) {
                                                                if (path[i] != 0) {
                                                                        mark[i] = fid; // MUST AVOID FINDING NEW
                                                                        // PIXELS AS NEW REGION
                                                                        newLumenMask[i] = fid;
                                                                        pixels[i] = fid;
                                                                }
                                                        }
                                                }
                                                // Calculate and store distance from lumen
                                                short[] dmap = ImageFunctions.distance2d(ImageFunctions
                                                                .invertMask(newLumenMask, false), width,
                                                                height, bounds);
                                                Lumen newLumen = new Lumen(newLumenMask, dmap, bounds,
                                                                count, fid);
                                                lumens.add(newLumen);
                                        }
                                }
                        }
                        // Copy lumens into a separate file (only used for visual inspection
                        // later)
                        short[] lumenOutput = new short[pixels.length];
                        for (int i = 0; i < lumenOutput.length; i++) {
                                // copy only non-white pixels to output
                                if (!isReserved((short) i)) {
                                        lumenOutput[i] = pixels[i];
                                }
                        }
                        lumenlog.putSlice(lumenOutput, z);
                        short[] waterShed = getWatershed(z);
                        for (int lum = 0; lum < lumens.size(); lum++) {
                                Lumen lumen = lumens.get(lum);
                                Rectangle bounds = lumen.getBounds();
                                int startX = bounds.x;
                                int stopX = bounds.width + bounds.x;
                                int startY = bounds.y;
                                int stopY = bounds.height + bounds.y;
                                short[] currentWall = waterShed.clone();
                                short[] searchMask = new short[width * height];
                                searchMask = ImageFunctions.dilate(lumen.getLumenMask(), width,
                                                height);
                                for (int q = 0; q < pixels.length; q++) {
                                        if (lumen.getLumenMask()[q] != 0)
                                                searchMask[q] = 0;
                                }
                                // Join regions inside lumens
                                for (int x = 0; x < width; x++) {
                                        for (int y = 0; y < height; y++) {
                                                int index = x + y * width;
                                                if (currentWall[index] == 0 && searchMask[index] != 0) {
                                                        ImageFunctions
                                                        .floodFill2D(currentWall.clone(),
                                                                        currentWall, x, y, width, height,
                                                                        (short) 1);
                                                }
                                        }
                                }
                                // Fill in watershed borders
                                for (int x = 0; x < pixels.length; x++) {
                                        if (currentWall[x] == (short) 255) {
                                                currentWall[x] = 0;
                                        }
                                }
                                // dilate to fill borders
                                currentWall = ImageFunctions.dilate(currentWall.clone(),
                                                currentWall, width, height, lumen.getBounds());
                                currentWall = ImageFunctions.dilate(currentWall.clone(),
                                                currentWall, width, height, lumen.getBounds());
                                currentWall = ImageFunctions.dilate(currentWall.clone(),
                                                currentWall, width, height, lumen.getBounds());
                                // Obtain the original pixels ( without connecting paths written
                                // to them)
                                pixels = control.getSlice(z);
                                // get median wallsize
                                short ff = (short) 0xff;
                                ArrayList<Integer> sizeList = new ArrayList<Integer>();
                                for (int x = startX; x < stopX; x++) {
                                        for (int y = startY; y < stopY; y++) {
                                                int i = x + y * width;
                                                if (pixels[i] == ff && lumen.getDistanceMap()[i] < 9) {
                                                        sizeList.add((int) map[i]);
                                                }
                                        }
                                }
                                Collections.sort(sizeList);
                                int medianWallSize = 0;
                                if (sizeList.size() == 0) {
                                        System.err.println("warning sizelist empty!");
                                        break;
                                }
                                medianWallSize = sizeList.get(sizeList.size() / 2);
                                medianWallSize = medianWallSize * 2 + 1;
                                int maxWallSize = (int) (medianWallSize * 1.5);
                                // Mask out erroneously dilated pixels
                                double avgx = 0;
                                double avgy = 0;
                                double count = 0;
                                for (int x = startX; x < stopX; x++) {
                                        for (int y = startY; y < stopY; y++) {
                                                int i = x + y * width;
                                                if (currentWall[i] != 0 && pixels[i] != (short) 255 // not
                                                                // a
                                                                // foreground
                                                                // pixel
                                                                || lumen.getDistanceMap()[i] > maxWallSize) { // outside
                                                        // max
                                                        // wallsize
                                                        currentWall[i] = 0;
                                                }
                                                // remove lumenmark from data
                                                if (lumen.getLumenMask()[i] != 0) {
                                                        count++;
                                                        avgx += x;
                                                        avgy += y;
                                                        pixels[i] = 0;
                                                }
                                                if (currentWall[i] != 0) {
                                                        lumen.getLumenMask()[i] = (short) 255;
                                                }
                                        }
                                }
                                sout[(int) (avgx / count) + (int) (avgy / count) * width] = lumen
                                .getId();
                        }
                        skeletonLog.putSlice(sout, z);
                        for (int i = 0; i < pixels.length; i++) {
                                int min = Integer.MAX_VALUE;
                                short minId = 0;
                                if (pixels[i] != 0) {
                                        for (int j = 0; j < lumens.size(); j++) {
                                                Lumen candidate = lumens.get(j);
                                                if (candidate.getLumenMask()[i] == (short) 255) {
                                                        int val = candidate.getDistanceMap()[i];
                                                        if (val < min) {
                                                                min = val;
                                                                minId = candidate.getId();
                                                        }
                                                }
                                        }
                                        if (minId != 0) {
                                                pixels[i] = minId; // set type to type of closest lumen
                                        }
                                }
                        }
                        // Copy finished pixels and remove tracked walls from data source
                        short[] wallOutput = new short[pixels.length];
                        for (int i = 0; i < wallOutput.length; i++) {
                                if (pixels[i] != control.INVALID && pixels[i] != (short) 255
                                                && pixels[i] != 0) {
                                        wallOutput[i] = pixels[i];
                                        // pixels[i] = 0;
                                }
                        }
                        walls.putSlice(wallOutput, z);
                        control.updateSlice(z);
                        System.out.println("Total time for slice " + z + " : "
                                        + (System.currentTimeMillis() - time) / 1000.0);
                } // end z-loop
                control.updateImage();
                walls.flush();
                lumenlog.flush();
        }

        private void checkMerges(Point[][] skel) {
                int width = control.getStack().getWidth();
                int height = control.getStack().getHeight();
                int depth = control.getStack().getDepth();
                int[] mergeList = new int[control.MAX_FIBERS];
                for (int z = 1; z < depth - 1; z++) {
                        short[] slice = control.getSlice(z);
                        short[] next = control.getSlice(z + 1);
                        short[] prev = control.getSlice(z - 1);
                        for (int x = 1; x < width - 1; x++) {
                                for (int y = 1; y < height - 1; y++) {
                                        int index = x + y * width;
                                        if (slice[index] == currentId) {
                                                // merge only fibres that touch from above or below
                                                // mergeList[slice[index + 1] ] = 1;
                                                // mergeList[slice[index - 1] ] = 1;
                                                // mergeList[slice[index + width] ] = 1;
                                                // mergeList[slice[index - width] ] = 1;
                                                mergeList[next[index]] = 1;
                                                mergeList[prev[index]] = 1;
                                        }
                                }
                        }
                }
                for (short color = 1; color < control.MAX_FIBERS; color++) {
                        if (isReserved(color))
                                continue;
                        if (mergeList[color] != 0 && color != currentId) {
                                double avgDist = 0;
                                double count = 0;
                                for (int i = 0; i < depth; i++) {
                                        if (skel[currentId][i] != null && skel[color][i] != null) {
                                                double a = skel[currentId][i].x - skel[color][i].x;
                                                double b = skel[currentId][i].y - skel[color][i].y;
                                                double dist = Math.sqrt(a * a + b * b);
                                                // System.out.println(dist);
                                                avgDist += dist;
                                                count++;
                                        }
                                }
                                if (count > 0)
                                        avgDist /= count;
                                System.out.println("Overlaps " + count);
                                System.out.println("Avgdist " + avgDist);
                                // Only merge regions of they have an average distance
                                // across all coinciding lumen sections below this threshhold
                                if (avgDist < 50 || count == 0) {
                                        System.out.println("Merging " + color + " with "
                                                        + (currentId));
                                        control.releaseFiberId(color);
                                        for (int z = 0; z < depth; z++) {
                                                short[] slice = control.getSlice(z);
                                                for (int x = 1; x < width - 1; x++) {
                                                        for (int y = 1; y < height - 1; y++) {
                                                                int index = x + y * width;
                                                                if (slice[index] == color) {
                                                                        slice[index] = currentId;
                                                                }
                                                        }
                                                }
                                                skel[currentId][z] = ImageFunctions.getAveragePoint(
                                                                slice, width, height, currentId);
                                                skel[color][z] = null;
                                        }
                                } else {
                                        System.out.println("Rejected merge!");
                                }
                        }
                }
        }

        private short[] getWatershed(int z) {
                int width = control.getStack().getWidth();
                int height = control.getStack().getHeight();
                //int depth = control.getStack().getDepth();
                short[] waterShed = control.getSlice(z).clone();
                // Filling tracked lumen with white
                for (int i = 0; i < waterShed.length; i++) {
                        if (waterShed[i] != 0) {
                                waterShed[i] = (short) 0xff;
                        }
                }
                // make distancemap
                EDM edm = new EDM();
                edm.run(waterShed, width, height);
                for (int i = 0; i < waterShed.length; i++) {
                        waterShed[i] = (short) (waterShed[i] >> 2);
                }
                // invert map
                for (int i = 0; i < waterShed.length; i++) {
                        waterShed[i] = (short) (255 - waterShed[i]);
                }
                // set pixels originally of value 0 (backgrund) to 0, because we don't
                // want the watershed to grow into the background
                for (int i = 0; i < waterShed.length; i++) {
                        if (control.getSlice(z)[i] == 0)
                                waterShed[i] = (short) 0;
                }
                Watershed_Algorithm wc = new Watershed_Algorithm();
                wc.run(waterShed, width, height);
                int regionID = 1;
                int maxCount = Integer.MIN_VALUE;
                Point p = null;
                short[] mark = new short[width * height];
                int countregs = 0;
                for (int xx = 0; xx < width; xx++) {
                        for (int yy = 0; yy < height; yy++) {
                                if (mark[xx + yy * width] == 0
                                                && waterShed[xx + yy * width] == 0) {
                                        countregs++;
                                        regionID++;
                                        int count = ImageFunctions.floodFill2D(xx, yy, width,
                                                        height, waterShed, mark, (short) 1);
                                        if (count > maxCount) {
                                                p = new Point(xx, yy);
                                                maxCount = count;
                                        }
                                }
                        }
                }
                ImageFunctions.floodFill2D(p.x, p.y, width, height, waterShed.clone(),
                                waterShed, (short) 255);
                return waterShed;
        }

        @Override
        public void reset() {
        }

        @Override
        public void stop() {
        }
}