ODEsolver.java

package com.github.celldynamics.quimp.plugin.ecmm;

import com.github.celldynamics.quimp.Vert;
import com.github.celldynamics.quimp.geom.ExtendedVector2d;
import com.github.celldynamics.quimp.utils.QuimpToolsCollection;

import ij.IJ;
import ij.process.ImageProcessor;

/**
 * ODE Solver.
 * 
 * @author rtyson
 *
 */
public class ODEsolver {

  private static boolean inside;

  /**
   * Default constructor.
   */
  public ODEsolver() {
  }

  /**
   * Euler solver.
   * 
   * @param v vertex
   * @param s sector
   * @return ?
   */
  public static ExtendedVector2d euler(Vert v, Sector s) {
    // Vect2d[] history = new Vect2d[ECMp.maxIter];
    int x;
    int y;
    int lastSampleX = -1;
    int lastSampleY = -1; // store where last sample was
    double dist = 0; // distance migrated
    double tempFlu;
    Vert edge;
    ExtendedVector2d p;
    ExtendedVector2d pp;

    v.snapped = false;

    if (ECMp.ANA) { // sample at boundary
      p = v.getPoint();
      x = (int) Math.round(p.getX());
      y = (int) Math.round(p.getY());
      lastSampleX = x;
      lastSampleY = y;
      tempFlu = ODEsolver.sampleFluo(ECMp.image, x, y);
      v.fluores[0].intensity = tempFlu;
      v.fluores[0].x = x; // store in first slot
      v.fluores[0].y = y;
    }

    if (ECMp.plot) {
      ECMM_Mapping.plot.setColor(0, 0, 0);
    }

    p = new ExtendedVector2d(v.getX(), v.getY());
    pp = new ExtendedVector2d(v.getX(), v.getY()); // previouse position

    // history[0] = new Vect2d(p.getX(), p.getY());

    boolean maxHit = false;
    int i = 1;
    ExtendedVector2d k;

    for (; i < ECMp.maxIter - 1; i++) {
      // IJ.log("\tIt " + i); //debug
      if (ODEsolver.proximity(p, s) || (ECMp.ANA && dist >= (ECMp.anaMigDist)) || maxHit) {
        // stop when within d of the target segment or
        // if migrated more than the ana set cortex width (in pixels)
        pp.setX(p.getX());
        pp.setY(p.getY());

        // if(!ECMp.ANA) { // no need to snap ana result. landing coord
        // not needed
        edge = ODEsolver.snap(p, s);
        dist += ExtendedVector2d.lengthP2P(pp, p);
        v.distance = QuimpToolsCollection.speedToScale(dist, ECMp.scale, ECMp.frameInterval);
        // if (s.expanding && !ECMp.ANA) {
        v.setLandingCoord(p, edge);
        // }
        // }

        if (ECMp.plot && ECMp.drawPaths) {
          ECMM_Mapping.plot.setColor(0, 0, 0);
          ECMM_Mapping.plot.drawLine(pp, p);
        }

        v.snapped = true;
        // System.out.println("iterations: " + i);
        break;
      }

      k = ODEsolver.dydt(p, s);
      k.multiply(ECMp.h);

      pp.setX(p.getX());
      pp.setY(p.getY());
      p.setX(p.getX() + k.getX());
      p.setY(p.getY() + k.getY());
      dist += ExtendedVector2d.lengthP2P(pp, p);

      if (ECMp.plot && ECMp.drawPaths) {
        // ECMM_Mapping.plot.setColor(1, 0, 0);
        ECMM_Mapping.plot.drawLine(pp, p);
      }
      // history[i] = new Vect2d(p.getX(), p.getY());

      if (ECMp.ANA) { // sample
        x = (int) Math.round(p.getX());
        y = (int) Math.round(p.getY());
        if (!(x == lastSampleX && y == lastSampleY)) { // on sample new locations
          lastSampleX = x;
          lastSampleY = y;
          tempFlu = ODEsolver.sampleFluo(ECMp.image, x, y);

          if (tempFlu > v.fluores[0].intensity) { // store first one
            // if((tempFlu / v.fluores[0].intensity)<1.1){
            // maxHit = true;
            // }
            v.fluores[0].intensity = tempFlu;
            v.fluores[0].x = x;
            v.fluores[0].y = y;

          }
        }
      }

      ECMp.its++;
    }

    if (ECMp.plot && !v.snapped && ECMp.drawFails) { // mark the start point of failed nodes
      ECMM_Mapping.plot.setColor(1, 0, 0);
      // p.print(v.getTrackNum() + "p: ");
      // pp.print(v.getTrackNum() + "pp: ");
      ECMM_Mapping.plot.drawCircle(v.getPoint(), 5);
    }

    return p;
  }

