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/../../../../com/github/celldynamics/quimp/geom/ExtendedVector2d.html#ExtendedVector2d">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/../../../com/github/celldynamics/quimp/geom/ExtendedVector2d.html#ExtendedVector2d">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/../com/github/celldynamics/quimp/geom/ExtendedVector2d.html#ExtendedVector2d">ExtendedVector2d fieldAtPoints(ExtendedVector2d p, Sector s) {
     // calc the field size at p according to to migrating and target charges
     ExtendedVector2deom/ExtendedVector2d.html#ExtendedVector2d">ExtendedVector2d field = new ExtendedVector2d(0, 0);
     ExtendedVector2dom/ExtendedVector2d.html#ExtendedVector2d">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../../../../../com/github/celldynamics/quimp/geom/ExtendedVector2d.html#ExtendedVector2d">ExtendedVector2d./../../../com/github/celldynamics/quimp/geom/ExtendedVector2d.html#ExtendedVector2d">ExtendedVector2d force, ExtendedVector2d../../../../../com/github/celldynamics/quimp/geom/ExtendedVector2d.html#ExtendedVector2d">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./../com/github/celldynamics/quimp/geom/ExtendedVector2d.html#ExtendedVector2d">ExtendedVector2d fieldAtLines(ExtendedVector2d p, Sector s) {
     // calc the field size at p according to to migrating and target charges
     ExtendedVector2deom/ExtendedVector2d.html#ExtendedVector2d">ExtendedVector2d field = new ExtendedVector2d(0, 0);
     ExtendedVector2dom/ExtendedVector2d.html#ExtendedVector2d">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../../../../../com/github/celldynamics/quimp/geom/ExtendedVector2d.html#ExtendedVector2d">ExtendedVector2d./../../../com/github/celldynamics/quimp/geom/ExtendedVector2d.html#ExtendedVector2d">ExtendedVector2d force, ExtendedVector2d../../../../../com/github/celldynamics/quimp/geom/ExtendedVector2d.html#ExtendedVector2d">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;
   }
 }