package com.srbenoit.microscopy;

import java.awt.AlphaComposite;
import java.awt.Color;
import java.awt.Composite;
import java.awt.Graphics2D;
import java.awt.RenderingHints;
import java.awt.geom.Path2D;
import java.awt.image.BufferedImage;
import java.io.BufferedWriter;
import java.io.File;
import java.io.FileOutputStream;
import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintStream;
import java.util.ArrayList;
import java.util.List;
import java.util.Map;
import java.util.logging.Level;
import javax.imageio.ImageIO;
import com.srbenoit.color.Gradient;
import com.srbenoit.filter.AbstractFilter;
import com.srbenoit.filter.FilterException;
import com.srbenoit.filter.FilterInput;
import com.srbenoit.filter.FilterTreeExecutor;
import com.srbenoit.filter.Pipe;
import com.srbenoit.filter.items.ImageArrayPipeItem;
import com.srbenoit.filter.items.ImagePoint;
import com.srbenoit.filter.items.Trajectory;
import com.srbenoit.filter.items.TrajectoryListPipeItem;
import com.srbenoit.math.delaunay.GraphEdge;
import com.srbenoit.math.delaunay.Vertex;
import com.srbenoit.math.delaunay.Voronoi;
import com.srbenoit.math.grapher.Grapher;
import com.srbenoit.math.grapher.PointListGraph;
import com.srbenoit.math.optimizers.FailedToConvergeException;
import com.srbenoit.math.optimizers.Optimizable;
import com.srbenoit.math.optimizers.Optimizer;
import com.srbenoit.math.optimizers.OutputValue;

/**
 * A filter that performs a diffusion analysis on a set of trajectories.
 */
public class DiffusionAnalysisFilter extends AbstractFilter {

    /** version number for serialization */
    private static final long serialVersionUID = 9105886485480586656L;

    /** PNG file format */
    private static final String PNG = "PNG";

    /** the set of trajectories loaded */
    private final transient List<float[][]> paths;

    /** site data corresponding to each trajectory */
    private final transient List<double[]> siteData;

    /** an offscreen image in which to draw the large graph */
    private final transient Grapher large;

    /** an offscreen image in which to draw the small graph */
    private final transient Grapher small;

    /** a gradient over hue with 100 steps */
    private final Gradient gradient;

    /**
     * Constructs a new <code>DiffusionAnalysisFilter</code>.
     */
    public DiffusionAnalysisFilter() {

        super("Diffusion Analysis", DiffusionAnalysisFilter.class.getName());

        this.inputs.add(new FilterInput(ImageArrayPipeItem.class, "Motion-compensated images"));
        this.inputs.add(new FilterInput(TrajectoryListPipeItem.class, "Trajectories to analyze"));
        makeRenderer();

        this.paths = new ArrayList<float[][]>(50);
        this.siteData = new ArrayList<double[]>(100);
        this.large = new Grapher(800, 600);
        this.small = new Grapher(400, 300);

        this.gradient = new Gradient(360, 1);
    }

    /**
     * Duplicates the filter including all of its settings, but returns an independent object.
     *
     * @return  the duplicated object
     */
    @Override public AbstractFilter duplicate() {

        return new DiffusionAnalysisFilter();
    }

    /**
     * Performs the filter operation.
     *
     * @param   executor  the <code>FilterTreeExecutor</code> that is executing the filter
     * @param   pipe      a pipe containing the input data items
     * @throws  FilterException  if the filter cannot complete
     */
    @Override public void filter(final FilterTreeExecutor executor, final Pipe pipe)
        throws FilterException {

        ImageArrayPipeItem images;
        TrajectoryListPipeItem trajectories;

        validateInputs(pipe);
        executor.indicateProgress(1);

        images = (ImageArrayPipeItem) pipe.get(this.inputs.get(0).getKey());
        trajectories = (TrajectoryListPipeItem) pipe.get(this.inputs.get(1).getKey());
        executor.indicateProgress(2);

        runFilter(executor, pipe, images, trajectories);

        executor.indicateProgress(100);
    }

