Constrictor.java

package com.github.celldynamics.quimp;

import java.io.FileWriter;
import java.io.PrintWriter;

import com.github.celldynamics.quimp.geom.ExtendedVector2d;

import ij.process.ImageProcessor;

/**
 * Calculate forces that affect the snake. Active Contour segmentation.
 * 
 * <p>This procedure is aware of possible overlapping and try to counteract them.
 * 
 * @author rtyson
 */
public class Constrictor {

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

  /**
   * Compute force power and moves nodes by predefined step.
   * 
   * <p>This routine should be called in loop. Each call moves nodes by
   * {@link BOAState.BOAp#delta_t}.
   * 
   * @param snake Processed snake
   * @param ip Original image
   * @return status of snake (true if it is frozen)
   * @see BOA_#tightenSnake
   */
  public boolean constrict(final Snake snake, final ImageProcessor ip) {

    ExtendedVector2d tempF; // temp vectors for forces
    ExtendedVector2d tempV = new ExtendedVector2d();

    Node n = snake.getHead();
    do { // compute F_total

      if (!n.isFrozen()) {

        // compute F_central
        tempV.setX(n.getNormal().getX() * BOA_.qState.segParam.f_central);
        tempV.setY(n.getNormal().getY() * BOA_.qState.segParam.f_central);
        n.setF_total(tempV);

        // compute F_contract
        tempF = contractionForce(n);
        tempV.setX(tempF.getX() * BOA_.qState.segParam.f_contract);
        tempV.setY(tempF.getY() * BOA_.qState.segParam.f_contract);
        n.addF_total(tempV);

        // compute F_image and F_friction
        tempF = imageForce(n, ip);
        tempV.setX(tempF.getX() * BOA_.qState.segParam.f_image);// - n.getVel().getX() *
        // boap.f_friction);
        tempV.setY(tempF.getY() * BOA_.qState.segParam.f_image);// - n.getVel().getY() *
        // boap.f_friction);
        n.addF_total(tempV);

        // compute new velocities of the node
        tempV.setX(BOA_.qState.boap.delta_t * n.getF_total().getX());
        tempV.setY(BOA_.qState.boap.delta_t * n.getF_total().getY());
        n.addVel(tempV);

        // store the prelimanary point to move the node to
        tempV.setX(BOA_.qState.boap.delta_t * n.getVel().getX());
        tempV.setY(BOA_.qState.boap.delta_t * n.getVel().getY());
        n.setPrelim(tempV); // normal
        // if (BOA_.qState.segParam.contractingDirection == true) {
        // n.setPrelim(tempV); // normal
        // } else {
        // if (n.isFrozen()) {
        // n.setPrelim(new ExtendedVector2d(0, 0)); // expanding
        // } else {
        // n.setPrelim(tempV); // normal
        // }
        // }

        // add some friction
        n.getVel().multiply(BOA_.qState.boap.f_friction);

        // freeze node if vel is below velCrit
        if (n.getVel().length() < BOA_.qState.segParam.vel_crit) {
          snake.freezeNode(n);
        }
      }
      n = n.getNext(); // move to next node
    } while (!n.isHead()); // continue while have not reached tail

    // update all nodes to new positions
    n = snake.getHead();
    do {
      n.update(); // use preliminary variables
      n = n.getNext();
    } while (!n.isHead());

    snake.updateNormals(BOA_.qState.segParam.expandSnake);

    return snake.isFrozen(); // true if all nodes frozen
  }