  /**
   * Get first derivative.
   * 
   * @param p ExtendedVector2d
   * @param s Sector
   * @return dy/dt
   */
  public static ExtendedVector2d dydt(ExtendedVector2d p, Sector s) {
    ExtendedVector2d result = fieldAt(p, s);
    result.multiply(ECMp.mobileQ);

    if (true) { // Math.abs(result.length()) > ECMp.maxVertF) {
      // IJ.log("!WARNING-max force exceeded: " +
      // Math.abs(result.length()));
      result.makeUnit();
      result.multiply(ECMp.maxVertF);
    }
    return result;
  }

  /**
   * Proximity.
   * 
   * @param p ExtendedVector2d
   * @param s Sector
   * @return ?
   */
  public static boolean proximity(ExtendedVector2d p, Sector s) {
    // could test against the chrages or the actual contour.
    // if using polar lines can use actual contour
    // Vert v = s.getTarStart();
    // if(true) return false;
    // Vert v = s.tarCharges.getHead();
    Vert v = s.getTarStart();
    do {
      double d = ExtendedVector2d.distPointToSegment(p, v.getPoint(), v.getNext().getPoint());
      // IJ.log("\t\tprox to: " + d); //debug
      if (d <= ECMp.d) {
        return true;
      }
      v = v.getNext();
    } while (!v.isIntPoint());
    return false;
  }

  private static Vert snap(ExtendedVector2d p, Sector s) {
    // snap p to the closest segment of target contour
    ExtendedVector2d current;
    Vert closestEdge;
    double distance;// = ECMp.d + 1; // must be closer then d+1, good starting value
    double tempDis;

    Vert v = s.getTarStart().getPrev(); // include the edge to the starting intersect pount
    distance = ExtendedVector2d.distPointToSegment(p, v.getPoint(), v.getNext().getPoint());
    v = v.getNext();
    closestEdge = v;
    do {
      current = ExtendedVector2d.pointToSegment(p, v.getPoint(), v.getNext().getPoint());
      tempDis = ExtendedVector2d.lengthP2P(p, current);

      if (tempDis < distance) {
        closestEdge = v;
        distance = tempDis;
      }
      v = v.getNext();
    } while (!v.isIntPoint());

    // p.setX(closest.getX());
    // p.setY(closest.getY());

    return closestEdge;
  }

  private static ExtendedVector2d fieldAt(ExtendedVector2d p, Sector s) {

    // Use line charges or point charges. remove if for speed
    // return fieldAtLines(p, s);
    if (ECMp.lineCharges) {
      return fieldAtLines(p, s);
    } else {
      return fieldAtPoints(p, s);
    }
  }

  private static ExtendedVector2d fieldAtPoints(ExtendedVector2d p, Sector s) {
    // calc the field size at p according to to migrating and target charges
    ExtendedVector2d field = new ExtendedVector2d(0, 0);
    ExtendedVector2d totalF = new ExtendedVector2d(0, 0);

    Vert v = s.migCharges.getHead();
    do {

      forceP(field, p, v.getPoint(), ECMp.migQ, ECMp.migPower);
      totalF.addVec(field);
      // totalF.print("\ttotlaF = ");
      v = v.getNext();
    } while (!v.getPrev().isIntPoint() || v.getPrev().isHead());

    v = s.tarCharges.getHead();
    do {
      forceP(field, p, v.getPoint(), ECMp.tarQ, ECMp.tarPower);
      totalF.addVec(field);
      v = v.getNext();
    } while (!v.getPrev().isIntPoint() || v.getPrev().isHead());

    return totalF;
  }