    /**
     * Runs the filter, reading the source Metamorph TIF files and extracting an array of images,
     * the first dimension of which is time, and the second dimension of which is z plane.
     *
     * @param  executor      the <code>FilterTreeExecutor</code> that is executing the filter
     * @param  dir           the directory where filter output files should be written
     * @param  images        the set of images with maxima marked
     * @param  trajectories  the set of trajectories to process
     */
    private void runFilter(final FilterTreeExecutor executor, final Pipe pipe,
        final ImageArrayPipeItem images, final TrajectoryListPipeItem trajectories) {

        File dir;

        dir = new File(pipe.getDir(), "analysis");

        if (!dir.exists()) {
            dir.mkdir();
        }

        try {
            extractTrajectories(trajectories);
            executor.indicateProgress(30);

            if (!executor.isCancelled()) {
                genReport(dir);
                executor.indicateProgress(50);
                genImages(images, dir);
                executor.indicateProgress(90);
            }
        } catch (IOException e) {
            LOG.log(Level.WARNING, "Failure in analysis", e);
        }
    }

    /**
     * Performs the analysis and generates a report of findings.
     *
     * @param   dir  the directory where the report is to be saved
     * @throws  IOException  if there is an error reading trajectories or writing the report
     */
    private void genReport(final File dir) throws IOException {

        File report;
        File summary;
        BufferedWriter outRep;
        BufferedWriter outSum;
        Optimizer opt;
        double[] scale;
        double[] guess;
        double[] result;
        double[] site;
        double avg;
        float sse;
        float sst;
        float curve;
        float rSquared;
        float deltaA;
        float deltaB;
        float angle;
        int angleBucket;
        int track;
        float dist;
        float speed;

        report = new File(dir, "analysis_report.csv");
        summary = new File(dir, "analysis_summary.csv");

        report.getParentFile().mkdirs();
        makeReadme(dir);

        outRep = new BufferedWriter(new FileWriter(report));
        outSum = new BufferedWriter(new FileWriter(summary));

        try {
            outRep.write("Trajectory analysis" + Pipe.CRLF);
            outSum.write("Summary data" + Pipe.CRLF + Pipe.CRLF);
            outSum.write(
                "track, start frame, length, x(0), y(0), K, alpha, SSE, SST, r^2, angle, angle bucket, distance, avg. speed"
                + Pipe.CRLF);

            // Each trajectory is an array of arrays of floats, each entry
            // containing: frame, x, y, <x>, <y>, x - <x>, y - <y>, d, <d>,
            // d - <d>, (d - <d>)^2, and <(d - <d>)^2>

            // Curve fit each trajectory to traj[i][11] = 4 K t^alpha
            // where K is the diffusion constant, and alpha is the diffusion
            // exponent. To do this, we minimize the following:
            // [ sum_{i=1}^n (4 K i^{2 alpha} - i^alpha traj[i][11]) ]^2 +
            // [ sum_{i=1}^n (4 K i^{2 alpha} (ln i)
            // - traj[i][11] i^alpha (ln i)) ]^2

            scale = new double[] { 0.1, 0.1 };
            guess = new double[] { 1, 1 };

            track = 1;

            for (float[][] traj : this.paths) {
                opt = new Optimizer(new OptimizableFxn(traj), scale, false); // NOPMD SRB

                try {
                    result = opt.optimize(guess, Level.WARNING);

                    avg = 0;

                    for (int i = 0; i < traj.length; i++) {
                        avg += traj[i][10];
                    }

                    avg /= traj.length;

                    sse = 0;
                    sst = 0;
                    outRep.write(Pipe.CRLF + "Track ");
                    outRep.write(Integer.toString(track));
                    outRep.write(Pipe.CRLF
                        + "t, frame, x, y, <x>, <y>, x - <x>, y - <y>, d, <d>, d - <d>, (d - <d>)^2, <(d - <d>)^2>, Curve Fit"
                        + Pipe.CRLF);

                    for (int i = 0; i < traj.length; i++) {
                        curve = (float) (4 * result[0] * Math.pow(i, result[1]));
                        outRep.write("  " + Integer.toString(i) + ", " + Float.toString(traj[i][0])
                            + ", " + Float.toString(traj[i][1]) + ", " + Float.toString(traj[i][2])
                            + ", " + Float.toString(traj[i][3]) + ", " + Float.toString(traj[i][4])
                            + ", " + Float.toString(traj[i][5]) + ", " + Float.toString(traj[i][6])
                            + ", " + Float.toString(traj[i][7]) + ", " + Float.toString(traj[i][8])
                            + ", " + Float.toString(traj[i][9]) + ", "
                            + Float.toString(traj[i][10]) + ", " + Float.toString(traj[i][11])
                            + ", " + Float.toString(curve) + Pipe.CRLF);
                        sse += (curve - traj[i][11]) * (curve - traj[i][11]);
                        sst += (avg - traj[i][11]) * (avg - traj[i][11]);
                    }

                    rSquared = 1 - (sse / sst);

                    deltaA = traj[traj.length - 1][1] - traj[0][1];
                    deltaB = traj[traj.length - 1][2] - traj[0][2];
                    angle = (float) (180 * Math.atan2(deltaB, deltaA) / Math.PI);
                    angleBucket = (int) Math.floor((angle + 12) / 24);
                    dist = (float) Math.sqrt((deltaA * deltaA) + (deltaB * deltaB));
                    speed = dist / traj.length;

                    outRep.write(
                        "start frame, length, x(0), y(0), K, alpha, SSE, SST, r^2, angle, angle bucket, distance, avg. speed"
                        + Pipe.CRLF);
                    outRep.write(Float.toString(traj[0][0]) + ", " + traj.length + ", "
                        + Float.toString(traj[0][1]) + ", " + Float.toString(traj[0][2]) + ", "
                        + result[0] + ", " + result[1] + ", " + sse + ", " + sst + ", " + rSquared
                        + ", " + angle + ", " + angleBucket + ", " + dist + ", " + speed
                        + Pipe.CRLF);
                    site = new double[13]; // NOPMD SRB
                    site[0] = traj[0][0]; // Start frame
                    site[1] = traj.length; // Length
                    site[2] = traj[0][1]; // x(0)
                    site[3] = traj[0][2]; // y(0)
                    site[4] = result[0]; // K
                    site[5] = result[1]; // Alpha
                    site[6] = sse; // SSE
                    site[7] = sst; // SST
                    site[8] = rSquared; // R^2 for curve fit
                    site[9] = angle; // angle of motion
                    site[10] = angleBucket; // angle bucket of motion
                    site[11] = dist; // distance moved
                    site[12] = speed; // average speed
                    this.siteData.add(site);

                    outSum.write(Integer.toString(track) + ", " + Float.toString(traj[0][0]) + ", "
                        + traj.length + ", " + Float.toString(traj[0][1]) + ", "
                        + Float.toString(traj[0][2]) + ", " + result[0] + ", " + result[1] + ", "
                        + sse + ", " + sst + ", " + rSquared + ", " + angle + ", " + angleBucket
                        + ", " + dist + ", " + speed + Pipe.CRLF);

                    drawGraphs(dir, traj, result, track);
                } catch (FailedToConvergeException e) {
                    LOG.warning("Optimizer failed to converge");
                }

                track++;
            }

        } finally {
            outRep.close();
            outSum.close();
        }
    }