  /**
   * constrictWrite.
   * 
   * @param snake snake
   * @param ip ip
   * @return true on success
   * @deprecated Strictly related to absolute paths on disk. Probably for testing only.
   */
  public boolean constrictWrite(final Snake snake, final ImageProcessor ip) {
    // for writing forces at each frame
    try {
      PrintWriter pw = new PrintWriter(
              new FileWriter("/Users/rtyson/Documents/phd/tmp/test/forcesWrite/forces.txt"), true);
      ExtendedVector2d tempF; // temp vectors for forces
      ExtendedVector2d tempV = new ExtendedVector2d();

      Node n = snake.getHead();
      do { // compute F_total

        // if (!n.isFrozen()) {

        // compute F_central
        tempV.setX(n.getNormal().getX() * BOA_.qState.segParam.f_central);
        tempV.setY(n.getNormal().getY() * BOA_.qState.segParam.f_central);
        pw.print("\n" + n.getTrackNum() + "," + tempV.length() + ",");
        n.setF_total(tempV);

        // compute F_contract
        tempF = contractionForce(n);
        if (n.getCurvatureLocal() > 0) {
          pw.print(tempF.length() + ",");
        } else {
          pw.print((tempF.length() * -1) + ",");
        }
        tempV.setX(tempF.getX() * BOA_.qState.segParam.f_contract);
        tempV.setY(tempF.getY() * BOA_.qState.segParam.f_contract);
        n.addF_total(tempV);

        // compute F_image and F_friction
        tempF = imageForce(n, ip);
        pw.print((tempF.length() * -1) + ",");
        tempV.setX(tempF.getX() * BOA_.qState.segParam.f_image);// - n.getVel().getX()*
        // boap.f_friction);
        tempV.setY(tempF.getY() * BOA_.qState.segParam.f_image);// - n.getVel().getY()*
        // boap.f_friction);
        n.addF_total(tempV);
        pw.print(n.getF_total().length() + "");

        // compute new velocities of the node
        tempV.setX(BOA_.qState.boap.delta_t * n.getF_total().getX());
        tempV.setY(BOA_.qState.boap.delta_t * n.getF_total().getY());
        n.addVel(tempV);

        // add some friction
        n.getVel().multiply(BOA_.qState.boap.f_friction);

        // store the prelimanary point to move the node to
        tempV.setX(BOA_.qState.boap.delta_t * n.getVel().getX());
        tempV.setY(BOA_.qState.boap.delta_t * n.getVel().getY());
        n.setPrelim(tempV);

        // freeze node if vel is below velCrit
        if (n.getVel().length() < BOA_.qState.segParam.vel_crit) {
          snake.freezeNode(n);
        }
        // }
        n = n.getNext(); // move to next node
      } while (!n.isHead()); // continue while have not reached tail

      // update all nodes to new positions
      n = snake.getHead();
      do {
        n.update(); // use preliminary variables
        n = n.getNext();
      } while (!n.isHead());

      snake.updateNormals(BOA_.qState.segParam.expandSnake);

      pw.close();
      return snake.isFrozen(); // true if all nodes frozen
    } catch (Exception e) {
      return false;
    }
  }

  /**
   * Contraction force for node.
   * 
   * @param n Node
   * @return force vector
   */
  ExtendedVector2d contractionForce(final Node n) {

    ExtendedVector2d resultR;
    ExtendedVector2d resultL;
    ExtendedVector2d force = new ExtendedVector2d();

    // compute the unit vector pointing to the left neighbor (static method)
    resultL = ExtendedVector2d.unitVector(n.getPoint(), n.getPrev().getPoint());

    // compute the unit vector pointing to the right neighbor
    resultR = ExtendedVector2d.unitVector(n.getPoint(), n.getNext().getPoint());

    force.setX((resultR.getX() + resultL.getX()) * 0.5); // combine vector to left and right
    force.setY((resultR.getY() + resultL.getY()) * 0.5);

    return (force);
  }

