00001 package theba.core.math;
00002
00003 import ij.IJ;
00004 import ij.ImagePlus;
00005 import ij.ImageStack;
00006 import ij.gui.ImageWindow;
00007 import ij.process.ImageProcessor;
00008
00013 public class EDM {
00014
00015 ImagePlus imp;
00016
00017 String arg;
00018
00019 int maxEDM;
00020
00021 short[] xCoordinate, yCoordinate;
00022
00023 int[] levelStart;
00024
00025 int[] levelOffset;
00026
00027 int[] histogram;
00028
00029 int slice;
00030
00031 int count;
00032
00033 boolean watershed;
00034
00035 ImageWindow win;
00036
00037 boolean canceled;
00038
00039 ImageStack movie;
00040
00041 boolean debug = IJ.debugMode;
00042
00043 boolean invertImage;
00044
00045 static boolean whiteBackground = true;
00046
00047 public void run(short[] pixels, int width, int height) {
00048
00049 makeEDM(pixels, width, height);
00050
00051 }
00052
00059 public void makeEDM(short[] pixels, int width, int height) {
00060 int xmax, ymax;
00061 int offset, rowsize;
00062
00063 rowsize = width;
00064 xmax = width - 3;
00065 ymax = height - 3;
00066
00067 short[] image16 = new short[pixels.length];
00068
00069 for (int i = 0; i < pixels.length; i++)
00070 image16[i] = (short) (pixels[i] & 0xff);
00071
00072 for (int y = 0; y < height; y++) {
00073 for (int x = 0; x < width; x++) {
00074 offset = x + y * rowsize;
00075 if (image16[offset] > 0) {
00076 if ((x < 2) || (x > xmax) || (y < 2) || (y > ymax))
00077 setEdgeValue(offset, rowsize, image16, x, y, xmax, ymax);
00078 else
00079 setValue(offset, rowsize, image16);
00080 }
00081 }
00082 }
00083
00084 for (int y = height - 1; y >= 0; y--) {
00085 for (int x = width - 1; x >= 0; x--) {
00086 offset = x + y * rowsize;
00087 if (image16[offset] > 0) {
00088 if ((x < 2) || (x > xmax) || (y < 2) || (y > ymax))
00089 setEdgeValue(offset, rowsize, image16, x, y, xmax, ymax);
00090 else
00091 setValue(offset, rowsize, image16);
00092 }
00093 }
00094 }
00095
00096 int max = Integer.MIN_VALUE;
00097 int min = Integer.MAX_VALUE;
00098
00099 for (int i = 0; i < pixels.length; i++) {
00100 if ((image16[i] & 0xffff) > max)
00101 max = image16[i] & 0xffff;
00102 if ((image16[i] & 0xffff) < min)
00103 min = image16[i] & 0xffff;
00104 }
00105
00106 for (int i = 0; i < pixels.length; i++) {
00107 pixels[i] = (byte) ((int) (255.0 * ((double) (image16[i] & 0xffff) / (double) max)));
00108 }
00109
00110 }
00111
00112 void setValue(int offset, int rowsize, short[] image16) {
00113 final int one = 41;
00114 final int sqrt2 = 58;
00115 final int sqrt5 = 92;
00116 int v;
00117 int r1 = offset - rowsize - rowsize - 2;
00118 int r2 = r1 + rowsize;
00119 int r3 = r2 + rowsize;
00120 int r4 = r3 + rowsize;
00121 int r5 = r4 + rowsize;
00122 int min = 32767;
00123
00124 v = image16[r2 + 2] + one;
00125 if (v < min)
00126 min = v;
00127 v = image16[r3 + 1] + one;
00128 if (v < min)
00129 min = v;
00130 v = image16[r3 + 3] + one;
00131 if (v < min)
00132 min = v;
00133 v = image16[r4 + 2] + one;
00134 if (v < min)
00135 min = v;
00136
00137 v = image16[r2 + 1] + sqrt2;
00138 if (v < min)
00139 min = v;
00140 v = image16[r2 + 3] + sqrt2;
00141 if (v < min)
00142 min = v;
00143 v = image16[r4 + 1] + sqrt2;
00144 if (v < min)
00145 min = v;
00146 v = image16[r4 + 3] + sqrt2;
00147 if (v < min)
00148 min = v;
00149
00150 v = image16[r1 + 1] + sqrt5;
00151 if (v < min)
00152 min = v;
00153 v = image16[r1 + 3] + sqrt5;
00154 if (v < min)
00155 min = v;
00156 v = image16[r2 + 4] + sqrt5;
00157 if (v < min)
00158 min = v;
00159 v = image16[r4 + 4] + sqrt5;
00160 if (v < min)
00161 min = v;
00162 v = image16[r5 + 3] + sqrt5;
00163 if (v < min)
00164 min = v;
00165 v = image16[r5 + 1] + sqrt5;
00166 if (v < min)
00167 min = v;
00168 v = image16[r4] + sqrt5;
00169 if (v < min)
00170 min = v;
00171 v = image16[r2] + sqrt5;
00172 if (v < min)
00173 min = v;
00174
00175 image16[offset] = (short) min;
00176
00177 }
00178
00179 void setEdgeValue(int offset, int rowsize, short[] image16, int x, int y,
00180 int xmax, int ymax) {
00181 final int one = 41;
00182 final int sqrt2 = 58;
00183 final int sqrt5 = 92;
00184 int v;
00185 int r1 = offset - rowsize - rowsize - 2;
00186 int r2 = r1 + rowsize;
00187 int r3 = r2 + rowsize;
00188 int r4 = r3 + rowsize;
00189 int r5 = r4 + rowsize;
00190 int min = 32767;
00191 int offimage = image16[r3 + 2];
00192
00193 if (y < 2)
00194 v = offimage + one;
00195 else
00196 v = image16[r2 + 2] + one;
00197 if (v < min)
00198 min = v;
00199
00200 if (x < 2)
00201 v = offimage + one;
00202 else
00203 v = image16[r3 + 1] + one;
00204 if (v < min)
00205 min = v;
00206
00207 if (x > xmax)
00208 v = offimage + one;
00209 else
00210 v = image16[r3 + 3] + one;
00211 if (v < min)
00212 min = v;
00213
00214 if (y > ymax)
00215 v = offimage + one;
00216 else
00217 v = image16[r4 + 2] + one;
00218 if (v < min)
00219 min = v;
00220
00221 if ((x < 2) || (y < 2))
00222 v = offimage + sqrt2;
00223 else
00224 v = image16[r2 + 1] + sqrt2;
00225 if (v < min)
00226 min = v;
00227
00228 if ((x > xmax) || (y < 2))
00229 v = offimage + sqrt2;
00230 else
00231 v = image16[r2 + 3] + sqrt2;
00232 if (v < min)
00233 min = v;
00234
00235 if ((x < 2) || (y > ymax))
00236 v = offimage + sqrt2;
00237 else
00238 v = image16[r4 + 1] + sqrt2;
00239 if (v < min)
00240 min = v;
00241
00242 if ((x > xmax) || (y > ymax))
00243 v = offimage + sqrt2;
00244 else
00245 v = image16[r4 + 3] + sqrt2;
00246 if (v < min)
00247 min = v;
00248
00249 if ((x < 2) || (y < 2))
00250 v = offimage + sqrt5;
00251 else
00252 v = image16[r1 + 1] + sqrt5;
00253 if (v < min)
00254 min = v;
00255
00256 if ((x > xmax) || (y < 2))
00257 v = offimage + sqrt5;
00258 else
00259 v = image16[r1 + 3] + sqrt5;
00260 if (v < min)
00261 min = v;
00262
00263 if ((x > xmax) || (y < 2))
00264 v = offimage + sqrt5;
00265 else
00266 v = image16[r2 + 4] + sqrt5;
00267 if (v < min)
00268 min = v;
00269
00270 if ((x > xmax) || (y > ymax))
00271 v = offimage + sqrt5;
00272 else
00273 v = image16[r4 + 4] + sqrt5;
00274 if (v < min)
00275 min = v;
00276
00277 if ((x > xmax) || (y > ymax))
00278 v = offimage + sqrt5;
00279 else
00280 v = image16[r5 + 3] + sqrt5;
00281 if (v < min)
00282 min = v;
00283
00284 if ((x < 2) || (y > ymax))
00285 v = offimage + sqrt5;
00286 else
00287 v = image16[r5 + 1] + sqrt5;
00288 if (v < min)
00289 min = v;
00290
00291 if ((x < 2) || (y > ymax))
00292 v = offimage + sqrt5;
00293 else
00294 v = image16[r4] + sqrt5;
00295 if (v < min)
00296 min = v;
00297
00298 if ((x < 2) || (y < 2))
00299 v = offimage + sqrt5;
00300 else
00301 v = image16[r2] + sqrt5;
00302 if (v < min)
00303 min = v;
00304
00305 image16[offset] = (short) min;
00306
00307 }
00308
00309 void convertToBytes(int width, int height, short[] image16, byte[] image8) {
00310 final int one = 41;
00311 int v, offset;
00312 int round = one / 2;
00313
00314 maxEDM = 0;
00315 for (int y = 0; y < height; y++) {
00316 for (int x = 0; x < width; x++) {
00317 offset = x + y * width;
00318 v = (image16[offset] + round) / one;
00319 if (v > 255)
00320 v = 255;
00321 if (v > maxEDM)
00322 maxEDM = v;
00323 image8[offset] = (byte) v;
00324 }
00325 }
00326 }
00327
00332 public void findUltimatePoints(ImageProcessor ip) {
00333 int rowsize, offset, count, x, y;
00334 int CoordOffset, xmax, ymax;
00335 boolean setPixel;
00336
00337 IJ.showStatus("Finding ultimate points");
00338 if (debug) {
00339 movie = new ImageStack(ip.getWidth(), ip.getHeight());
00340 movie.addSlice("EDM", ip.duplicate());
00341 }
00342 if (watershed) {
00343 filterEDM(ip, true);
00344 filterEDM(ip, false);
00345 }
00346 if (debug)
00347 movie.addSlice("Filtered EDM", ip.duplicate());
00348 makeCoordinateArrays(ip);
00349 byte[] image = (byte[]) ip.getPixels();
00350 ImageProcessor ip2 = null;
00351 if (watershed)
00352 ip2 = ip.duplicate();
00353 int width = ip.getWidth();
00354 int height = ip.getHeight();
00355 rowsize = width;
00356 xmax = width - 1;
00357 ymax = height - 1;
00358 for (int level = maxEDM - 1; level >= 1; level--) {
00359 do {
00360 count = 0;
00361 for (int i = 0; i < histogram[level]; i++) {
00362 CoordOffset = levelStart[level] + i;
00363 x = xCoordinate[CoordOffset];
00364 y = yCoordinate[CoordOffset];
00365 offset = x + y * rowsize;
00366 if ((image[offset] & 255) != 255) {
00367 setPixel = false;
00368 if ((x > 0)
00369 && (y > 0)
00370 && ((image[offset - rowsize - 1] & 255) > level))
00371 setPixel = true;
00372 if ((y > 0)
00373 && ((image[offset - rowsize] & 255) > level))
00374 setPixel = true;
00375 if ((x < xmax)
00376 && (y > 0)
00377 && ((image[offset - rowsize + 1] & 255) > level))
00378 setPixel = true;
00379 if ((x < xmax) && ((image[offset + 1] & 255) > level))
00380 setPixel = true;
00381 if ((x < xmax)
00382 && (y < ymax)
00383 && ((image[offset + rowsize + 1] & 255) > level))
00384 setPixel = true;
00385 if ((y < ymax)
00386 && ((image[offset + rowsize] & 255) > level))
00387 setPixel = true;
00388 if ((x > 0)
00389 && (y < ymax)
00390 && ((image[offset + rowsize - 1] & 255) > level))
00391 setPixel = true;
00392 if ((x > 0) && ((image[offset - 1] & 255) > level))
00393 setPixel = true;
00394 if (setPixel) {
00395 image[offset] = (byte) 255;
00396 count++;
00397 }
00398 }
00399 }
00400 } while (count != 0);
00401 }
00402
00403 if (watershed) {
00404 byte[] image2 = (byte[]) ip2.getPixels();
00405 for (int i = 0; i < width * height; i++) {
00406 if (((image[i] & 255) > 0) && ((image[i] & 255) < 255))
00407 image2[i] = (byte) 0xff;
00408 }
00409 ip.insert(ip2, 0, 0);
00410 } else {
00411 for (int i = 0; i < width * height; i++) {
00412 if ((image[i] & 255) == 255)
00413 image[i] = (byte) 0;
00414 }
00415 }
00416 }
00417
00418 void filterEDM(ImageProcessor edm, boolean smooth) {
00419 int rowsize, offset;
00420 int xmax, ymax;
00421
00422 byte[] image = (byte[]) edm.getPixels();
00423 ImageProcessor ip2 = edm.duplicate();
00424 byte[] image2 = (byte[]) ip2.getPixels();
00425 int width = edm.getWidth();
00426 int height = edm.getHeight();
00427 rowsize = width;
00428 xmax = width - 1;
00429 ymax = height - 1;
00430 int p0, p1, p2, p3, p4, p5, p6, p7, p8;
00431 int v;
00432 for (int y = 0; y < height; y++) {
00433 for (int x = 0; x < width; x++) {
00434 offset = x + y * rowsize;
00435 p0 = image2[offset];
00436 if (p0 > 1) {
00437 p1 = x > 0 && y > 0 ? image2[offset - rowsize - 1] : get(
00438 x - 1, y - 1, image2, width, height);
00439 p2 = y > 0 ? image2[offset - rowsize] : get(x, y - 1,
00440 image2, width, height);
00441 p3 = x < xmax && y > 0 ? image2[offset - rowsize + 1]
00442 : get(x + 1, y - 1, image2, width, height);
00443 p4 = x < xmax ? image2[offset + 1] : get(x + 1, y, image2,
00444 width, height);
00445 p5 = x < xmax && y < ymax ? image2[offset + rowsize + 1]
00446 : get(x + 1, y + 1, image2, width, height);
00447 p6 = y < ymax ? image2[offset + rowsize] : get(x, y + 1,
00448 image2, width, height);
00449 p7 = x > 0 && y < ymax ? image2[offset + rowsize - 1]
00450 : get(x - 1, y + 1, image2, width, height);
00451 p8 = x > 0 ? image2[offset - 1] : get(x - 1, y, image2,
00452 width, height);
00453 v = p0 - 1;
00454 if (smooth)
00455 image[offset] = (byte) ((p0 + p1 + p2 + p3 + p4 + p5
00456 + p6 + p7 + p8) / 9);
00457 else {
00458 if (p2 == v
00459 && p4 == v
00460 && p6 == v
00461 && p8 == v
00462 && ((p1 == p0 && p3 == v && p5 == v && p7 == v)
00463 || (p3 == p0 && p1 == v && p5 == v && p7 == v)
00464 || (p5 == p0 && p1 == v && p3 == v && p7 == v) || (p7 == p0
00465 && p1 == v && p3 == v && p5 == v)))
00466 image[offset] = (byte) v;
00467 }
00468 }
00469 }
00470 }
00471 }
00472
00473 int get(int x, int y, byte[] pixels, int width, int height) {
00474 if (x <= 0)
00475 x = 0;
00476 if (x >= width)
00477 x = width - 1;
00478 if (y <= 0)
00479 y = 0;
00480 if (y >= height)
00481 y = height - 1;
00482 return pixels[x + y * width];
00483 }
00484
00491 void makeCoordinateArrays(ImageProcessor edm) {
00492 int rowsize, offset, v, ArraySize;
00493 int width = edm.getWidth();
00494 int height = edm.getHeight();
00495 histogram = edm.getHistogram();
00496
00497 ArraySize = 0;
00498 for (int i = 0; i < maxEDM - 1; i++)
00499 ArraySize += histogram[i];
00500 xCoordinate = new short[ArraySize];
00501 yCoordinate = new short[ArraySize];
00502 byte[] image = (byte[]) edm.getPixels();
00503 offset = 0;
00504 levelStart = new int[256];
00505 for (int i = 0; i < 256; i++) {
00506 levelStart[i] = offset;
00507 if ((i > 0) && (i < maxEDM))
00508 offset += histogram[i];
00509 }
00510
00511 levelOffset = new int[256];
00512 rowsize = width;
00513 for (int y = 0; y < height; y++) {
00514 for (int x = 0; x < width; x++) {
00515 v = image[x + y * rowsize] & 255;
00516 if ((v > 0) && (v < maxEDM)) {
00517 offset = levelStart[v] + levelOffset[v];
00518 xCoordinate[offset] = (short) x;
00519 yCoordinate[offset] = (short) y;
00520 levelOffset[v] += 1;
00521 }
00522 }
00523 }
00524 }
00525
00526
00527
00528
00529
00530
00531
00532
00533 void doWatershed(ImageProcessor ip1) {
00534 int[] table = makeFateTable();
00535 if (debug)
00536 movie.addSlice("EDM+UEPs", ip1.duplicate());
00537 IJ.showStatus("Watershed (press esc to cancel)");
00538 ImageProcessor ip2 = ip1.duplicate();
00539 for (int level = maxEDM - 1; level >= 1; level--) {
00540 IJ.showProgress(maxEDM - level, maxEDM - 1);
00541 do {
00542 count = 0;
00543 processLevel(level, 1, ip1, ip2, table);
00544 processLevel(level, 3, ip1, ip2, table);
00545 processLevel(level, 2, ip1, ip2, table);
00546 processLevel(level, 4, ip1, ip2, table);
00547 } while (count > 0);
00548 if (debug)
00549 movie.addSlice("level " + level, ip1.duplicate());
00550 if (win.running != true) {
00551 canceled = true;
00552 IJ.beep();
00553 break;
00554 }
00555 }
00556 if (!canceled)
00557 postProcess(ip1);
00558 if (debug) {
00559 movie.addSlice("Post-processed", ip1.duplicate());
00560 new ImagePlus("The movie", movie).show();
00561 }
00562 }
00563
00564 int[] makeFateTable() {
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574 int[] table =
00575
00576
00577 { 0, 0, 4, 4, 0, 0, 4, 4, 8, 0, 12, 12, 8, 0, 12, 12, 0, 0, 0, 0, 0, 0,
00578 0, 0, 8, 0, 12, 12, 8, 0, 12, 12, 1, 0, 0, 0, 0, 0, 0, 0, 8, 0,
00579 15, 15, 8, 0, 15, 15, 1, 0, 0, 0, 0, 0, 0, 0, 9, 0, 15, 15, 9,
00580 0, 15, 15, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
00581 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
00582 0, 0, 9, 0, 15, 15, 9, 0, 15, 15, 1, 0, 0, 0, 0, 0, 0, 0, 9, 0,
00583 15, 15, 9, 0, 15, 15, 2, 2, 6, 6, 0, 0, 6, 6, 0, 0, 15, 15, 0,
00584 0, 15, 15, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15, 15, 0, 0, 15, 15,
00585 1, 3, 15, 15, 0, 0, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 3,
00586 3, 15, 15, 0, 0, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 2, 2,
00587 6, 6, 0, 0, 6, 6, 0, 0, 15, 15, 0, 0, 15, 15, 0, 0, 0, 0, 0, 0,
00588 0, 0, 0, 0, 15, 15, 0, 0, 15, 15, 3, 3, 15, 15, 0, 0, 15, 15,
00589 15, 15, 15, 15, 15, 15, 15, 15, 3, 3, 15, 15, 0, 0, 15, 15, 15,
00590 15, 15, 15, 15, 15, 15, 15 };
00591 return table;
00592 }
00593
00594 void processLevel(int level, int pass, ImageProcessor ip1,
00595 ImageProcessor ip2, int[] table) {
00596 int rowSize = ip1.getWidth();
00597 int height = ip1.getHeight();
00598 int xmax = rowSize - 1;
00599 int ymax = ip1.getHeight() - 1;
00600 byte[] pixels1 = (byte[]) ip1.getPixels();
00601 byte[] pixels2 = (byte[]) ip2.getPixels();
00602 System.arraycopy(pixels1, 0, pixels2, 0, rowSize * height);
00603
00604 for (int i = 0; i < histogram[level]; i++) {
00605 int coordOffset = levelStart[level] + i;
00606 int x = xCoordinate[coordOffset];
00607 int y = yCoordinate[coordOffset];
00608 int offset = x + y * rowSize;
00609 int index;
00610 if ((pixels2[offset] & 255) != 255) {
00611 index = 0;
00612 if (x > 0 && y > 0
00613 && (pixels2[offset - rowSize - 1] & 255) == 255)
00614 index ^= 1;
00615 if (y > 0 && (pixels2[offset - rowSize] & 255) == 255)
00616 index ^= 2;
00617 if (x < xmax && y > 0
00618 && (pixels2[offset - rowSize + 1] & 255) == 255)
00619 index ^= 4;
00620 if (x < xmax && (pixels2[offset + 1] & 255) == 255)
00621 index ^= 8;
00622 if (x < xmax && y < ymax
00623 && (pixels2[offset + rowSize + 1] & 255) == 255)
00624 index ^= 16;
00625 if (y < ymax && (pixels2[offset + rowSize] & 255) == 255)
00626 index ^= 32;
00627 if (x > 0 && y < ymax
00628 && (pixels2[offset + rowSize - 1] & 255) == 255)
00629 index ^= 64;
00630 if (x > 0 && (pixels2[offset - 1] & 255) == 255)
00631 index ^= 128;
00632 switch (pass) {
00633 case 1:
00634 if ((table[index] & 1) == 1) {
00635 pixels1[offset] = (byte) 255;
00636 count++;
00637 }
00638 break;
00639 case 2:
00640 if ((table[index] & 2) == 2) {
00641 pixels1[offset] = (byte) 255;
00642 count++;
00643 }
00644 break;
00645 case 3:
00646 if ((table[index] & 4) == 4) {
00647 pixels1[offset] = (byte) 255;
00648 count++;
00649 }
00650 break;
00651 case 4:
00652 if ((table[index] & 8) == 8) {
00653 pixels1[offset] = (byte) 255;
00654 count++;
00655 }
00656 break;
00657 }
00658 }
00659 }
00660 }
00661
00662 void postProcess(ImageProcessor ip) {
00663 byte[] pixels = (byte[]) ip.getPixels();
00664 int size = ip.getWidth() * ip.getHeight();
00665 for (int i = 0; i < size; i++) {
00666 if ((pixels[i] & 255) < 255)
00667 pixels[i] = (byte) 0;
00668 }
00669 }
00670
00671 }