    /**
     * Extracts trajectories into arrays of values.
     *
     * @param   trajectories  the trajectories to extract
     * @throws  IOException  if there is an error reading trajectories
     */
    private void extractTrajectories(final TrajectoryListPipeItem trajectories)
        throws IOException {

        Trajectory trajectory;
        ImagePoint imgPoint;
        float[][] traj;
        float totalX;
        float totalY;
        float totalD;
        float total;

        for (int i = 0; i < trajectories.getNumTrajectories(); i++) {

            trajectory = trajectories.getTrajectory(i);

            traj = new float[trajectory.numPoints()][12]; // NOPMD SRB

            for (int j = 0; j < trajectory.numPoints(); j++) {
                imgPoint = trajectory.getPoint(j);

                traj[j][0] = trajectory.getTimePoint(j);
                traj[j][1] = imgPoint.getXPos();
                traj[j][2] = -imgPoint.getYPos(); // Change sign

                totalX = 0;
                totalY = 0;

                for (int k = 0; k <= j; k++) {
                    totalX += traj[k][1];
                    totalY += traj[k][2];
                }

                traj[j][3] = totalX / (j + 1); // [3] is <x> so far
                traj[j][4] = totalY / (j + 1); // [4] is <y> so far

                traj[j][5] = traj[j][1] - traj[j][3]; // [5] is x-<x>
                traj[j][6] = traj[j][2] - traj[j][4]; // [6] is y-<y>
                traj[j][7] = (float) Math.sqrt(((traj[j][1] - traj[0][1])
                            * (traj[j][1] - traj[0][1]))
                        + ((traj[j][2] - traj[0][2]) * (traj[j][2] - traj[0][2]))); // [7] is d

                totalD = 0;

                for (int k = 0; k <= j; k++) {
                    totalD += traj[k][7];
                }

                traj[j][8] = totalD / (j + 1); // [8] is <d> so far
                traj[j][9] = traj[j][7] - traj[j][8]; // [9] is d-<d>
                traj[j][10] = traj[j][9] * traj[j][9]; // [10] is (d-<d>)^2

                total = 0;

                for (int k = 0; k <= j; k++) {
                    total += traj[k][10];
                }

                traj[j][11] = total / (j + 1); // [11] is <(d - <d>)^2>
            }

            this.paths.add(traj);
        }
    }