  /**
   * Calculate image force.
   * 
   * @param n node
   * @param ip image
   * @return image force at node
   */
  ExtendedVector2d imageForce(final Node n, final ImageProcessor ip) {
    ExtendedVector2d result = new ExtendedVector2d();
    ExtendedVector2d tan = n.getTangent(); // Tangent at node
    int i;
    int j; // loop vars

    double a = 0.75; // subsampling factor
    double deltaI; // intensity contrast
    double x;
    double y; // co-ordinates of the norm
    double xt;
    double yt; // co-ordinates of the tangent
    int insideI = 0;
    int outsideI = 0; // Intensity of neighbourhood of a node (insde/outside of the chain)
    int inI = 0;
    int outI = 0; // number of pixels in the neighbourhood of a node

    // determine num pixels and total intensity of neighbourhood: a rectangle with sampleTan x
    // sampleNorm
    for (i = 0; i <= 1. / a * BOA_.qState.segParam.sample_tan; i++) {
      // determine points on the tangent
      xt = n.getPoint().getX() + (a * i - BOA_.qState.segParam.sample_tan / 2) * tan.getX();
      yt = n.getPoint().getY() + (a * i - BOA_.qState.segParam.sample_tan / 2) * tan.getY();

      for (j = 0; j <= 1. / a * BOA_.qState.segParam.sample_norm / 2; ++j) {
        x = xt + a * j * n.getNormal().getX();
        y = yt + a * j * n.getNormal().getY();

        insideI += ip.getPixel((int) x, (int) y);
        inI++;

        x = xt - a * j * n.getNormal().getX();
        y = yt - a * j * n.getNormal().getY();

        outsideI += ip.getPixel((int) x, (int) y);
        outI++;
      }
    }

    deltaI = ((double) insideI / inI - (double) outsideI / outI) / 255.;
    if (deltaI > 0.) { // else remains at zero
      result.setX(-Math.sqrt(deltaI) * n.getNormal().getX());
      result.setY(-Math.sqrt(deltaI) * n.getNormal().getY());
    }

    return (result);
  }

  /**
   * Expand all snakes while preventing overlaps adequately to blowup parameter.
   * 
   * <p>Dead snakes are ignored. Count snakes on frame.
   * 
   * <p>This method blows up live snakes (e.g. from previous frame) and prepare for contracting at
   * current one next frame. snakes are expanded in small steps and at every step their nodes are
   * tested for proximity {@link BOAState.BOAp#proxFreeze}. Only snakes whose centroids are closer
   * than {@link BOAState.BOAp#proximity} are tested for overlapping. Note that this actions happen
   * before contracting snakes around objects. This is preparatory step and if we assume contracting
   * direction ({@link BOAState.SegParam#contractingDirection} is true) so then after this step
   * snakes will not overlap during actual segmentation
   * {@link Constrictor#constrict(Snake, ImageProcessor)}
   * 
   * @param nest nest
   * @param frame frame
   * @throws BoaException on snake scale
   */
  public void loosen(final Nest nest, int frame) throws BoaException {
    int nestSize = nest.size();
    Snake snakeA;
    Snake snakeB;

    double[][] prox = computeProxMatrix(nest);
    // will be negative if blowup is <0
    double stepSize = 0.1 * Math.signum(BOA_.qState.segParam.blowup);
    double steps = (double) BOA_.qState.segParam.blowup / stepSize; // always positive

    for (int i = 0; i < steps; i++) {
      // check for contacts, freeze nodes in contact.
      // Ignore snakes that begin after 'frame'
      for (int si = 0; si < nestSize; si++) {
        // if (nest.getHandler(si).isSnakeHandlerFrozen()) {
        // continue;
        // }
        snakeA = nest.getHandler(si).getLiveSnake();
        if (!snakeA.alive || frame < nest.getHandler(si).getStartFrame()) {
          continue;
        }
        for (int sj = si + 1; sj < nestSize; sj++) {
          snakeB = nest.getHandler(sj).getLiveSnake();
          if (!snakeB.alive || frame < nest.getHandler(si).getStartFrame()) {
            continue;
          }
          // proximity is computed for centroids, this is limit below we test for contact.
          // if snake is big enough it can be not tested even if interact with other
          if (prox[si][sj] > BOA_.qState.boap.proximity) {
            continue; // snakes far away, assume no chance that they will interact
          }
          freezeProx(snakeA, snakeB);
        }

      }

      // scale up all snakes by one step (if node not frozen, or dead) unless they start at this
      // frame or after
      for (int s = 0; s < nestSize; s++) {
        // if (nest.getHandler(s).isSnakeHandlerFrozen()) {
        // continue;
        // }
        snakeA = nest.getHandler(s).getLiveSnake();
        if (snakeA.alive && frame > nest.getHandler(s).getStartFrame()) {
          snakeA.scaleSnake(stepSize, Math.abs(stepSize), true);
        }
      }

    }
  }

