Mise en œuvre de la segmentation des bassins versants en Java

J’essaie d’écrire ma propre implémentation de la segmentation des bassins versants pour un projet. J’ai une version qui renvoie quelque chose qui ressemble à la segmentation correcte en donnant des images vraiment sortingviales. Malheureusement, c’est super lent / inefficace et cela peut s’arrêter dans tous les cas.

Je travaille à partir de la description dans “Traitement d’images numériques” de Woods et Gonzales et de la page Wikipedia de Watershed. L’algorithme général est codé et inclus ci-dessous, mais j’ai l’impression que je boucle sur beaucoup de choses que je n’ai pas besoin d’être. J’apprécie toute aide ici, merci d’avance.

public static Set watershed(double[][] im) { //Get the Gradient Magnitude to use as our "topographic map." double t1 = System.currentTimeMillis(); double[][] edges = EdgeDetection.sobelMagnitude(im); //Only iterate over values in the gradient magnitude to speed up. double[] histogram = Image.getHistogram(edges); Image.drawHistogram(histogram, Color.red); int minPixelValue = 0; for (int i = 0; i  0) { minPixelValue = i; break; } } int h = im.length; int w = im[0].length; //SE is a 3x3 structuring element for morphological dilation. boolean[][] SE = {{true, true, true}, {true, true, true}, {true, true, true}}; //Keeping track of last iteration's components to see if two flooded together. ArrayList<Set> lastComponents = connectedComponents(getSet(EdgeDetection.underThreshold(edges, minPixelValue + 1))); ArrayList<Set> components; Set boundary = new HashSet(); for (int i = minPixelValue + 1; i < 256; i++) { if (histogram[i] != 0) { System.out.println("BEHHH " + i); t1 = System.currentTimeMillis(); ArrayList damLocations = new ArrayList(); HashMap<Integer, ArrayList> correspondingSets = new HashMap<Integer, ArrayList>(); //Figure out which of the old sets the new sets incorporated. //Here is where we check if sets flooded into eachother. //System.out.println("Checking for flooding"); components = connectedComponents(getSet(EdgeDetection.underThreshold(edges, i))); for (int nc = 0; nc < components.size(); nc++) { //System.out.println("Checking component " + nc); Set newComponent = components.get(nc); for (int oc = 0; oc < lastComponents.size(); oc++) { //System.out.println(" Against component " + oc); Set oldComponent = lastComponents.get(oc); if (numberInCommon(newComponent, oldComponent) > 0) { //System.out.println(" In there."); ArrayList oldSetsContained; if (correspondingSets.containsKey(nc)) { oldSetsContained = correspondingSets.get(nc); damLocations.add(nc); } else { //System.out.println(" Nope."); oldSetsContained = new ArrayList(); } oldSetsContained.add(oc); correspondingSets.put(nc, oldSetsContained); } } } System.out.println("Calculating overlapping sets: " + (System.currentTimeMillis() - t1)); //System.out.println("Check done."); for (int key : correspondingSets.keySet()) { Integer[] cs = new Integer[correspondingSets.get(key).size()]; correspondingSets.get(key).toArray(cs); if (cs.length == 1) { //System.out.println("Set " + cs[0] + " has grown without flooding."); } else { //System.out.println("The following sets have flooded together: " + Arrays.toSsortingng(cs)); } } //Build Damns to prevent flooding for (int c : damLocations) { System.out.println("Building dam for component " + c); Set bigComponent = components.get(c); System.out.println("Total size: " + bigComponent.size()); ArrayList<Set> littleComponent = new ArrayList<Set>(); for (int lcindex : correspondingSets.get(c)) { littleComponent.add(lastComponents.get(lcindex)); } Set unionSet = new HashSet(boundary); for (Set lc : littleComponent) { unionSet = union(unionSet, lc); } System.out.println("Building union sets: " + (System.currentTimeMillis() - t1)); while (intersection(unionSet, bigComponent).size() < bigComponent.size()) { for (int lIndex = 0; lIndex < littleComponent.size(); lIndex++) { Set lc = littleComponent.get(lIndex); Set lcBoundary = extractBoundaries(lc, SE, h, w); Set toAdd = new HashSet(); Set otherComponents = new HashSet(unionSet); otherComponents.removeAll(lc); otherComponents.removeAll(boundary); otherComponents = extractBoundaries(otherComponents, SE, h, w); for (Point pt : lcBoundary) { Set eightNbrs = get8Neighborhood(pt); for (Point nbr : eightNbrs) { if (bigComponent.contains(nbr) & !boundary.contains(nbr)) { Set neighborNbr = get8Neighborhood(nbr); if (intersection(neighborNbr, otherComponents).size() > 0) { boundary.add(nbr); edges[nbr.y][nbr.x] = 256; break; } else if (!lc.contains(nbr)) { toAdd.add(nbr); //if(i==65)System.out.println("Adding point " + nbr.y + " " + nbr.x); } else { //if(i==65)System.out.println("Already in here " + nbr.y + " " + nbr.x); } } } } t1 = System.currentTimeMillis(); lc.addAll(toAdd); toAdd.removeAll(toAdd); littleComponent.set(lIndex, lc); unionSet = new HashSet(boundary); for (Set ltc : littleComponent) { unionSet = union(unionSet, ltc); } System.out.println("This is a donk " + intersection(unionSet, bigComponent).size()); otherComponents = new HashSet(unionSet); otherComponents.removeAll(lc); otherComponents.removeAll(boundary); } } } } } boundary = close(boundary,h,w); Image.drawSet(boundary, h, w); return boundary; } 

    L’algorithme tel qu’il semble est au maximum O (n ^ 2). Vous avez beaucoup de boucles nestedes. Je n’ai pas pu trouver la description de Woods. Ce code de Christopher Mei implémente l’algorithme: il s’agit d’une implémentation très simple.

    WatershedPixel.java

     /* * Watershed algorithm * * Copyright (c) 2003 by Christopher Mei (christopher.mei@sophia.inria.fr) * * This plugin is free software; you can redissortingbute it and/or modify * it under the terms of the GNU General Public License version 2 * as published by the Free Software Foundation. * * This program is dissortingbuted in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this plugin; if not, write to the Free Software * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ import java.lang.*; import java.util.*; import ij.*; /** * The aim of WatershedPixel is to enable * sorting the pixels of an Image according * to their grayscale value. * * This is the first step of the Vincent * and Soille Watershed algorithm (1991) * **/ public class WatershedPixel implements Comparable { /** Value used to initialise the image */ final static int INIT = -1; /** Value used to indicate the new pixels that * are going to be processed (intial value * at each level) **/ final static int MASK = -2; /** Value indicating that the pixel belongs * to a watershed. **/ final static int WSHED = 0; /** Fictitious pixel **/ final static int FICTITIOUS = -3; /** x coordinate of the pixel **/ private int x; /** y coordinate of the pixel **/ private int y; /** grayscale value of the pixel **/ private byte height; /** Label used in the Watershed immersion algorithm **/ private int label; /** Distance used for working on pixels */ private int dist; /** Neighbours **/ private Vector neighbours; public WatershedPixel(int x, int y, byte height) { this.x = x; this.y = y; this.height = height; label = INIT; dist = 0; neighbours = new Vector(8); } public WatershedPixel() { label = FICTITIOUS; } public void addNeighbour(WatershedPixel neighbour) { /*IJ.write("In Pixel, adding :"); IJ.write(""+neighbour); IJ.write("Add done"); */ neighbours.add(neighbour); } public Vector getNeighbours() { return neighbours; } public Ssortingng toSsortingng() { return new Ssortingng("("+x+","+y+"), height : "+getIntHeight()+", label : "+label+", distance : "+dist); } public final byte getHeight() { return height; } public final int getIntHeight() { return (int) height&0xff; } public final int getX() { return x; } public final int getY() { return y; } /** Method to be able to use the Collections.sort static method. **/ public int compareTo(Object o) { if(!(o instanceof WatershedPixel)) throw new ClassCastException(); WatershedPixel obj = (WatershedPixel) o; if( obj.getIntHeight() < getIntHeight() ) return 1; if( obj.getIntHeight() > getIntHeight() ) return -1; return 0; } public void setLabel(int label) { this.label = label; } public void setLabelToINIT() { label = INIT; } public void setLabelToMASK() { label = MASK; } public void setLabelToWSHED() { label = WSHED; } public boolean isLabelINIT() { return label == INIT; } public boolean isLabelMASK() { return label == MASK; } public boolean isLabelWSHED() { return label == WSHED; } public int getLabel() { return label; } public void setDistance(int distance) { dist = distance; } public int getDistance() { return dist; } public boolean isFICTITIOUS() { return label == FICTITIOUS; } public boolean allNeighboursAreWSHED() { for(int i=0 ; i 

    Bassin hydrographiqueFIFO.java

     /* * Watershed plugin * * Copyright (c) 2003 by Christopher Mei (christopher.mei@sophia.inria.fr) * * This plugin is free software; you can redissortingbute it and/or modify * it under the terms of the GNU General Public License version 2 * as published by the Free Software Foundation. * * This program is dissortingbuted in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this plugin; if not, write to the Free Software * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ import java.util.*; import ij.*; /** This class implements a FIFO queue that * uses the same formalism as the Vincent * and Soille algorithm (1991) **/ public class WatershedFIFO { private LinkedList watershedFIFO; public WatershedFIFO() { watershedFIFO = new LinkedList(); } public void fifo_add(WatershedPixel p) { watershedFIFO.addFirst(p); } public WatershedPixel fifo_remove() { return (WatershedPixel) watershedFIFO.removeLast(); } public boolean fifo_empty() { return watershedFIFO.isEmpty(); } public void fifo_add_FICTITIOUS() { watershedFIFO.addFirst(new WatershedPixel()); } public Ssortingng toSsortingng() { SsortingngBuffer ret = new SsortingngBuffer(); for(int i=0; i 

    Bassin versantStructure.java

     /* * Watershed algorithm * * Copyright (c) 2003 by Christopher Mei (christopher.mei@sophia.inria.fr) * * This plugin is free software; you can redissortingbute it and/or modify * it under the terms of the GNU General Public License version 2 * as published by the Free Software Foundation. * * This program is dissortingbuted in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this plugin; if not, write to the Free Software * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ import java.lang.*; import java.util.*; import ij.process.*; import ij.*; import java.awt.*; /** * WatershedStructure contains the pixels * of the image ordered according to their * grayscale value with a direct access to their * neighbours. * **/ public class WatershedStructure { private Vector watershedStructure; public WatershedStructure(ImageProcessor ip) { byte[] pixels = (byte[])ip.getPixels(); Rectangle r = ip.getRoi(); int width = ip.getWidth(); int offset, topOffset, bottomOffset, i; watershedStructure = new Vector(r.width*r.height); /** The structure is filled with the pixels of the image. **/ for (int y=ry; y<(r.y+r.height); y++) { offset = y*width; IJ.showProgress(0.1+0.3*(yr.y)/(r.height)); for (int x=rx; x<(r.x+r.width); x++) { i = offset + x; int indiceY = yr.y; int indiceX = xr.x; watershedStructure.add(new WatershedPixel(indiceX, indiceY, pixels[i])); } } /** The WatershedPixels are then filled with the reference to their neighbours. **/ for (int y=0; y=0) currentPixel.addNeighbour((WatershedPixel)watershedStructure.get(x+1+bottomOffset)); if(y+1=0) { currentPixel.addNeighbour((WatershedPixel)watershedStructure.get(x-1+offset)); if(y-1>=0) currentPixel.addNeighbour((WatershedPixel)watershedStructure.get(x-1+bottomOffset)); if(y+1=0) currentPixel.addNeighbour((WatershedPixel)watershedStructure.get(x+bottomOffset)); if(y+1 

    Watershed_Algorithm.java

     /* * Watershed algorithm * * Copyright (c) 2003 by Christopher Mei (christopher.mei@sophia.inria.fr) * * This plugin is free software; you can redissortingbute it and/or modify * it under the terms of the GNU General Public License version 2 * as published by the Free Software Foundation. * * This program is dissortingbuted in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this plugin; if not, write to the Free Software * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ import ij.*; import ij.plugin.filter.PlugInFilter; import ij.process.*; import ij.gui.*; import ij.plugin.frame.PlugInFrame; import java.awt.*; import java.util.*; /** * This algorithm is an implementation of the watershed immersion algorithm * written by Vincent and Soille (1991). * * @Article{Vincent/Soille:1991, * author = "Lee Vincent and Pierre Soille", * year = "1991", * keywords = "IMAGE-PROC SKELETON SEGMENTATION GIS", * institution = "Harvard/Paris+Louvain", * title = "Watersheds in digital spaces: An efficient algorithm * based on immersion simulations", * journal = "IEEE PAMI, 1991", * volume = "13", * number = "6", * pages = "583--598", * annote = "Watershed lines (eg the continental divide) mark the * boundaries of catchment regions in a topographical map. * The height of a point on this map can have a direct * correlation to its pixel intensity. WIth this analogy, * the morphological operations of closing (or opening) * can be understood as smoothing the ridges (or filling * in the valleys). Develops a new algorithm for obtaining * the watershed lines in a graph, and then uses this in * developing a new segmentation approach based on the * {"}depth of immersion{"}.", * } * * A review of Watershed algorithms can be found at : * http://www.cs.rug.nl/~roe/publications/parwshed.pdf * * @Article{RoeMei00, * author = "Roerdink and Meijster", * title = "The Watershed Transform: Definitions, Algorithms and * Parallelization Strategies", * journal = "FUNDINF: Fundamenta Informatica", * volume = "41", * publisher = "IOS Press", * year = "2000", * } **/ public class Watershed_Algorithm implements PlugInFilter { private int threshold; final static int HMIN = 0; final static int HMAX = 256; public int setup(Ssortingng arg, ImagePlus imp) { if (arg.equals("about")) {showAbout(); return DONE;} return DOES_8G+DOES_STACKS+SUPPORTS_MASKING+NO_CHANGES; } public void run(ImageProcessor ip) { boolean debug = false; IJ.showStatus("Sorting pixels..."); IJ.showProgress(0.1); /** First step : the pixels are sorted according to increasing grey values **/ WatershedStructure watershedStructure = new WatershedStructure(ip); if(debug) IJ.write(""+watershedStructure.toSsortingng()); IJ.showProgress(0.8); IJ.showStatus("Start flooding..."); if(debug) IJ.write("Starting algorithm...\n"); /** Start flooding **/ WatershedFIFO queue = new WatershedFIFO(); int curlab = 0; int heightIndex1 = 0; int heightIndex2 = 0; for(int h=HMIN; h=0) {/*Initialise queue with neighbours at level h of current basins or watersheds*/ p.setDistance(1); queue.fifo_add(p); break; } // end if } // end for } // end for int curdist = 1; queue.fifo_add_FICTITIOUS(); while(true) /** extend basins **/{ WatershedPixel p = queue.fifo_remove(); if(p.isFICTITIOUS()) if(queue.fifo_empty()) break; else { queue.fifo_add_FICTITIOUS(); curdist++; p = queue.fifo_remove(); } if(debug) { IJ.write("\nWorking on :"); IJ.write(""+p); } Vector neighbours = p.getNeighbours(); for(int i=0 ; i0 || q.isLabelWSHED()) ) {*/ if( (q.getDistance() <= curdist) && (q.getLabel()>=0) ) { /* q belongs to an existing basin or to a watershed */ if(q.getLabel() > 0) { if( p.isLabelMASK() ) // Removed from original algorithm || p.isLabelWSHED() ) p.setLabel(q.getLabel()); else if(p.getLabel() != q.getLabel()) p.setLabelToWSHED(); } // end if lab>0 else if(p.isLabelMASK()) p.setLabelToWSHED(); } else if( q.isLabelMASK() && (q.getDistance() == 0) ) { if(debug) IJ.write("Adding value"); q.setDistance( curdist+1 ); queue.fifo_add( q ); } } // end for, end processing neighbours if(debug) { IJ.write("End processing neighbours"); IJ.write("New val :\n"+p); IJ.write("Queue :\n"+queue); } } // end while (loop) /* Detect and process new minima at level h */ for(int pixelIndex = heightIndex2 ; pixelIndex