ini/trakem2/vector/VectorString2D.java

   1 /**
   2
   3 TrakEM2 plugin for ImageJ(C).
   4 Copyright (C) 2006, 2007 Albert Cardona.
   5
   6 This program is free software; you can redistribute it and/or
   7 modify it under the terms of the GNU General Public License
   8 as published by the Free Software Foundation (http://www.gnu.org/licenses/gpl.txt )
   9
  10 This program is distributed in the hope that it will be useful,
  11 but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 GNU General Public License for more details.
  14
  15 You should have received a copy of the GNU General Public License
  16 along with this program; if not, write to the Free Software
  17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  18
  19 You may contact Albert Cardona at acardona at ini.phys.ethz.ch
  20 Institute of Neuroinformatics, University of Zurich / ETH, Switzerland.
  21 **/
  22
  23 package ini.trakem2.vector;
  24
  25 import ini.trakem2.utils.IJError;
  26 import ini.trakem2.utils.Utils;
  27 import ini.trakem2.utils.M;
  28 import ij.measure.Calibration;
  29
  30 /** String of vectors. */
  31 public class VectorString2D implements VectorString {
  32
  33         // all private to the package: accessible from Editions class
  34         private double[] x; // points
  35         private double[] y;
  36         private double[] v_x = null; // vectors, after resampling
  37         private double[] v_y = null;
  38         /** The length of the x,y and v_x, v_y resampled points (the actual arrays may be a bit longer) */
  39         private int length;
  40         /** The point interdistance after resampling. */
  41         private double delta = 0; // the delta used for resampling
  42         /** The Z coordinate of the entire planar curve represented by this VectorString2D. */
  43         private double z;
  44         private boolean closed;
  45
  46         private Calibration cal = null;
  47
  48         /** Construct a new String of Vectors from the given points and desired resampling point interdistance 'delta'. */
  49         public VectorString2D(double[] x, double[] y, double z, boolean closed) throws Exception {
  50                 if (!(x.length == y.length)) {
  51                         throw new Exception("x and y must be of the same length.");
  52                 }
  53                 this.length = x.length;
  54                 this.x = x;
  55                 this.y = y;
  56                 this.z = z;
  57                 this.closed = closed;
  58         }
  59
  60         /** Does NOT clone the vector arrays, which are initialized to NULL; only the x,y,delta,z. */
  61         public Object clone() {
  62                 double[] x2 = new double[length];
  63                 double[] y2 = new double[length];
  64                 System.arraycopy(x, 0, x2, 0, length);
  65                 System.arraycopy(y, 0, y2, 0, length);
  66                 try {
  67                         return new VectorString2D(x2, y2, z, closed);
  68                 } catch (Exception e) {
  69                         IJError.print(e);
  70                         return null;
  71                 }
  72         }
  73
  74         public int length() { return length; }
  75
  76         /** If not resampled, returns zero. */
  77         public double getDelta() { return delta; }
  78
  79         public double[] getPoints(final int dim) {
  80                 switch (dim) {
  81                         case 0: return x;
  82                         case 1: return y;
  83                 }
  84                 return null;
  85         }
  86         public double[] getVectors(final int dim) {
  87                 switch (dim) {
  88                         case 0: return v_x;
  89                         case 1: return v_y;
  90                 }
  91                 return null;
  92         }
  93         public double getPoint(final int dim, final int i) {
  94                 switch (dim) {
  95                         case 0: return x[i];
  96                         case 1: return y[i];
  97                         case 2: return z;
  98                 }
  99                 return 0;
 100         }
 101         public double getVector(final int dim, final int i) {
 102                 switch (dim) {
 103                         case 0: return v_x[i];
 104                         case 1: return v_y[i];
 105                 }
 106                 return 0;
 107         }
 108         public boolean isClosed() { return closed; }
 109
 110         public double getAverageDelta() { // equivalent to C's getAndSetAveragePointInterdistance function
 111                 double d = 0;
 112                 for (int i=x.length -1; i>0; i--) {
 113                         d += Math.sqrt( (x[i] - x[i-1])*(x[i] - x[i-1]) + (y[i] - y[i-1])*(y[i] - y[i-1])); // pytagorian rule for 2D
 114                 }
 115                 return d / x.length;
 116         }
 117
 118         /** Same as resample(delta). */
 119         public void resample(double delta, boolean with_source) {
 120                 resample(delta);
 121         }
 122
 123         /** Homogenize the average point interdistance to 'delta'. */ // There are problems with the last point.
 124         public void resample(double delta) {
 125                 if (Math.abs(delta - this.delta) < 0.0000001) {
 126                         // delta is the same
 127                         return;
 128                 }
 129                 this.delta = delta; // store for checking purposes
 130                 this.resample();
 131         }
 132
 133         /** As in the resampling method for CurveMorphing_just_C.c but for 2D. Uses the assigned 'delta'. Will reorder to counter-clock wise if necessary. */
 134         private void resample() {
 135                 final int MAX_AHEAD = 6;
 136                 final double MAX_DISTANCE = 2.5 * delta;
 137                 double[] ps_x = new double[length];
 138                 double[] ps_y = new double[length];
 139                 double[] v_x = new double[length];
 140                 double[] v_y = new double[length];
 141                 final int p_length = this.length; // to keep my head cool
 142                 int ps_length = this.length; // the length of the resampled vectors
 143
 144                 // reorder to CCW if needed: (so all curves have the same orientation)
 145                 // find bounding box:
 146                 double x_max = 0;                       int x_max_i = 0;
 147                 double y_max = 0;                       int y_max_i = 0;
 148                 double x_min = Double.MAX_VALUE;                int x_min_i = 0;
 149                 double y_min = Double.MAX_VALUE;                int y_min_i = 0;
 150                 for (int i=0;i <p_length; i++) {
 151                         if (x[i] > x_max) { x_max = x[i]; x_max_i = i; } // this lines could be optimized, the p->x etc. are catched below
 152                         if (y[i] > y_max) { y_max = y[i]; y_max_i = i; }
 153                         if (x[i] < x_min) { x_min = x[i]; x_min_i = i; }
 154                         if (y[i] < y_min) { y_min = y[i]; y_min_i = i; }
 155                 }
 156                 int collect = 0;
 157                 if (y_min_i - x_max_i >= 0) collect++;
 158                 if (x_min_i - y_min_i >= 0) collect++;
 159                 if (y_max_i - x_min_i >= 0) collect++;
 160                 if (x_max_i - y_max_i >= 0) collect++;
 161                 //if (3 == collect)
 162                 if (3 != collect) { // this should be '3 == collect', but then we are looking at the curves from the other side relative to what ImageJ is showing. In any case as long as one or the other gets rearranged, they'll be fine.
 163                         // Clockwise! Reorder to CCW by reversing the arrays in place
 164                         int n = p_length;
 165                         double tmp;
 166                         for (int i=0; i< p_length /2; i++) {
 167
 168                                 tmp = x[i];
 169                                 x[i] = x[n-i-1];
 170                                 x[n-i-1] = tmp;
 171
 172                                 tmp = y[i];
 173                                 y[i] = y[n-i-1];
 174                                 y[n-i-1] = tmp;
 175                         }
 176                 }
 177
 178                 // first resampled point is the same as 0
 179                 ps_x[0] = x[0];
 180                 ps_y[0] = y[0];
 181                 // the index over x,y
 182                 int i = 1;
 183                 // the index over ps (the resampled points)
 184                 int j = 1;
 185                 // some vars:
 186                 double dx, dy, sum;
 187                 double angle, angleXY, dist1, dist2;
 188                 int[] ahead = new int[MAX_AHEAD];
 189                 double[] distances = new double[MAX_AHEAD];
 190                 int n_ahead = 0;
 191                 int u, ii, k;
 192                 int t, s;
 193                 int prev_i = i;
 194                 int iu;
 195                 double dist_ahead;
 196                 double[] w = new double[MAX_AHEAD];
 197
 198                 // start infinite loop
 199                 for (;prev_i <= i;) {
 200                         if (prev_i > i) {//the loop has completed one round, since 'i' can go up to MAX_POINTS ahead of the last point into the points at the beggining of the array. Whenever the next point 'i' to start exploring is set beyond the length of the array, then the condition is met.
 201                                 break;
 202                         }
 203                         // check ps and v array lengths
 204                         if (j >= ps_length) {
 205                                 // must enlarge
 206                                 ps_length += 20;
 207                                 v_x = Utils.copy(v_x, ps_length);
 208                                 v_y = Utils.copy(v_y, ps_length);
 209                                 ps_x = Utils.copy(ps_x, ps_length);
 210                                 ps_y = Utils.copy(ps_y, ps_length);
 211                         }
 212                         // get distances of MAX_POINTs ahead from the previous point
 213                         n_ahead = 0; // reset
 214                         for (t=0; t<MAX_AHEAD; t++) {
 215                                 s = i + t; // 'i' is the first to start inspecting from
 216                                 // fix 's' if it goes over the end // TODO this is problematic for sure for open curves.
 217                                 if (s >= p_length) {
 218                                         if (this.closed) s = s - p_length;
 219                                         else s = p_length -1; // the last
 220                                 }
 221                                 dist_ahead = Math.sqrt((x[s] - ps_x[j-1])*(x[s] - ps_x[j-1]) + (y[s] - ps_y[j-1])*(y[s] - ps_y[j-1]));
 222                                 if (dist_ahead < MAX_DISTANCE) {
 223                                         ahead[n_ahead] = s;
 224                                         distances[n_ahead] = dist_ahead;
 225                                         n_ahead++;
 226                                 }
 227                         }
 228                         //
 229                         if (0 == n_ahead) { // no points ahead found within MAX_DISTANCE
 230                                 // all MAX_POINTS ahead lay under delta.
 231                                 // ...
 232                                 // simpler version: use just the next point
 233                                 dist1 = Math.sqrt((x[i] - ps_x[j-1])*(x[i] - ps_x[j-1]) + (y[i] - ps_y[j-1])*(y[i] - ps_y[j-1]));
 234                                 angleXY = M.getAngle(x[i] - ps_x[j-1], y[i] - ps_y[j-1]);
 235
 236                                 dx = Math.cos(angleXY) * delta;
 237                                 dy = Math.sin(angleXY) * delta;
 238                                 ps_x[j] = ps_x[j-1] + dx;
 239                                 ps_y[j] = ps_y[j-1] + dy;
 240                                 v_x[j] = dx;
 241                                 v_y[j] = dy;
 242
 243                                 //correct for point overtaking the not-close-enough point ahead in terms of 'delta_p' as it is represented in MAX_DISTANCE, but overtaken by the 'delta' used for subsampling:
 244                                 if (dist1 <= delta) {
 245                                         //look for a point ahead that is over distance delta from the previous j, so that it will lay ahead of the current j
 246                                         for (u=i; u<=p_length; u++) {
 247                                                 dist2 = Math.sqrt((x[u] - ps_x[j-1])*(x[u] - ps_x[j-1]) + (y[u] - ps_y[j-1])*(y[u] - ps_y[j-1]));
 248                                                 if (dist2 > delta) {
 249                                                         prev_i = i;
 250                                                         i = u;
 251                                                         break;
 252                                                 }
 253                                         }
 254                                 }
 255                         } else {
 256                                 w[0] = distances[0] / MAX_DISTANCE;
 257                                 double largest = w[0];
 258                                 for (u=1; u<n_ahead; u++) {
 259                                         w[u] = 1 - (distances[u] / MAX_DISTANCE);
 260                                         if (w[u] > largest) {
 261                                                 largest = w[u];
 262                                         }
 263                                 }
 264                                 // normalize weights: divide by largest
 265                                 sum = 0;
 266                                 for (u=0; u<n_ahead; u++) {
 267                                         w[u] = w[u] / largest;
 268                                         sum += w[u];
 269                                 }
 270                                 // correct error. The closest point gets the extra
 271                                 if (sum < 1.0) {
 272                                         w[0] += 1.0 - sum;
 273                                 } else {
 274                                         w = recalculate(w, n_ahead, sum);
 275                                 }
 276                                 // calculate the new point using the weights
 277                                 dx = 0.0;
 278                                 dy = 0.0;
 279                                 for (u=0; u<n_ahead; u++) {
 280                                         iu = i+u;
 281                                         if (iu >= p_length) iu -= p_length;
 282                                         angleXY = M.getAngle(x[iu] - ps_x[j-1], y[iu] - ps_y[j-1]);
 283                                         dx += w[u] * Math.cos(angleXY);
 284                                         dy += w[u] * Math.sin(angleXY);
 285                                 }
 286                                 // make vector and point:
 287                                 dx = dx * delta;
 288                                 dy = dy * delta;
 289                                 ps_x[j] = ps_x[j-1] + dx;
 290                                 ps_y[j] = ps_y[j-1] + dy;
 291                                 v_x[j] = dx;
 292                                 v_y[j] = dy;
 293
 294                                 // find the first point that is right ahead of the newly added point
 295                                 // so: loop through points that lay within MAX_DISTANCE, and find the first one that is right past delta.
 296                                 ii = i;
 297                                 for (k=0; k<n_ahead; k++) {
 298                                         if (distances[k] > delta) {
 299                                                 ii = ahead[k];
 300                                                 break;
 301                                         }
 302                                 }
 303                                 // correct for the case of unseen point (because all MAX_POINTS ahead lay under delta):
 304                                 prev_i = i;
 305                                 if (i == ii) {
 306                                         i = ahead[n_ahead-1] +1; //the one after the last.
 307                                         if (i >= p_length) i = i - p_length;
 308                                 } else {
 309                                         i = ii;
 310                                 }
 311                         }
 312                         //advance index in the new points
 313                         j += 1;
 314                 } // end of for (;;) loop
 315
 316                 // see whether the subsampling terminated too early, and fill with a line of points.
 317                 // TODO this is sort of a patch. Why didn't j overcome the last point is not resolved.
 318                 double lastx = x[0];
 319                 double lasty = y[0];
 320                 if (!closed) {
 321                         lastx = x[p_length -1];
 322                         lasty = y[p_length -1];
 323                 }
 324                 dist1 = Math.sqrt((lastx - ps_x[j-1])*(lastx - ps_x[j-1]) + (lasty - ps_y[j-1])*(lasty - ps_y[j-1]));
 325                 if (dist1 > delta*1.2) {
 326                         // TODO needs revision
 327                         // System.out.println("resampling terminated too early. Why?");
 328                         angleXY = M.getAngle(lastx - ps_x[j-1], lasty - ps_y[j-1]);
 329                         dx = Math.cos(angleXY) * delta;
 330                         dy = Math.sin(angleXY) * delta;
 331                         while (dist1 > delta*1.2) {//added 1.2 to prevent last point from being generated too close to the first point
 332                                 // check ps and v array lengths
 333                                 if (j >= ps_length) {
 334                                         // must enlarge.
 335                                         ps_length += 20;
 336                                         v_x = Utils.copy(v_x, ps_length);
 337                                         v_y = Utils.copy(v_y, ps_length);
 338                                         ps_x = Utils.copy(ps_x, ps_length);
 339                                         ps_y = Utils.copy(ps_y, ps_length);
 340                                 }
 341                                 //add a point
 342                                 ps_x[j] = ps_x[j-1] + dx;
 343                                 ps_y[j] = ps_y[j-1] + dy;
 344                                 v_x[j] = dx;
 345                                 v_y[j] = dy;
 346                                 j++;
 347                                 dist1 = Math.sqrt((lastx - ps_x[j-1])*(lastx - ps_x[j-1]) + (lasty - ps_y[j-1])*(lasty - ps_y[j-1]));
 348                         }
 349                 }
 350                 // set vector 0 to be the vector from the last point to the first // TODO also for non-closed?
 351                 angleXY = M.getAngle(lastx - ps_x[j-1], lasty - ps_y[j-1]);
 352                 //v_x[0] = Math.cos(angle) * delta;
 353                 //v_y[0] = Math.sin(angle) * delta; // can't use delta, it may be too long and thus overtake the first point!
 354                 v_x[0] = Math.cos(angleXY) * dist1;
 355                 v_y[0] = Math.sin(angleXY) * dist1;
 356
 357                 // assign the new subsampled points
 358                 this.x = ps_x;
 359                 this.y = ps_y;
 360                 // assign the vectors
 361                 this.v_x = v_x;
 362                 this.v_y = v_y;
 363                 // j is now the length of the ps_x, ps_y, v_x and v_y arrays (i.e. the number of resampled points).
 364                 this.length = j;
 365
 366                 // done!
 367         }
 368
 369         private final double[] recalculate(final double[] w, final int length, final double sum_) {
 370                 double sum = 0;
 371                 int q;
 372                 for (q=0; q<length; q++) {
 373                         w[q] = w[q] / sum_;
 374                         sum += w[q];
 375                 }
 376                 double error = 1.0 - sum;
 377                 // make it be an absolute value
 378                 if (error < 0.0) {
 379                         error = -error;
 380                 }
 381                 if (error < 0.005) {
 382                         w[0] += 1.0 - sum;
 383                 } else if (sum > 1.0) {
 384                         return recalculate(w, length, sum);
 385                 }
 386                 return w;
 387         }
 388
 389         public void reorder(final int new_zero) { // this function is optimized for speed: no array duplications beyond minimally necessary, and no superfluous method calls.
 390                 int i, j;
 391                 // copying
 392                 double[] tmp = new double[this.length];
 393                 double[] src;
 394                 // x
 395                 src = x;
 396                 for (i=0, j=new_zero; j<length; i++, j++) { tmp[i] = src[j]; }
 397                 for (j=0; j<new_zero; i++, j++) { tmp[i] = src[j]; }
 398                 x = tmp;
 399                 tmp = src;
 400                 // y
 401                 src = y;
 402                 for (i=0, j=new_zero; j<length; i++, j++) { tmp[i] = src[j]; }
 403                 for (j=0; j<new_zero; i++, j++) { tmp[i] = src[j]; }
 404                 y = tmp;
 405                 tmp = src;
 406                 // v_x
 407                 src = v_x;
 408                 for (i=0, j=new_zero; j<length; i++, j++) { tmp[i] = src[j]; }
 409                 for (j=0; j<new_zero; i++, j++) { tmp[i] = src[j]; }
 410                 v_x = tmp;
 411                 tmp = src;
 412                 // v_y
 413                 src = v_y;
 414                 for (i=0, j=new_zero; j<length; i++, j++) { tmp[i] = src[j]; }
 415                 for (j=0; j<new_zero; i++, j++) { tmp[i] = src[j]; }
 416                 v_y = tmp;
 417                 tmp = src;
 418
 419                 // equivalent would be:
 420                 //System.arraycopy(src, new_zero, tmp, 0, m - new_zero);
 421                 //System.arraycopy(src, 0, tmp, m - new_zero, new_zero);
 422                 //sv2.x = tmp;
 423                 //tmp = src;
 424
 425                 // release
 426                 tmp = null;
 427         }
 428
 429         /** Subtracts vs2 vector j to this vector i and returns its length, without changing any data. */
 430         public double getDiffVectorLength(final int i, final int j, final VectorString vs2) {
 431                 final VectorString2D vs = (VectorString2D)vs2;
 432                 final double dx = v_x[i] - vs.v_x[j];
 433                 final double dy = v_y[i] - vs.v_y[j];
 434                 return Math.sqrt(dx*dx + dy*dy);
 435         }
 436
 437         /** Distance from point i in this to point j in vs2. */
 438         public double distance(int i, VectorString vs, int j) {
 439                 VectorString2D vs2 = (VectorString2D)vs;
 440                 return distance(x[i], y[i],
 441                                 vs2.x[i], vs2.y[i]);
 442         }
 443
 444         static public double distance(double x1, double y1,
 445                                       double x2, double y2) {
 446                 return Math.sqrt(Math.pow(x1 - x2, 2)
 447                                + Math.pow(y1 - y2, 2));
 448         }
 449
 450         /** Create a new VectorString for the given range. If last &lt; first, it will be created as reversed. */
 451         public VectorString subVectorString(int first, int last) throws Exception {
 452                 boolean reverse = false;
 453                 if (last < first) {
 454                         int tmp = first;
 455                         first = last;
 456                         last = tmp;
 457                         reverse = true;
 458                 }
 459                 int len = last - first + 1;
 460                 double[] x = new double[len];
 461                 double[] y = new double[len];
 462                 System.arraycopy(this.x, first, x, 0, len);
 463                 System.arraycopy(this.y, first, y, 0, len);
 464                 final VectorString2D vs = new VectorString2D(x, y, this.z, this.closed);
 465                 if (reverse) vs.reverse();
 466                 if (null != this.v_x) {
 467                         // this is resampled, so:
 468                         vs.delta = this.delta;
 469                         // create vectors
 470                         vs.v_x = new double[len];
 471                         vs.v_y = new double[len];
 472                         for (int i=1; i<len; i++) {
 473                                 vs.v_x[i] = vs.x[i] - vs.x[i-1];
 474                                 vs.v_y[i] = vs.y[i] - vs.y[i-1];
 475                         }
 476                 }
 477                 return vs;
 478         }
 479
 480         /** Invert the order of points. Will clear all vector arrays if any! */
 481         public void reverse() {
 482                 Utils.reverse(x);
 483                 Utils.reverse(y);
 484                 delta = 0;
 485                 if (null != v_x) v_x = v_y = null;
 486         }
 487
 488         public int getDimensions() { return 2; }
 489
 490         /** Scale to match cal.pixelWidth, cal.pixelHeight and computed depth. If cal is null, returns immediately. Will make all vectors null, so you must call resample(delta) again after calibrating. So it brings all values to cal.units, such as microns. */
 491         public void calibrate(final Calibration cal) {
 492                 if (null == cal) return;
 493                 this.cal = cal;
 494                 final int sign = cal.pixelDepth < 0 ? - 1 :1;
 495                 for (int i=0; i<x.length; i++) {
 496                         x[i] *= cal.pixelWidth;
 497                         y[i] *= cal.pixelHeight; // should be the same as pixelWidth
 498                 }
 499                 z *= cal.pixelWidth * sign; // since layer Z is in pixels.
 500                 // reset vectors
 501                 if (null != v_x) v_x = v_y = null;
 502                 delta = 0;
 503         }
 504
 505         public boolean isCalibrated() {
 506                 return null != this.cal;
 507         }
 508         public Calibration getCalibrationCopy() {
 509                 return null == this.cal ? null : this.cal.copy();
 510         }
 511 }