    /**
     * Generates images based on colorings of Voronoi cells that correspond to cell movement
     * parameters.
     *
     * @param   images  the set of images with maxima marked
     * @param   dir     the directory where images are to be saved
     * @throws  IOException  if there is an error reading trajectories or writing the report
     */
    private void genImages(final ImageArrayPipeItem images, final File dir) throws IOException {

        BufferedImage base;

        // Load the base image that we will use for overlays
        base = images.getImage(0, 0);

        if (base != null) {
            makeImages(base, dir);
        }
    }

    /**
     * Builds a series of overlay images.
     *
     * @param   base  the base image on which to overlay
     * @param   dir   the directory where images are to be saved
     * @throws  IOException  if there is an error writing an image file
     */
    private void makeImages(final BufferedImage base, final File dir) throws IOException {

        Vertex[] vertices;
        double[] site;
        Voronoi vor;
        List<GraphEdge> edges;
        BufferedImage overlay;
        Graphics2D grx;
        int index;
        Color color;
        Composite orig;
        double min;
        double max;
        Map<Vertex, Path2D> map;
        double sat;

        // Build the vertex list based on the initial points of all
        // trajectories, with Y coordinate negated
        vertices = new Vertex[this.siteData.size()];

        for (int i = 0; i < this.siteData.size(); i++) {
            site = this.siteData.get(i);
            vertices[i] = new Vertex(site[2], -site[3]); // NOPMD SRB
        }

        vor = new Voronoi(1);
        edges = vor.generateVoronoi(vertices);
        map = vor.makeVoronoiPolygons(vertices, edges);

        overlay = new BufferedImage(base.getWidth(), base.getHeight(), BufferedImage.TYPE_INT_RGB);
        grx = (Graphics2D) overlay.getGraphics();
        grx.setRenderingHint(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON);

        // Image 1: Just the Voronoi cell boundaries
        grx.drawImage(base, 0, 0, null);
        grx.setColor(Color.YELLOW);

        for (GraphEdge edge : edges) {
            grx.drawLine((int) edge.xPos1, (int) edge.yPos1, (int) edge.xPos2, (int) edge.yPos2);
        }

        grx.setColor(Color.RED);

        for (Vertex vert : vertices) {
            grx.fillOval((int) vert.xPos - 1, (int) vert.yPos - 1, 3, 3);
        }

        ImageIO.write(overlay, PNG, new File(dir, "Voronoi_boundaries.png"));

        //
        //
        //

        // Image 2: Color hue based on angle, saturation based on K
        min = this.siteData.get(0)[4];
        max = min;

        for (int i = 1; i < this.siteData.size(); i++) {
            site = this.siteData.get(i);

            if (site[4] < min) {
                min = site[4];
            }

            if (site[4] > max) {
                max = site[4];
            }
        }

        grx.drawImage(base, 0, 0, null);
        orig = grx.getComposite();
        grx.setComposite(AlphaComposite.getInstance(AlphaComposite.SRC_OVER, 0.5f));

        for (Vertex vert : map.keySet()) {

            for (int i = 0; i < this.siteData.size(); i++) {
                site = this.siteData.get(i);

                if ((vert.xPos == site[2]) && (vert.yPos == -site[3])) {

                    // Get hue based on angle
                    index = (int) site[9];

                    while (index < 0) {
                        index += 360;
                    }

                    while (index >= 360) {
                        index -= 360;
                    }

                    color = this.gradient.getColor(index);

                    // Adjust saturation based on alpha
                    sat = (site[4] - min) / (max - min);

                    color = new Color((int) (color.getRed() * sat), // NOPMD SRB
                            (int) (color.getGreen() * sat), (int) (color.getBlue() * sat));

                    grx.setColor(color);
                    grx.fill(map.get(vert));

                    break;
                }
            }
        }

        grx.setComposite(orig);
        grx.setColor(Color.RED);

        for (Vertex vert : vertices) {
            grx.fillOval((int) vert.xPos - 1, (int) vert.yPos - 1, 3, 3);
        }

        drawColorWheel(grx);
        ImageIO.write(overlay, PNG, new File(dir, "Voronoi_k.png"));

        //
        //
        //

        // Image 3: Color hue based on angle, saturation based on ALPHA
        min = this.siteData.get(0)[5];
        max = min;

        for (int i = 1; i < this.siteData.size(); i++) {
            site = this.siteData.get(i);

            if (site[5] < min) {
                min = site[5];
            }

            if (site[5] > max) {
                max = site[5];
            }
        }

        grx.drawImage(base, 0, 0, null);
        orig = grx.getComposite();
        grx.setComposite(AlphaComposite.getInstance(AlphaComposite.SRC_OVER, 0.5f));

        for (Vertex vert : map.keySet()) {

            for (int i = 0; i < this.siteData.size(); i++) {
                site = this.siteData.get(i);

                if ((vert.xPos == site[2]) && (vert.yPos == -site[3])) {

                    // Get hue based on angle
                    index = (int) site[9];

                    while (index < 0) {
                        index += 360;
                    }

                    while (index >= 360) {
                        index -= 360;
                    }

                    color = this.gradient.getColor(index);

                    // Adjust saturation based on alpha
                    sat = (site[5] - min) / (max - min);

                    color = new Color((int) (color.getRed() * sat), // NOPMD SRB
                            (int) (color.getGreen() * sat), (int) (color.getBlue() * sat));

                    grx.setColor(color);
                    grx.fill(map.get(vert));

                    break;
                }
            }
        }

        grx.setComposite(orig);
        grx.setColor(Color.RED);

        for (Vertex vert : vertices) {
            grx.fillOval((int) vert.xPos - 1, (int) vert.yPos - 1, 3, 3);
        }

        drawColorWheel(grx);
        ImageIO.write(overlay, PNG, new File(dir, "Voronoi_alpha.png"));

        //
        //
        //

        // Image 4: Color hue based on angle, saturation based on dist
        min = this.siteData.get(0)[11];
        max = min;

        for (int i = 1; i < this.siteData.size(); i++) {
            site = this.siteData.get(i);

            if (site[11] < min) {
                min = site[11];
            }

            if (site[11] > max) {
                max = site[11];
            }
        }

        grx.drawImage(base, 0, 0, null);
        orig = grx.getComposite();
        grx.setComposite(AlphaComposite.getInstance(AlphaComposite.SRC_OVER, 0.5f));

        for (Vertex vert : map.keySet()) {

            for (int i = 0; i < this.siteData.size(); i++) {
                site = this.siteData.get(i);

                if ((vert.xPos == site[2]) && (vert.yPos == -site[3])) {

                    // Get hue based on angle
                    index = (int) site[9];

                    while (index < 0) {
                        index += 360;
                    }

                    while (index >= 360) {
                        index -= 360;
                    }

                    color = this.gradient.getColor(index);

                    // Adjust saturation based on alpha
                    sat = (site[11] - min) / (max - min);

                    color = new Color((int) (color.getRed() * sat), // NOPMD SRB
                            (int) (color.getGreen() * sat), (int) (color.getBlue() * sat));

                    grx.setColor(color);
                    grx.fill(map.get(vert));

                    break;
                }
            }
        }

        grx.setComposite(orig);
        grx.setColor(Color.RED);

        for (Vertex vert : vertices) {
            grx.fillOval((int) vert.xPos - 1, (int) vert.yPos - 1, 3, 3);
        }

        drawColorWheel(grx);
        ImageIO.write(overlay, PNG, new File(dir, "Voronoi_dist.png"));

        //
        //
        //

        // Image 5: Color hue based on angle, saturation based on speed
        min = this.siteData.get(0)[12];
        max = min;

        for (int i = 1; i < this.siteData.size(); i++) {
            site = this.siteData.get(i);

            if (site[12] < min) {
                min = site[12];
            }

            if (site[12] > max) {
                max = site[12];
            }
        }

        grx.drawImage(base, 0, 0, null);
        orig = grx.getComposite();
        grx.setComposite(AlphaComposite.getInstance(AlphaComposite.SRC_OVER, 0.5f));

        for (Vertex vert : map.keySet()) {

            for (int i = 0; i < this.siteData.size(); i++) {
                site = this.siteData.get(i);

                if ((vert.xPos == site[2]) && (vert.yPos == -site[3])) {

                    // Get hue based on angle
                    index = (int) site[9];

                    while (index < 0) {
                        index += 360;
                    }

                    while (index >= 360) {
                        index -= 360;
                    }

                    color = this.gradient.getColor(index);

                    // Adjust saturation based on alpha
                    sat = (site[12] - min) / (max - min);

                    color = new Color((int) (color.getRed() * sat), // NOPMD SRB
                            (int) (color.getGreen() * sat), (int) (color.getBlue() * sat));

                    grx.setColor(color);
                    grx.fill(map.get(vert));

                    break;
                }
            }
        }

        grx.setComposite(orig);
        grx.setColor(Color.RED);

        for (Vertex vert : vertices) {
            grx.fillOval((int) vert.xPos - 1, (int) vert.yPos - 1, 3, 3);
        }

        drawColorWheel(grx);
        ImageIO.write(overlay, PNG, new File(dir, "Voronoi_speed.png"));

        // site[0] = traj[0][0]; // Start frame
        // site[1] = traj.length; // Length
        // site[2] = traj[0][1]; // x(0)
        // site[3] = traj[0][2]; // y(0)
        // site[4] = result[0]; // K
        // site[5] = result[1]; // Alpha
        // site[6] = sse; // SSE
        // site[7] = sst; // SST
        // site[8] = rSquared; // R^2 for curve fit
        // site[9] = angle; // angle of motion
        // site[10] = angleBucket; // angle bucket of motion
        // site[11] = dist; // distance moved
        // site[12] = speed; // average speed

    }

