View Javadoc
1   package com.github.celldynamics.quimp.utils.graphics;
2   
3   import java.awt.Rectangle;
4   import java.io.BufferedOutputStream;
5   import java.io.FileOutputStream;
6   import java.io.IOException;
7   import java.io.PrintWriter;
8   import java.util.Arrays;
9   import java.util.List;
10  
11  import org.apache.commons.lang3.ArrayUtils;
12  import org.scijava.vecmath.Point2d;
13  import org.scijava.vecmath.Vector2d;
14  import org.slf4j.Logger;
15  import org.slf4j.LoggerFactory;
16  
17  import com.github.celldynamics.quimp.geom.ExtendedVector2d;
18  import com.github.celldynamics.quimp.plugin.qanalysis.STmap;
19  import com.github.celldynamics.quimp.plugin.utils.IPadArray;
20  import com.github.celldynamics.quimp.utils.QuimPArrayUtils;
21  import com.github.celldynamics.quimp.utils.graphics.svg.SVGwritter;
22  
23  /**
24   * Generate polar plots of motility speed along cell perimeter.
25   * 
26   * <p>Use basic mapping - location of point on polar plot depends on its position on cell outline.
27   * 
28   * @author p.baniukiewicz
29   *
30   */
31  public class PolarPlot {
32  
33    /**
34     * The Constant LOGGER.
35     */
36    static final Logger LOGGER = LoggerFactory.getLogger(PolarPlot.class.getName());
37    private STmap mapCell;
38    private Point2d gradientcoord;
39    /**
40     * Size and position of left upper corner of plot area. Most code below assumes that area is
41     * square and centered in 0,0.
42     */
43    protected Rectangle plotArea;
44    /**
45     * Distance of polar plot from 0. Rescaling factor.
46     */
47    public double kscale;
48    /**
49     * Distance of plot from edge. Rescaling factor.
50     */
51    public double uscale;
52  
53    /**
54     * Whether to plot rho and theta labels.
55     */
56    public boolean labels = false;
57  
58    /**
59     * Create PolarPlot object with default plot size.
60     * 
61     * @param mapCell mapCell
62     * @param gradientcoord coordinates of gradient point. Gradient is a feeding of cell put into
63     *        cell environment.
64     */
65    public PolarPlot(STmap mapCell, Point2d gradientcoord) {
66      this.mapCell = mapCell;
67      this.gradientcoord = gradientcoord;
68      plotArea = new Rectangle(-3, -3, 6, 6);
69      kscale = 0.1;
70      uscale = 0.05;
71    }
72  
73    /**
74     * Compute shift for every frame. Shift value indicates which index of outline point should be
75     * first in maps. This point is closest to <tt>gradientcoord</tt>.
76     * 
77     * @return Indexes of first points (x-coordinate) for map for every frame (y-cordinate)
78     */
79    int[] getShift() {
80      int[] ret = new int[mapCell.getT()]; // shift for every frame
81      for (int f = 0; f < mapCell.getT(); f++) { // along frames
82        double dist = Double.MAX_VALUE; // closest point for current frame
83        for (int i = 0; i < mapCell.getRes(); i++) { // along points
84          Point2d p = new Point2d(mapCell.getxMap()[f][i], mapCell.getyMap()[f][i]); // outline point
85          double disttmp = p.distance(gradientcoord); // distance from gradinet point
86          if (disttmp < dist) { // we have closer point
87            dist = disttmp;
88            ret[f] = i; // remember index of closer point
89          }
90        }
91      }
92      return ret;
93    }
94  
95    /**
96     * Compute mass centres for every frame.
97     * 
98     * @return Vector of mass centers for every frame.
99     */
100   Point2d[] getMassCentre() {
101     Point2d[] ret = new Point2d[mapCell.getT()];
102     double[] xmeans = QuimPArrayUtils.getMeanR(mapCell.getxMap());
103     double[] ymeans = QuimPArrayUtils.getMeanR(mapCell.getyMap());
104     for (int f = 0; f < mapCell.getT(); f++) {
105       ret[f] = new Point2d(xmeans[f], ymeans[f]);
106     }
107     return ret;
108   }
109 
110   /**
111    * Compute vectors for one frame between mass centre and outline point. Vectors are in order
112    * starting from closest point. This is representation of outline as vectors.
113    * 
114    * @param f
115    * @param mass
116    * @param shift
117    * 
118    * @return List of vectors starting from closes to gradientcoord.
119    */
120   Vector2d[] getVectors(int f, Point2d[] mass, int[] shift) {
121     Vector2d[] ret = new Vector2d[mapCell.getRes()];
122     int start = shift[f];
123     Point2d mc = mass[f];
124     int l = 0; // output index
125     for (int i = start; i < mapCell.getRes() + start; i++) { // first point is that shifted
126       // true array index
127       int index = IPadArray.getIndex(mapCell.getRes(), i, IPadArray.CIRCULARPAD);
128       // outline point
129       Point2d p = new Point2d(mapCell.getxMap()[f][index], mapCell.getyMap()[f][index]);
130       p.sub(mc); // p = p-mc - vector from centre to point
131       ret[l++] = new Vector2d(p); // put [index] as first
132     }
133     return ret;
134   }
135 
136   /**
137    * Get values from selected map shifting it according to shift.
138    * 
139    * @param f Frame to get.
140    * @param shift Shift value
141    * @param map
142    * @return Vector of map values with first value closest to gradientcoord
143    */
144   double[] getRadius(int f, int shift, double[][] map) {
145     double[] ret = new double[map[f].length];
146     int l = 0; // output index
147     for (int i = shift; i < mapCell.getRes() + shift; i++) {
148       int index = IPadArray.getIndex(map[f].length, i, IPadArray.CIRCULARPAD);
149       ret[l++] = map[f][index];
150     }
151     return ret;
152   }
153 
154   /**
155    * Compute angles between reference vector and vectors.
156    * 
157    * @param vectors array of vectors (in correct order)
158    * @param ref reference vector
159    * @return angles between reference vector and all <tt>vectors</tt>
160    */
161   double[] getAngles(Vector2d[] vectors, Vector2d ref) {
162     double[] ret = new double[vectors.length];
163     for (int i = 0; i < vectors.length; i++) {
164       double a1 = Math.atan2(vectors[i].y, vectors[i].x);
165       double a2 = Math.atan2(ref.y, ref.x);
166       ret[i] = -a1 + a2;
167       // convert to 4-squares angle (left comment to comp with matlab plotPolarPlot)
168       // ret[i] = (ret[i] < 0) ? (ret[i] + 2 * Math.PI) : ret[i];
169     }
170     return ret;
171 
172   }
173 
174   /**
175    * Polar plot of one frame.
176    * 
177    * @param filename filename
178    * @param frame frame
179    * @throws IOException on file error
180    * @see #generatePlot
181    */
182   public void generatePlotFrame(String filename, int frame) throws IOException {
183     int[] shifts = getShift(); // calculate shifts of points according to gradientcoord
184     // shift motility
185     double[] magn = getRadius(frame, shifts[frame], mapCell.getMotMap());
186     // generate vector of arguments
187     double[] angles = getUniformAngles(magn.length);
188     // remove negative values (shift)
189     double min = QuimPArrayUtils.arrayMin(magn);
190     double max = QuimPArrayUtils.arrayMax(magn);
191     for (int i = 0; i < magn.length; i++) {
192       magn[i] -= (min + kscale * min); // k*min - distance from 0
193     }
194 
195     if (labels) {
196       polarPlotPoints(filename, angles, magn, new double[] { min, max });
197     } else {
198       polarPlotPoints(filename, angles, magn, null);
199     } // generate plot
200 
201   }
202 
203   /**
204    * Polar plot of mean of motility along frames.
205    * 
206    * <p>The position of point on polar plot depends on its position on cell outline, First point
207    * after shifting is that closest to selected gradient. It is plotted
208    * on angle 0.
209    * 
210    * @param filename Name of the svg file.
211    * @throws IOException on file save
212    */
213   public void generatePlot(String filename) throws IOException {
214     int[] shifts = getShift(); // calculate shifts of points according to gradientcoord
215     // contains magnitudes of polar plot (e.g. motility) shifted, so first point is that
216     // related to gradientcoord, for every frame [frames][outline points]
217     double[][] magnF = new double[mapCell.getT()][];
218     for (int f = 0; f < mapCell.getT(); f++) {
219       // shift motility for every frame to have gradientcoord related point on 0 index
220       magnF[f] = getRadius(f, shifts[f], mapCell.getMotMap());
221     }
222     // generate vector of arguments
223     double[] angles = getUniformAngles(magnF[0].length);
224 
225     double[] magn = QuimPArrayUtils.getMeanC(magnF); // compute mean for map on frames
226     // rescale to remove negative values and make distance from 0 point
227     double min = QuimPArrayUtils.arrayMin(magn);
228     double max = QuimPArrayUtils.arrayMax(magn);
229     for (int i = 0; i < magn.length; i++) {
230       magn[i] -= (min + kscale * min); // k*min - distance from 0
231     }
232     if (labels) {
233       polarPlotPoints(filename, angles, magn, new double[] { min, max });
234     } else {
235       polarPlotPoints(filename, angles, magn, null);
236     } // generate plot
237     LOGGER.debug("Polar plot saved: " + filename);
238 
239   }
240 
241   /**
242    * Generate svg plot of points.
243    * 
244    * @param filename name of svg file
245    * @param angles vector of arguments (angles) generated {@link #getUniformAngles(int)}
246    * @param magn vector of values related to <tt>angles</tt>
247    * @param valMinMax array of [minVal maxVal] of magn before scaling. Note that this method expect
248    *        magn already scaled to positive values. If null or empty axis legend will not show
249    * @throws IOException on file save
250    */
251   private void polarPlotPoints(String filename, double[] angles, double[] magn, double[] valMinMax)
252           throws IOException {
253     // scale for plotting (6/2) - half of plot area size as plotted from centre)
254     double plotscale = plotArea.getWidth() / 2 / QuimPArrayUtils.arrayMax(magn);
255     plotscale -= plotscale * uscale; // move a little from edge
256     // generate svg
257     BufferedOutputStream out;
258     out = new BufferedOutputStream(new FileOutputStream(filename)); // TODO change to FileWriter
259     PrintWriter osw = new PrintWriter(out, true);
260     Rectangle extendedPlotArea = (Rectangle) plotArea.clone(); // extend square for axis labels
261     extendedPlotArea.grow(1, 1);
262     SVGwritter.writeHeader(osw, extendedPlotArea); // write header with sizes
263     // plot axes
264     SVGwritter.QPolarAxes qaxes = new SVGwritter.QPolarAxes(plotArea);
265     qaxes.angleLabel = labels;
266     if (valMinMax != null && valMinMax.length != 0) {
267       List<Double> ret =
268               ExtendedVector2d.linspace(valMinMax[0], valMinMax[1], qaxes.numofIntCircles + 2);
269       // +2 because of middle and outer (max)
270       // remove first and alst (just to not display it on inner and outer circle)
271       ret.remove(0);
272       ret.remove(ret.size() - 1);
273       // qaxes.radiusLabels = new double[] { 1e-5, 2, 3, 4, 5 };
274       qaxes.radiusLabels = ArrayUtils.toPrimitive(ret.toArray(new Double[0]));
275     }
276     qaxes.draw(osw);
277     // plot points
278     for (int i = 0; i < angles.length; i++) {
279       // x1;y1, x2;y2 two points that define line segment
280       double x = Math.cos(angles[i]) * magn[i];
281       double y = Math.sin(angles[i]) * magn[i];
282       double x1 = Math.cos(angles[(i + 1) % angles.length]) * magn[(i + 1) % angles.length];
283       double y1 = Math.sin(angles[(i + 1) % angles.length]) * magn[(i + 1) % angles.length];
284       // LOGGER.trace("Point coords:" + x + " " + y + " Polar coords:"
285       // + Math.toDegrees(angles[i]) + " " + magn[i]);
286       SVGwritter.Qline ql =
287               new SVGwritter.Qline(x * plotscale, y * plotscale, x1 * plotscale, y1 * plotscale);
288       ql.thickness = 0.01;
289       ql.draw(osw);
290     }
291 
292     osw.write("</svg>\n");
293     osw.close();
294   }
295 
296   /**
297    * Generate uniformly distributed angles for given resolution.
298    * 
299    * <p>Generate angles for polar plot (related to plot, not for position of outline points) assume
300    * basic mapping - first outline point at angle 0, second at angle delta, etc CCW system is
301    * used, but angles are counted in CW. IV and III quarter are negative, then II and I are
302    * positive.
303    * 
304    * @param res Number of angles to generate.
305    * @return Array of angles in radians counted CW, IV quarter is first and negative, II quarter
306    *         is positive, e.g. -1,-2,...-90,...-180,179,178...,90,....0
307    */
308   private double[] getUniformAngles(int res) {
309     double[] angles = new double[res];
310     double delta = (2 * Math.PI - 0) / (res - 1);
311     for (int i = 0; i < angles.length; i++) {
312       angles[i] = -(0 + delta * i);
313       if (angles[i] < -Math.PI) {
314         angles[i] = Math.PI + Math.PI + angles[i]; // negative in II q is changed to positv
315       }
316     }
317     return angles;
318   }
319 
320   /**
321    * http://www.java2s.com/Code/Java/Collections-Data-Structure/LinearInterpolation.htm
322    * 
323    * @param x coordinates
324    * @param y coordinates
325    * @param xi value to interpolate yi
326    * @return linear interpolation of y=f(xi)
327    * @throws IllegalArgumentException wrong input arrays size
328    */
329   public static double[] interpLinear(double[] x, double[] y, double[] xi)
330           throws IllegalArgumentException {
331 
332     if (x.length != y.length) {
333       throw new IllegalArgumentException("X and Y must be the same length");
334     }
335     if (x.length == 1) {
336       throw new IllegalArgumentException("X must contain more than one value");
337     }
338     double[] dx = new double[x.length - 1];
339     double[] dy = new double[x.length - 1];
340     double[] slope = new double[x.length - 1];
341     double[] intercept = new double[x.length - 1];
342 
343     // Calculate the line equation (i.e. slope and intercept) between each point
344     for (int i = 0; i < x.length - 1; i++) {
345       dx[i] = x[i + 1] - x[i];
346       if (dx[i] == 0) {
347         throw new IllegalArgumentException(
348                 "X must be montotonic. A duplicate " + "x-value was found");
349       }
350       if (dx[i] < 0) {
351         throw new IllegalArgumentException("X must be sorted");
352       }
353       dy[i] = y[i + 1] - y[i];
354       slope[i] = dy[i] / dx[i];
355       intercept[i] = y[i] - x[i] * slope[i];
356     }
357     // Perform the interpolation here
358     double[] yi = new double[xi.length];
359     for (int i = 0; i < xi.length; i++) {
360       if ((xi[i] > x[x.length - 1]) || (xi[i] < x[0])) {
361         yi[i] = Double.NaN;
362       } else {
363         int loc = Arrays.binarySearch(x, xi[i]);
364         if (loc < -1) {
365           loc = -loc - 2;
366           yi[i] = slope[loc] * xi[i] + intercept[loc];
367         } else {
368           yi[i] = y[loc];
369         }
370       }
371     }
372 
373     return yi;
374   }
375 
376 }