TrakEM2_/src/main/java/ini/trakem2/imaging/BinaryInterpolation2D.java

   1 /** Released under the GPL. Use at your own risk.
   2  *
   3  * @author Albert Cardona, Johannes Schindelin
   4  */
   5 package ini.trakem2.imaging;
   6
   7 import java.util.concurrent.Callable;
   8 import java.util.concurrent.ExecutionException;
   9 import java.util.concurrent.ExecutorService;
  10 import java.util.concurrent.Executors;
  11 import java.util.concurrent.Future;
  12
  13 import mpicbg.imglib.algorithm.OutputAlgorithm;
  14 import mpicbg.imglib.container.array.ArrayContainerFactory;
  15 import mpicbg.imglib.cursor.Cursor;
  16 import mpicbg.imglib.cursor.LocalizableByDimCursor;
  17 import mpicbg.imglib.image.Image;
  18 import mpicbg.imglib.image.ImageFactory;
  19 import mpicbg.imglib.type.logic.BitType;
  20 import mpicbg.imglib.type.numeric.integer.IntType;
  21
  22
  23 /** Given two binary images of the same dimensions,
  24  * generate an interpolated image that sits somewhere
  25  * in between, as specified by the weight.
  26  *
  27  * For each binary image, the edges are found
  28  * and then each pixel is assigned a distance to the nearest edge.
  29  * Inside, distance values are positive; outside, negative.
  30  * Then both processed images are compared, and wherever
  31  * the weighted sum is larger than zero, the result image
  32  * gets a pixel set to true (or white, meaning inside).
  33  *
  34  * A weight of zero means that the first image is not present at all
  35  * in the interpolated image;
  36  * a weight of one means that the first image is present exclusively.
  37  *
  38  * The code was originally created by Johannes Schindelin
  39  * in the VIB's vib.BinaryInterpolator class, for ij.ImagePlus.
  40  */
  41 public class BinaryInterpolation2D implements OutputAlgorithm<BitType>
  42 {
  43
  44         final private Image<BitType> img1, img2;
  45         private float weight;
  46         private Image<BitType> interpolated;
  47         private String errorMessage;
  48         private IDT2D idt1, idt2;
  49
  50         public BinaryInterpolation2D(final Image<BitType> img1, final Image<BitType> img2, final float weight) {
  51                 this.img1 = img1;
  52                 this.img2 = img2;
  53                 this.weight = weight;
  54         }
  55
  56         /** NOT thread safe, stateful. */
  57         private final class IDT2D {
  58                 final Image<IntType> result;
  59                 final int w, h;
  60                 final int[] position = new int[2];
  61                 final LocalizableByDimCursor<BitType> csrc;
  62                 final LocalizableByDimCursor<IntType> cout;
  63
  64                 IDT2D(final Image<BitType> img) {
  65                         this.w = img.getDimension(0);
  66                         this.h = img.getDimension(1);
  67                         ImageFactory<IntType> f = new ImageFactory<IntType>(new IntType(), new ArrayContainerFactory());
  68                         this.result = f.createImage(new int[]{w, h});
  69
  70                         // Set all result pixels to infinity
  71                         final int infinity = (w + h) * 9;
  72                         for (final IntType v : this.result) {
  73                                 v.set(infinity);
  74                         }
  75
  76                         // init result pixels with those of the image:
  77                         this.csrc = img.createLocalizableByDimCursor();
  78                         this.cout = result.createLocalizableByDimCursor();
  79
  80                         int count = 0;
  81                         for (int y = 0; y < h; y++) {
  82                                 for (int x = 0; x < w; x++) {
  83                                         if (isBoundary(x, y)) {
  84                                                 setOutValueAt(x, y, 0);
  85                                                 count++;
  86                                         } else if (isJustOutside(x, y)) {
  87                                                 setOutValueAt(x, y, -1);
  88                                         }
  89                                 }
  90                         }
  91
  92                         if (count > 0) {
  93                                 propagate();
  94                         }
  95
  96                         csrc.close();
  97                         cout.close();
  98                 }
  99
 100                 private final void setPosition(final LocalizableByDimCursor<?> c, final int x, final int y) {
 101                         position[0] = x;
 102                         position[1] = y;
 103                         c.setPosition(position);
 104                 }
 105
 106                 private final void setOutValueAt(final int x, final int y, final int value) {
 107                         setPosition(cout, x, y);
 108                         cout.getType().set(value);
 109                 }
 110
 111                 private final int getSrcValueAt(final int x, final int y) {
 112                         setPosition(csrc, x, y);
 113                         return csrc.getType().get() ? 1 : 0;
 114                 }
 115
 116                 private final int getOutValueAt(final int x, final int y) {
 117                         setPosition(cout, x, y);
 118                         return cout.getType().get();
 119                 }
 120
 121                 // reads from result, writes to result
 122                 private final void idt(final int x, final int y, final int dx, final int dy) {
 123                         if (x + dx < 0 || y + dy < 0 ||
 124                                 x + dx >= w || y + dy >= h) {
 125                                 return;
 126                         }
 127                         int value = getOutValueAt(x + dx, y + dy);
 128                         final int distance = (dx == 0 || dy == 0 ? 3 : 4);
 129                         value += distance * (value < 0 ? -1 : 1);
 130                         setPosition(cout, x, y);
 131                         if (Math.abs(cout.getType().get()) > Math.abs(value)) {
 132                                 cout.getType().set(value);
 133                         }
 134                 }
 135
 136                 private final void propagate() {
 137                         for (int j = 0; j < h; j++)
 138                                 for (int i = 0; i < w; i++) {
 139                                         idt(i, j, -1, 0);
 140                                         idt(i, j, -1, -1);
 141                                         idt(i, j, 0, -1);
 142                                 }
 143
 144                         for (int j = h - 1; j >= 0; j--)
 145                                 for (int i = w - 1; i >= 0; i--) {
 146                                         idt(i, j, +1, 0);
 147                                         idt(i, j, +1, +1);
 148                                         idt(i, j, 0, +1);
 149                                 }
 150
 151                         for (int i = w - 1; i >= 0; i--)
 152                                 for (int j = h - 1; j >= 0; j--) {
 153                                         idt(i, j, +1, 0);
 154                                         idt(i, j, +1, +1);
 155                                         idt(i, j, 0, +1);
 156                                 }
 157
 158                         for (int i = 0; i < w; i++)
 159                                 for (int j = 0; j < h; j++) {
 160                                         idt(i, j, -1, 0);
 161                                         idt(i, j, -1, -1);
 162                                         idt(i, j, 0, -1);
 163                                 }
 164                 }
 165
 166                 private final boolean isBoundary(final int x, final int y) {
 167                         if (getSrcValueAt(x, y) == 0)
 168                                 return false;
 169                         if (x <= 0 || getSrcValueAt(x - 1, y) == 0)
 170                                 return true;
 171                         if (x >= w - 1 || getSrcValueAt(x + 1, y) == 0)
 172                                 return true;
 173                         if (y <= 0 || getSrcValueAt(x, y - 1) == 0)
 174                                 return true;
 175                         if (y >= h - 1 || getSrcValueAt(x, y + 1) == 0)
 176                                 return true;
 177                         if (x <= 0 || y <= 0 || getSrcValueAt(x - 1, y - 1) == 0)
 178                                 return true;
 179                         if (x <= 0 || y >= h - 1 || getSrcValueAt(x - 1, y + 1) == 0)
 180                                 return true;
 181                         if (x >= w - 1 || y <= 0 || getSrcValueAt(x + 1, y - 1) == 0)
 182                                 return true;
 183                         if (x >= w - 1 || y >= h - 1 || getSrcValueAt(x + 1, y + 1) == 0)
 184                                 return true;
 185                         return false;
 186                 }
 187
 188                 private final boolean isJustOutside(final int x, final int y) {
 189                         if (getSrcValueAt(x, y) != 0)
 190                                 return false;
 191                         if (x > 0 && getSrcValueAt(x - 1, y) != 0)
 192                                 return true;
 193                         if (x < w - 1 && getSrcValueAt(x + 1, y) != 0)
 194                                 return true;
 195                         if (y > 0 && getSrcValueAt(x, y - 1) != 0)
 196                                 return true;
 197                         if (y < h - 1 && getSrcValueAt(x, y + 1) != 0)
 198                                 return true;
 199                         if (x > 0 && y > 0 && getSrcValueAt(x - 1, y - 1) != 0)
 200                                 return true;
 201                         if (x > 0 && y < h - 1 && getSrcValueAt(x - 1, y + 1) != 0)
 202                                 return true;
 203                         if (x < w - 1 && y > 0 && getSrcValueAt(x + 1, y - 1) != 0)
 204                                 return true;
 205                         if (x < w - 1 && y < h - 1 && getSrcValueAt(x + 1, y + 1) != 0)
 206                                 return true;
 207                         return false;
 208                 }
 209         }
 210
 211         @Override
 212         public Image<BitType> getResult() {
 213                 return interpolated;
 214         }
 215
 216         @Override
 217         public boolean checkInput() {
 218                 if (img1.getNumDimensions() < 2 || img2.getNumDimensions() < 2) {
 219                         errorMessage = "Need at least 2 dimensions";
 220                         return false;
 221                 }
 222                 if (img1.getDimension(0) != img2.getDimension(0) || img1.getDimension(1) != img2.getDimension(1)) {
 223                         errorMessage = "Dimensions do not match";
 224                         return false;
 225                 }
 226                 if (weight < 0 || weight > 1) {
 227                         errorMessage = "Weight must be between 0 and 1, both inclusive.";
 228                         return false;
 229                 }
 230                 return true;
 231         }
 232
 233         @Override
 234         public String getErrorMessage() {
 235                 return errorMessage;
 236         }
 237
 238         private final class NewITD2D implements Callable<IDT2D> {
 239                 private final Image<BitType> img;
 240                 NewITD2D(final Image<BitType> img) {
 241                         this.img = img;
 242                 }
 243                 @Override
 244                 public IDT2D call() throws Exception {
 245                         return new IDT2D(img);
 246                 }
 247         }
 248
 249         /** After changing the weight, it's totally valid to call process() again,
 250          * and then getResult(). */
 251         public void setWeight(final float weight) throws IllegalArgumentException {
 252                 if (weight < 0 || weight > 1) {
 253                         throw new IllegalArgumentException("Weight must be between 0 and 1, both inclusive.");
 254                 }
 255                 this.weight = weight;
 256         }
 257
 258         @Override
 259         public boolean process() {
 260                 this.interpolated = process(this.weight);
 261                 return null != this.interpolated;
 262         }
 263
 264         /** The first time, it will prepare the distance transform images, which are computed only once. */
 265         public Image<BitType> process(final float weight)
 266         {
 267                 synchronized (this) {
 268                         if (null == idt1 || null == idt2) {
 269                                 ExecutorService exec = Executors.newFixedThreadPool(Math.min(2, Runtime.getRuntime().availableProcessors()));
 270                                 Future<IDT2D> fu1 = exec.submit(new NewITD2D(img1));
 271                                 Future<IDT2D> fu2 = exec.submit(new NewITD2D(img2));
 272                                 exec.shutdown();
 273
 274                                 try {
 275                                         this.idt1 = fu1.get();
 276                                         this.idt2 = fu2.get();
 277                                 } catch (InterruptedException ie) {
 278                                         throw new RuntimeException(ie);
 279                                 } catch (ExecutionException e) {
 280                                         throw new RuntimeException(e);
 281                                 }
 282                         }
 283                 }
 284
 285                 // Cannot just img1.createNewImage() because the container may not be able to receive data,
 286                 // such as the ShapeList container.
 287                 final ImageFactory<BitType> f = new ImageFactory<BitType>(new BitType(), new ArrayContainerFactory());
 288                 final Image<BitType> interpolated = f.createImage(new int[]{img1.getDimension(0), img1.getDimension(1)});
 289
 290                 if (img1.getContainer().compareStorageContainerCompatibility(img2.getContainer())) {
 291                         final Cursor<IntType> c1 = idt1.result.createCursor();
 292                         final Cursor<IntType> c2 = idt2.result.createCursor();
 293                         final Cursor<BitType> ci = interpolated.createCursor();
 294
 295                         while (ci.hasNext()) {
 296                                 c1.fwd();
 297                                 c2.fwd();
 298                                 ci.fwd();
 299
 300                                 if ((c1.getType().get() * weight) + (c2.getType().get() * (1 - weight)) > 0) {
 301                                         ci.getType().set(true);
 302                                 }
 303                         }
 304
 305                         c1.close();
 306                         c2.close();
 307                         ci.close();
 308                 } else {
 309                         System.out.println("using option 2");
 310                         final LocalizableByDimCursor<IntType> c1 = idt1.result.createLocalizableByDimCursor();
 311                         final LocalizableByDimCursor<IntType> c2 = idt2.result.createLocalizableByDimCursor();
 312                         final LocalizableByDimCursor<BitType> ci = interpolated.createLocalizableByDimCursor();
 313
 314                         while (ci.hasNext()) {
 315                                 ci.fwd();
 316                                 c1.setPosition(ci);
 317                                 c2.setPosition(ci);
 318
 319                                 if (0 <= c1.getType().get() * weight + c2.getType().get() * (1 - weight)) {
 320                                         ci.getType().set(true);
 321                                 }
 322                         }
 323
 324                         c1.close();
 325                         c2.close();
 326                         ci.close();
 327                 }
 328
 329                 return interpolated;
 330         }
 331 }