  private static void forceP(ExtendedVector2d force, ExtendedVector2d p, ExtendedVector2d pq,
          double q, double power) {
    double r = ExtendedVector2d.lengthP2P(pq, p);
    // System.out.println("\t r = " + r);
    if (r == 0) {
      force.setX(250);
      force.setY(250);
      IJ.log("!WARNING-FORCE INFINITE");
      return;
    }
    r = Math.abs(Math.pow(r, power));
    ExtendedVector2d unitV = ExtendedVector2d.unitVector(pq, p);
    double multiplier = (ECMp.k * (q / r));
    force.setX(unitV.getX() * multiplier);
    force.setY(unitV.getY() * multiplier);
  }

  private static ExtendedVector2d fieldAtLines(ExtendedVector2d p, Sector s) {
    // calc the field size at p according to to migrating and target charges
    ExtendedVector2d field = new ExtendedVector2d(0, 0);
    ExtendedVector2d totalF = new ExtendedVector2d(0, 0);
    double polarDir;

    // inside or outside sector?
    inside = s.insideCharges(p);

    if (!inside) {
      polarDir = -1;
      // System.out.println("switched");
    } else {
      polarDir = 1;
    }

    Vert v = s.migCharges.getHead();
    do {

      // forceL(field, p, v.getPoint(), v.getNext().getPoint(),
      // ECMp.migQ);

      /*
       * //times by the outerDirection to make lines polar sideDis =
       * Vect2d.distPoinToInfLine(p, v.getPoint(), v.getNext().getPoint()); if (sideDis < 0) {
       * polarDir = s.outerDirection * -1; } else { polarDir = s.outerDirection; }
       *
       */

      forceLpolar(field, p, v.getPoint(), v.getNext().getPoint(), ECMp.migQ, ECMp.migPower,
              polarDir);

      totalF.addVec(field);
      v = v.getNext();
    } while (!v.isIntPoint() || v.isHead());

    v = s.tarCharges.getHead();
    do {
      // forceL(field, p, v.getPoint(), v.getNext().getPoint(),
      // ECMp.tarQ);

      /*
       * sideDis = Vect2d.distPoinToInfLine(p, v.getPoint(), v.getNext().getPoint()); if
       * (sideDis < 0) { polarDir = s.outerDirection ; } else { polarDir = s.outerDirection *
       * -1; }
       *
       */

      forceLpolar(field, p, v.getPoint(), v.getNext().getPoint(), ECMp.tarQ, ECMp.tarPower,
              polarDir);

      totalF.addVec(field);
      v = v.getNext();
    } while (!v.isIntPoint() || v.isHead());

    return totalF;
  }

  private static void forceLpolar(ExtendedVector2d force, ExtendedVector2d p, ExtendedVector2d s1,
          ExtendedVector2d s2, double q, double power, double orientation) {
    double l = ExtendedVector2d.lengthP2P(s1, s2);
    ExtendedVector2d ru = ExtendedVector2d.unitVector(s2, p);
    double r = ExtendedVector2d.lengthP2P(s2, p);
    ExtendedVector2d rpU = ExtendedVector2d.unitVector(s1, p);
    double rp = ExtendedVector2d.lengthP2P(s1, p);

    double d = (((rp + r) * (rp + r)) - (l * l)) / (2 * l);
    // double d = ( Math.pow((rp + r), power) - (L * L)) / (2 * L);
    double multiplier = ((ECMp.k * q) / d);
    rpU.addVec(ru);

    force.setX(rpU.getX() * multiplier * orientation);
    force.setY(rpU.getY() * multiplier * orientation);
  }

  /**
   * Get intensity value for given coordinates.
   * 
   * @param ip image to sample intensities
   * @param x screen coordinate
   * @param y screen coordinate
   * @return mean intensity within 9-point stencil located at x,y
   */
  public static double sampleFluo(ImageProcessor ip, int x, int y) {
    double tempFlu = ip.getPixelValue(x, y) + ip.getPixelValue(x - 1, y)
            + ip.getPixelValue(x + 1, y) + ip.getPixelValue(x, y - 1) + ip.getPixelValue(x, y + 1)
            + ip.getPixelValue(x - 1, y - 1) + ip.getPixelValue(x + 1, y + 1)
            + ip.getPixelValue(x + 1, y - 1) + ip.getPixelValue(x - 1, y + 1);
    tempFlu = tempFlu / 9d;
    return tempFlu;
  }
}