    /**
     * Draws the color wheel.
     *
     * @param  grx  the <code>Graphics</code> to which to draw
     */
    private void drawColorWheel(final Graphics2D grx) {

        double angle;
        double cos;
        double sin;
        Color color;

        for (int index = 0; index < 360; index++) {
            color = this.gradient.getColor(index);

            angle = index * Math.PI / 180;
            cos = Math.cos(angle);
            sin = Math.sin(angle);

            grx.setColor(color);
            grx.drawLine((int) (30 - (15 * cos)), (int) (30 - (15 * sin)), (int) (30 - (20 * cos)),
                (int) (30 - (20 * sin)));
        }
    }

    /**
     * Computes twice the area of the oriented triangle <code>(xPos1, yPos1), (xPos2, y3), (x3,
     * y3)</code> (the area is positive if the triangle is oriented counterclockwise).
     *
     * @param   xPos1  the first point X coordinate
     * @param   yPos1  the first point Y coordinate
     * @param   xPos2  the second point X coordinate
     * @param   yPos2  the second point Y coordinate
     * @param   xPos3  the third point X coordinate
     * @param   yPos3  the third point Y coordinate
     * @return  twice the oriented area
     */
    public double triArea(final double xPos1, final double yPos1, final double xPos2,
        final double yPos2, final double xPos3, final double yPos3) {

        return ((xPos2 - xPos1) * (yPos3 - yPos1)) - ((yPos2 - yPos1) * (xPos3 - xPos1));
    }