  /**
   * Freeze snakes which are close to each other in nest.
   * 
   * @param nest nest to process.
   * @param frame current frame
   * @throws BoaException on error
   */
  public void freezeProxSnakes(final Nest nest, int frame) throws BoaException {
    int nestSize = nest.size();
    Snake snakeA;
    Snake snakeB;

    double[][] prox = computeProxMatrix(nest);

    // check for contacts, freeze nodes in contact.
    // Ignore snakes that begin after 'frame'
    for (int si = 0; si < nestSize; si++) {
      snakeA = nest.getHandler(si).getLiveSnake();
      if (!snakeA.alive || frame < nest.getHandler(si).getStartFrame()) {
        continue;
      }
      for (int sj = si + 1; sj < nestSize; sj++) {
        snakeB = nest.getHandler(sj).getLiveSnake();
        if (!snakeB.alive || frame < nest.getHandler(si).getStartFrame()) {
          continue;
        }
        // proximity is computed for centroids, this is limit below we test for contact.
        // if snake is big enough it can be not tested even if interact with other
        if (prox[si][sj] > BOA_.qState.boap.proximity) {
          continue; // snakes far away, assume no chance that they will interact
        }
        freezeProx(snakeA, snakeB);
      }

    }
  }

  /**
   * Compute distance between centroids of all snakes in nest.
   * 
   * @param nest nest to process
   * @return triangular matrix of centroids distances
   */
  private double[][] computeProxMatrix(final Nest nest) {
    int nestSize = nest.size();
    Snake snakeA;
    Snake snakeB;
    // dist between snake centroids, triangular
    double[][] prox = new double[nestSize][nestSize];
    for (int si = 0; si < nestSize; si++) {
      snakeA = nest.getHandler(si).getLiveSnake();
      snakeA.calcCentroid();
      for (int sj = si + 1; sj < nestSize; sj++) {
        snakeB = nest.getHandler(sj).getLiveSnake();
        snakeB.calcCentroid();
        prox[si][sj] = ExtendedVector2d.lengthP2P(snakeA.getCentroid(), snakeB.getCentroid());
      }
    }
    return prox;
  }

  /**
   * Freeze nodes that are close to each other in two snakes.
   * 
   * <p>This method is called for two snakes whose centroids are closer than
   * {@link BOAState.BOAp#proximity}
   * 
   * @param a snake
   * @param b snake
   * @see #loosen(Nest, int)
   */
  private void freezeProx(final Snake a, final Snake b) {

    Node bn;
    Node an = a.getHead();
    double prox;

    do {
      bn = b.getHead();
      do {
        if (an.isFrozen() && bn.isFrozen()) {
          // an = an.getNext();
          bn = bn.getNext();
          continue;
        }
        // test proximity and freeze
        prox = ExtendedVector2d.distPointToSegment(an.getPoint(), bn.getPoint(),
                bn.getNext().getPoint());
        if (prox < BOA_.qState.boap.proxFreeze) {
          a.freezeNode(an);
          b.freezeNode(bn);
          b.freezeNode(bn.getNext());
          // an.freeze();
          // bn.freeze(); // FIXME using this will exclude Snake.FREEZE from updating. Use from
          // Snake
          // bn.getNext().freeze();
          break;
        }
        bn = bn.getNext();
      } while (!bn.isHead());

      an = an.getNext();
    } while (!an.isHead());

  }

  /**
   * Implode nest.
   * 
   * @param nest nest
   * @param f frame number
   * @throws BoaException on snake.implode
   * @see BOAState.SegParam#expandSnake
   */
  public void implode(final Nest nest, int f) throws BoaException {
    // System.out.println("imploding snake");
    SnakeHandler snakeH;
    Snake snake;
    for (int s = 0; s < nest.size(); s++) {
      snakeH = nest.getHandler(s);
      if (nest.getHandler(s).isSnakeHandlerFrozen()) {
        continue;
      }
      snake = snakeH.getLiveSnake();
      if (snake.alive && f > snakeH.getStartFrame()) {
        snake.implode();
      }
    }
  }
}