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
25
26
27
28
29
30
31 public class PolarPlot {
32
33
34
35
36 static final Logger LOGGER = LoggerFactory.getLogger(PolarPlot.class.getName());
37 private STmap mapCell;
38 private Point2d gradientcoord;
39
40
41
42
43 protected Rectangle plotArea;
44
45
46
47 public double kscale;
48
49
50
51 public double uscale;
52
53
54
55
56 public boolean labels = false;
57
58
59
60
61
62
63
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
75
76
77
78
79 int[] getShift() {
80 int[] ret = new int[mapCell.getT()];
81 for (int f = 0; f < mapCell.getT(); f++) {
82 double dist = Double.MAX_VALUE;
83 for (int i = 0; i < mapCell.getRes(); i++) {
84 Point2d p = new Point2d(mapCell.getxMap()[f][i], mapCell.getyMap()[f][i]);
85 double disttmp = p.distance(gradientcoord);
86 if (disttmp < dist) {
87 dist = disttmp;
88 ret[f] = i;
89 }
90 }
91 }
92 return ret;
93 }
94
95
96
97
98
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
112
113
114
115
116
117
118
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;
125 for (int i = start; i < mapCell.getRes() + start; i++) {
126
127 int index = IPadArray.getIndex(mapCell.getRes(), i, IPadArray.CIRCULARPAD);
128
129 Point2d p = new Point2d(mapCell.getxMap()[f][index], mapCell.getyMap()[f][index]);
130 p.sub(mc);
131 ret[l++] = new Vector2d(p);
132 }
133 return ret;
134 }
135
136
137
138
139
140
141
142
143
144 double[] getRadius(int f, int shift, double[][] map) {
145 double[] ret = new double[map[f].length];
146 int l = 0;
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
156
157
158
159
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
168
169 }
170 return ret;
171
172 }
173
174
175
176
177
178
179
180
181
182 public void generatePlotFrame(String filename, int frame) throws IOException {
183 int[] shifts = getShift();
184
185 double[] magn = getRadius(frame, shifts[frame], mapCell.getMotMap());
186
187 double[] angles = getUniformAngles(magn.length);
188
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);
193 }
194
195 if (labels) {
196 polarPlotPoints(filename, angles, magn, new double[] { min, max });
197 } else {
198 polarPlotPoints(filename, angles, magn, null);
199 }
200
201 }
202
203
204
205
206
207
208
209
210
211
212
213 public void generatePlot(String filename) throws IOException {
214 int[] shifts = getShift();
215
216
217 double[][] magnF = new double[mapCell.getT()][];
218 for (int f = 0; f < mapCell.getT(); f++) {
219
220 magnF[f] = getRadius(f, shifts[f], mapCell.getMotMap());
221 }
222
223 double[] angles = getUniformAngles(magnF[0].length);
224
225 double[] magn = QuimPArrayUtils.getMeanC(magnF);
226
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);
231 }
232 if (labels) {
233 polarPlotPoints(filename, angles, magn, new double[] { min, max });
234 } else {
235 polarPlotPoints(filename, angles, magn, null);
236 }
237 LOGGER.debug("Polar plot saved: " + filename);
238
239 }
240
241
242
243
244
245
246
247
248
249
250
251 private void polarPlotPoints(String filename, double[] angles, double[] magn, double[] valMinMax)
252 throws IOException {
253
254 double plotscale = plotArea.getWidth() / 2 / QuimPArrayUtils.arrayMax(magn);
255 plotscale -= plotscale * uscale;
256
257 BufferedOutputStream out;
258 out = new BufferedOutputStream(new FileOutputStream(filename));
259 PrintWriter osw = new PrintWriter(out, true);
260 Rectangle extendedPlotArea = (Rectangle) plotArea.clone();
261 extendedPlotArea.grow(1, 1);
262 SVGwritter.writeHeader(osw, extendedPlotArea);
263
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
270
271 ret.remove(0);
272 ret.remove(ret.size() - 1);
273
274 qaxes.radiusLabels = ArrayUtils.toPrimitive(ret.toArray(new Double[0]));
275 }
276 qaxes.draw(osw);
277
278 for (int i = 0; i < angles.length; i++) {
279
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
285
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
298
299
300
301
302
303
304
305
306
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];
315 }
316 }
317 return angles;
318 }
319
320
321
322
323
324
325
326
327
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
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
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 }