    /**
     * Tests whether the points of a triangle are in counterclockwise order.
     *
     * @param   xPos1  the first point X coordinate
     * @param   yPos1  the first point Y coordinate
     * @param   xPos2  the second point X coordinate
     * @param   yPos2  the second point Y coordinate
     * @param   xPos3  the third point X coordinate
     * @param   yPos3  the third point Y coordinate
     * @return  <code>true</code> if triangle <code>(xPos1, yPos1), (xPos2, y3), (x3, y3)</code> is
     *          in counterclockwise order
     */
    public boolean isCcw(final double xPos1, final double yPos1, final double xPos2,
        final double yPos2, final double xPos3, final double yPos3) {

        return triArea(xPos1, yPos1, xPos2, yPos2, xPos3, yPos3) > 0;
    }

    /**
     * Generates the "read me" file in the analysis directory.
     *
     * @param   dir  the analysis directory
     * @throws  IOException  if there is an error writing the file
     */
    private void makeReadme(final File dir) throws IOException {

        File readme;
        PrintStream out;

        readme = new File(dir, "analysis_readme.txt");

        out = new PrintStream(new FileOutputStream(readme));

        out.println("Summary of analysis files generated:");
        out.println("");
        out.println("analysis_report.csv:");
        out.println("    This comma-separated-value (CSV) file, which can be opened directly in");
        out.println("Excel or similar programs contains a complete analysis of the trajectories");
        out.println("extracted from the frames. The file is divided into sections, each labeled");
        out.println("'Track n' where n ranges from 1 to the number of trajectories. Within each");
        out.println("section are the following columns:");
        out.println("  t: the time since the trajectory's start");
        out.println("  frame: the frame number");
        out.println("  x: the X coordinate of the cell (corrected for tissue motion for t > 0)");
        out.println("  y: the Y coordinate of the cell (corrected for tissue motion for t > 0)");
        out.println("  a: the X coordinate, rotated through an angle");
        out.println("  b: the Y coordinate, rotated through an angle");
        out.println("  <a>: the mean of the 'a' coordinates so far");
        out.println("  <b>: the mean of the 'b' coordinates so far");
        out.println("  a - <a>: difference between <a> and a");
        out.println("  b - <b>: difference between <b> and b");
        out.println("  d - <d>: difference between distance from start (d) and mean distance");
        out.println("           from start (<d>) so far");
        out.println("  (d - <d>)^2: square of (d - <d>)");
        out.println("  <(d - <d>)^2>: mean squared differences of distance of the cell from");
        out.println("                 its starting point and mean distance from start");
        out.println("  Curve Fit: corresponding y value on best-fit curve of the form");
        out.println("             y = 2 d K t^alpha, where d is dimension (2)");
        out.println("             [Metzler R, Klafter J, \"The restaurant at the end of the");
        out.println("              random walk: Recent developments in the description of");
        out.println("              anomalous transport by fractional dynamics.\", 2004.");
        out.println("              J Phys A, 37(31):R161�R208]");
        out.println("At the bottom of each section is summary data:");
        out.println("  Start frame: the frame number where the trajectory began");
        out.println("  x(0): the X coordinate where the cell began");
        out.println("  y(0): the Y coordinate where the cell began");
        out.println("  K: the diffusion constant in the best-fit curve");
        out.println("  alpha: the diffusion exponent in the best-fit curve");
        out.println("  SSE: the sum of squared errors for the curve fit");
        out.println("  SST: the total variation of the curve fit");
        out.println("  r^2: the coefficient of determination of the curve fit");
        out.println("  angle: the angle of the trajectory (based on end points)");
        out.println("  angle bucket: assignment of angle into a group of similar angles (*)");
        out.println("");
        out.println("analysis_summary.csv:");
        out.println("    This comma-separated-value (CSV) file, aggregates just the summary");
        out.println("lines that appear in the analysis_report.csv  file.");
        out.println("and is provided as a convenience for generating spreadsheet charts.");
        out.println("");
        out.println("report.txt:");
        out.println("    As in all steps of the analysis process, this file records the amount");
        out.println("of time the processing step took, and its presence indicates that a given");
        out.println("phase of the analysis has been completed.  Deleting this file will cause");
        out.println("the alalysis program to re-run this phase of the analysis.");
        out.println("");
        out.println("Track_n_curvefit_lg.png");
        out.println(
            "    This series of large plots graph the <(d - <d>)^2> value (shown in blue)");
        out.println(
            "and the corresponding best-fit curve (shown in dark red) for each trajectory.");
        out.println("Labeling has intentionally been kept to a minimim to make it easier to add");
        out.println("custom labels to the graphs for inclusion in publications.");
        out.println("");
        out.println("Track_n_curvefit_sm.png");
        out.println("    These are simply smaller versions of the Track_n_curvefit_lg.png plots,");
        out.println(
            "but with the same font size, making them clearer when printed in small size.");
        out.println("");
        out.println("");
        out.println("");
        out.println("* Angles are measured from the positive X axis, with angle increasing");
        out.println("  counterclockwise. Angle bucket allows grouping of trajectories with");
        out.println("  others of similar direction.  Each bucket subtends 24 degrees of arc.");
        out.println("  Bucket 0 represents angles from -12 to +12 degrees.  Bucket +1 goes");
        out.println("  from +12 to +36, and bucket -1 goes from -12 to -36.  Therefore, the");
        out.println("  possible range is from -8 to +8.");
        out.println("");
        out.println("  A useful analysis is to group trajectories by some measure (distance,");
        out.println("  angle, speed, or start position), then generate a histogram of how");
        out.println("  many cells fell into each bucket, correlating this count with the");
        out.println("  selected measure, as a way to correlate directional bias with that");
        out.println("  measure.");

        out.close();
    }

    /**
     * Renders the graphs and writes them out to PNG files.
     *
     * @param  dir       the directory in which to export graph files
     * @param  traj      the list of trajectory points
     * @param  result    the result of the curve fit (K, alpha)
     * @param  trackNum  the track number (matches the CSV file)
     */
    private void drawGraphs(final File dir, final float[][] traj, final double[] result,
        final int trackNum) {

        double[] xCoords;
        double[] actual;
        double[] fit;
        PointListGraph actCurve;
        PointListGraph fitCurve;
        File file;

        xCoords = new double[traj.length];
        actual = new double[traj.length];
        fit = new double[traj.length];

        for (int i = 0; i < traj.length; i++) {
            xCoords[i] = traj[i][0] - traj[0][0];
            actual[i] = traj[i][11];
            fit[i] = (float) (4 * result[0] * Math.pow(i, result[1]));
        }

        actCurve = new PointListGraph(xCoords, actual);
        fitCurve = new PointListGraph(xCoords, fit);

        file = new File(dir, "Track_" + trackNum + "_curvefit_lg.png");
        this.large.graph(actCurve, fitCurve);

        try {
            this.large.export(file);
        } catch (IOException e) {
            LOG.log(Level.WARNING, "Failed to export large graph", e);
        }

        file = new File(dir, "Track_" + trackNum + "_curvefit_sm.png");
        this.small.graph(actCurve, fitCurve);

        try {
            this.small.export(file);
        } catch (IOException e) {
            LOG.log(Level.WARNING, "Failed to export small graph", e);
        }
    }

    /**
     * A function that supports optimization.
     */
    private static class OptimizableFxn implements Optimizable {

        /** the data to which to fit the curve */
        private final transient float[][] data;

        /**
         * Constructs a new <code>OptimizableFxn</code>.
         *
         * @param  theData  the data to which to fit the curve
         */
        private OptimizableFxn(final float[][] theData) {

            this.data = theData.clone();
        }

        /**
         * Gets the dimension of the function.
         *
         * @return  the dimension
         */
        public int dimension() {

            return 2;
        }

        /**
         * Evaluates the function for a given set of parameters that are to be optimized.
         *
         * @param   parameters  the parameters ([0] is K, [1] is alpha)
         * @return  the function value based on the parameters
         */
        public OutputValue evaluate(final double[] parameters) {

            OutputValue output;
            double term;

            output = new OutputValue();

            for (int i = 1; i < this.data.length; i++) {
                term = (4 * parameters[0] * Math.pow(i, parameters[1])) - this.data[i][11];
                output.setValue(output.getValue() + (term * term));
            }

            if (parameters[1] < 0) {
                output.setOutOfRange(true);
                output.setValue(Double.POSITIVE_INFINITY);
            }

            return output;
        }

        /**
         * Called when the optimizer has found a new set of parameters that result in a lower value
         * during its search for the optimal parameters.
         *
         * @param  parameters  the newly found parameters
         */
        public void updatedParameters(final double[] parameters) {

            // No action
        }
    }

    /**
     * Generates the string representation of the filter.
     *
     * @return  the string representation
     */
    @Override public String toString() {

        return "TrajectoryFilter";
    }
}