List of usage examples for java.lang Math atan2
@HotSpotIntrinsicCandidate public static double atan2(double y, double x)
From source file:de.biomedical_imaging.ij.steger.Link.java
public void compute_contours(byte[] ismax, float[] eigval, float[] normx, float[] normy, float[] posx, float[] posy, float[] gradx, float[] grady, Lines contours, MutableInt num_result, double sigma, boolean extend_lines, int mode, double low, double high, int width, int height, Junctions junctions) { int i = 0, j = 0, k, l, it, pos, nextpos, nexti; int begin, end; int x, y;// w ww .ja va2 s . c om int octant, last_octant; int[] label; int num_cont, num_pnt; int size_cont, size_pnt; float[] row, col, trow, tcol; float[] angle, tangle; float[] resp, tresp; Junction[] junc; int num_junc, size_junc; Line[] cont; Line tmp_cont; LinesUtil.contour_class cls; double max; int maxx, maxy; int nextx, nexty; double nx, ny; double alpha, nextalpha, diff, mindiff, dist, mindist; double beta, last_beta, diff1, diff2; double px, py, nextpx = 0, nextpy = 0; double dx, dy; float tmp; int area; int[] indx; int indx_max; boolean nextismax; Crossref[] cross; int m = 0, max_line, num_add; MutableInt num_line = new MutableInt(); double length, response; Offset[] line; double mx, my, gx, gy, s, end_angle = 0, end_resp = 0; MutableDouble t = new MutableDouble(); float[] extx, exty; boolean add_ext; Region seg = new Region(); Chord[] rl; Width w = new Width(); /* * The image label contains information on the pixels that have been * processed by the linking algorithm. */ label = new int[(width * height)]; java.util.Arrays.fill(label, 0); // Vermutlich nicht notwendig, da // standardmig 0. /* * The image indx is an index into the table of all pixels that possibly * could be starting points for new lines. It is used to quickly * determine the next starting point of a line. */ indx = new int[(width * height)]; java.util.Arrays.fill(indx, 0); num_cont = 0; num_junc = 0; size_cont = LinesUtil.INITIAL_SIZE; size_pnt = LinesUtil.INITIAL_SIZE; size_junc = LinesUtil.INITIAL_SIZE; cont = new Line[size_cont]; for (int o = 0; o < cont.length; o++) { cont[o] = new Line(); } row = new float[size_pnt]; col = new float[size_pnt]; angle = new float[size_pnt]; resp = new float[size_pnt]; junc = new Junction[size_junc]; for (int o = 0; o < junc.length; o++) { junc[o] = new Junction(); } /* Select all pixels that can be starting points for lines. */ Threshold.threshold(ismax, 2, width, height, seg); /* Count the number of possible starting points. */ area = 0; for (i = 0; i < seg.num; i++) area += seg.rl[i].ce - seg.rl[i].cb + 1; /* Create the index of possible starting points. */ cross = new Crossref[area]; for (int o = 0; o < cross.length; o++) { cross[o] = new Crossref(); } k = 0; rl = seg.rl; for (i = 0; i < seg.num; i++) { x = rl[i].r; for (y = rl[i].cb; y <= rl[i].ce; y++) { pos = LinesUtil.LINCOOR(x, y, width); cross[k].x = (short) x; cross[k].y = (short) y; cross[k].value = eigval[pos]; cross[k].done = false; k++; } } java.util.Arrays.sort(cross); // qsort(cross,area,sizeof(*cross),compare_crossrefs); for (i = 0; i < area; i++) indx[LinesUtil.LINCOOR(cross[i].x, cross[i].y, width)] = i + 1; /* Link lines points. */ indx_max = 0; for (;;) { /* * Contour class unknown at this point; therefore assume both ends * free. */ cls = LinesUtil.contour_class.cont_no_junc; /* Search for next starting point. */ while (indx_max < area && cross[indx_max].done) indx_max++; /* Stop if no feasible starting point exists. */ if (indx_max == area) break; max = cross[indx_max].value; maxx = cross[indx_max].x; maxy = cross[indx_max].y; if (max == 0.0) break; /* Add starting point to the line. */ num_pnt = 0; pos = LinesUtil.LINCOOR(maxx, maxy, width); label[pos] = (num_cont + 1); if (!(indx[pos] == 0)) cross[(indx[pos] - 1)].done = true; row[num_pnt] = posx[pos]; col[num_pnt] = posy[pos]; /* Select line direction. */ nx = -normy[pos]; ny = normx[pos]; alpha = Math.atan2(ny, nx); if (alpha < 0.0) alpha += 2.0 * Math.PI; if (alpha >= Math.PI) alpha -= Math.PI; octant = (int) (Math.floor(4.0 / Math.PI * alpha + 0.5)) % 4; /* * Select normal to the line. The normal points to the right of the * line as the line is traversed from 0 to num-1. Since the points * are sorted in reverse order before the second iteration, the * first beta actually has to point to the left of the line! */ beta = alpha + Math.PI / 2.0; if (beta >= 2.0 * Math.PI) beta -= 2.0 * Math.PI; angle[num_pnt] = (float) beta; resp[num_pnt] = (float) interpolate_response(eigval, maxx, maxy, posx[pos], posy[pos], width, height); num_pnt++; /* Mark double responses as processed. */ for (i = 0; i < 2; i++) { nextx = maxx + cleartab[octant][i][0]; nexty = maxy + cleartab[octant][i][1]; if (nextx < 0 || nextx >= height || nexty < 0 || nexty >= width) continue; nextpos = LinesUtil.LINCOOR(nextx, nexty, width); if (ismax[nextpos] > 0) { nx = -normy[nextpos]; ny = normx[nextpos]; nextalpha = Math.atan2(ny, nx); if (nextalpha < 0.0) nextalpha += 2.0 * Math.PI; if (nextalpha >= Math.PI) nextalpha -= Math.PI; diff = Math.abs(alpha - nextalpha); if (diff >= Math.PI / 2.0) diff = Math.PI - diff; if (diff < MAX_ANGLE_DIFFERENCE) { label[nextpos] = (num_cont + 1); if (!(indx[nextpos] == 0)) cross[(indx[nextpos] - 1)].done = true; } } } for (it = 1; it <= 2; it++) { if (it == 1) { /* * Search along the initial line direction in the first * iteration. */ x = maxx; y = maxy; pos = LinesUtil.LINCOOR(x, y, width); nx = -normy[pos]; ny = normx[pos]; alpha = Math.atan2(ny, nx); if (alpha < 0.0) alpha += 2.0 * Math.PI; if (alpha >= Math.PI) alpha -= Math.PI; last_octant = (int) (Math.floor(4.0 / Math.PI * alpha + 0.5)) % 4; last_beta = alpha + Math.PI / 2.0; if (last_beta >= 2.0 * Math.PI) last_beta -= 2.0 * Math.PI; } else { /* Search in the opposite direction in the second iteration. */ x = maxx; y = maxy; pos = LinesUtil.LINCOOR(x, y, width); nx = -normy[pos]; ny = normx[pos]; alpha = Math.atan2(ny, nx); if (alpha < 0.0) alpha += 2.0 * Math.PI; if (alpha >= Math.PI) alpha -= Math.PI; last_octant = (int) (Math.floor(4.0 / Math.PI * alpha + 0.5)) % 4 + 4; last_beta = alpha + Math.PI / 2.0; if (last_beta >= 2.0 * Math.PI) last_beta -= 2.0 * Math.PI; } if (it == 2) { /* Sort the points found in the first iteration in reverse. */ for (i = 0; i < num_pnt / 2; i++) { tmp = row[i]; row[i] = row[(num_pnt - 1 - i)]; row[(num_pnt - 1 - i)] = tmp; tmp = col[i]; col[i] = col[(num_pnt - 1 - i)]; col[(num_pnt - 1 - i)] = tmp; tmp = angle[i]; angle[i] = angle[(num_pnt - 1 - i)]; angle[(num_pnt - 1 - i)] = tmp; tmp = resp[i]; resp[i] = resp[(num_pnt - 1 - i)]; resp[(num_pnt - 1 - i)] = tmp; } } /* Now start adding appropriate neighbors to the line. */ for (;;) { pos = LinesUtil.LINCOOR(x, y, width); nx = -normy[pos]; ny = normx[pos]; px = posx[pos]; py = posy[pos]; /* Orient line direction w.r.t. the last line direction. */ alpha = Math.atan2(ny, nx); if (alpha < 0.0) alpha += 2.0 * Math.PI; if (alpha >= Math.PI) alpha -= Math.PI; octant = (int) (Math.floor(4.0 / Math.PI * alpha + 0.5)) % 4; switch (octant) { case 0: if (last_octant >= 3 && last_octant <= 5) octant = 4; break; case 1: if (last_octant >= 4 && last_octant <= 6) octant = 5; break; case 2: if (last_octant >= 4 && last_octant <= 7) octant = 6; break; case 3: if (last_octant == 0 || last_octant >= 6) octant = 7; break; } last_octant = octant; /* Determine appropriate neighbor. */ nextismax = false; nexti = 1; mindiff = Double.MAX_VALUE; for (i = 0; i < 3; i++) { nextx = x + dirtab[octant][i][0]; nexty = y + dirtab[octant][i][1]; if (nextx < 0 || nextx >= height || nexty < 0 || nexty >= width) continue; nextpos = LinesUtil.LINCOOR(nextx, nexty, width); if (ismax[nextpos] == 0) continue; nextpx = posx[nextpos]; nextpy = posy[nextpos]; dx = nextpx - px; dy = nextpy - py; dist = Math.sqrt(dx * dx + dy * dy); nx = -normy[nextpos]; ny = normx[nextpos]; nextalpha = Math.atan2(ny, nx); if (nextalpha < 0.0) nextalpha += 2.0 * Math.PI; if (nextalpha >= Math.PI) nextalpha -= Math.PI; diff = Math.abs(alpha - nextalpha); if (diff >= Math.PI / 2.0) diff = Math.PI - diff; diff = dist + diff; if (diff < mindiff) { mindiff = diff; nexti = i; } if (!(ismax[nextpos] == 0)) nextismax = true; } /* Mark double responses as processed. */ for (i = 0; i < 2; i++) { nextx = x + cleartab[octant][i][0]; nexty = y + cleartab[octant][i][1]; if (nextx < 0 || nextx >= height || nexty < 0 || nexty >= width) continue; nextpos = LinesUtil.LINCOOR(nextx, nexty, width); if (ismax[nextpos] > 0) { nx = -normy[nextpos]; ny = normx[nextpos]; nextalpha = Math.atan2(ny, nx); if (nextalpha < 0.0) nextalpha += 2.0 * Math.PI; if (nextalpha >= Math.PI) nextalpha -= Math.PI; diff = Math.abs(alpha - nextalpha); if (diff >= Math.PI / 2.0) diff = Math.PI - diff; if (diff < MAX_ANGLE_DIFFERENCE) { label[nextpos] = (num_cont + 1); if (!(indx[nextpos] == 0)) cross[(indx[nextpos] - 1)].done = true; } } } /* Have we found the end of the line? */ if (!nextismax) break; /* If not, add the neighbor to the line. */ x += dirtab[octant][nexti][0]; y += dirtab[octant][nexti][1]; if (num_pnt >= size_pnt) { size_pnt = (int) Math.floor((double) (size_pnt * LinesUtil.REALLOC_FACTOR)); float[] newArr = new float[size_pnt]; for (int o = 0; o < row.length; o++) { newArr[o] = row[o]; } row = newArr; newArr = new float[size_pnt]; for (int o = 0; o < col.length; o++) { newArr[o] = col[o]; } col = newArr; newArr = new float[size_pnt]; for (int o = 0; o < angle.length; o++) { newArr[o] = angle[o]; } angle = newArr; newArr = new float[size_pnt]; for (int o = 0; o < resp.length; o++) { resp[o] = resp[o]; } resp = newArr; } pos = LinesUtil.LINCOOR(x, y, width); row[num_pnt] = posx[pos]; col[num_pnt] = posy[pos]; /* * Orient normal to the line direction w.r.t. the last * normal. */ nx = normx[pos]; ny = normy[pos]; beta = Math.atan2(ny, nx); if (beta < 0.0) beta += 2.0 * Math.PI; if (beta >= Math.PI) beta -= Math.PI; diff1 = Math.abs(beta - last_beta); if (diff1 >= Math.PI) diff1 = 2.0 * Math.PI - diff1; diff2 = Math.abs(beta + Math.PI - last_beta); if (diff2 >= Math.PI) diff2 = 2.0 * Math.PI - diff2; if (diff1 < diff2) { angle[num_pnt] = (float) beta; last_beta = beta; } else { angle[num_pnt] = (float) (beta + Math.PI); last_beta = beta + Math.PI; } resp[num_pnt] = (float) interpolate_response(eigval, x, y, posx[pos], posy[pos], width, height); num_pnt++; /* * If the appropriate neighbor is already processed a * junction point is found. */ if (label[pos] > 0) { if (num_junc >= size_junc) { size_junc = (int) Math.floor((double) (size_junc * LinesUtil.REALLOC_FACTOR)); Junction[] junch = new Junction[size_junc]; for (int o = 0; o < junch.length; o++) { if (o < junc.length) junch[o] = junc[o]; else junch[o] = new Junction(); } junc = junch; } /* Look for the junction point in the other line. */ k = label[pos] - 1; if (k == num_cont) { /* Line intersects itself. */ for (j = 0; j < num_pnt - 1; j++) { if (row[j] == posx[pos] && col[j] == posy[pos]) { if (j == 0) { /* Contour is closed. */ cls = LinesUtil.contour_class.cont_closed; for (i = 0; i < num_pnt / 2; i++) { tmp = row[i]; row[i] = row[(num_pnt - 1 - i)]; row[(num_pnt - 1 - i)] = tmp; tmp = col[i]; col[i] = col[(num_pnt - 1 - i)]; col[(num_pnt - 1 - i)] = tmp; tmp = angle[i]; angle[i] = angle[(num_pnt - 1 - i)]; angle[(num_pnt - 1 - i)] = tmp; tmp = resp[i]; resp[i] = resp[(num_pnt - 1 - i)]; resp[(num_pnt - 1 - i)] = tmp; } it = 2; } else { if (it == 2) { /* Determine contour class. */ if (cls == LinesUtil.contour_class.cont_start_junc) cls = LinesUtil.contour_class.cont_both_junc; else cls = LinesUtil.contour_class.cont_end_junc; /* Index j is the correct index. */ junc[num_junc].cont1 = num_cont; junc[num_junc].cont2 = num_cont; junc[num_junc].pos = j; junc[num_junc].x = posx[pos]; junc[num_junc].y = posy[pos]; num_junc++; } else { /* Determine contour class. */ cls = LinesUtil.contour_class.cont_start_junc; /* * Index num_pnt-1-j is the correct * index since the line is going to * be sorted in reverse. */ junc[num_junc].cont1 = num_cont; junc[num_junc].cont2 = num_cont; junc[num_junc].pos = num_pnt - 1 - j; junc[num_junc].x = posx[pos]; junc[num_junc].y = posy[pos]; num_junc++; } } break; } } /* * Mark this case as being processed for the * algorithm below. */ j = -1; } else { for (j = 0; j < cont[k].num; j++) { if (cont[k].row[j] == posx[pos] && cont[k].col[j] == posy[pos]) break; } /* * If no point can be found on the other line a * double response must have occured. In this case, * find the nearest point on the other line and add * it to the current line. */ if (j == cont[k].num) { mindist = Double.MAX_VALUE; j = -1; for (l = 0; l < cont[k].num; l++) { dx = posx[pos] - cont[k].row[l]; dy = posy[pos] - cont[k].col[l]; dist = Math.sqrt(dx * dx + dy * dy); if (dist < mindist) { mindist = dist; j = l; } } /* * Add the point with index j to the current * line. */ if (num_pnt >= size_pnt) { size_pnt = (int) Math.floor((double) (size_pnt * LinesUtil.REALLOC_FACTOR)); float[] newArr = new float[size_pnt]; for (int o = 0; o < row.length; o++) { newArr[o] = row[o]; } row = newArr; newArr = new float[size_pnt]; for (int o = 0; o < col.length; o++) { newArr[o] = col[o]; } col = newArr; newArr = new float[size_pnt]; for (int o = 0; o < angle.length; o++) { newArr[o] = angle[o]; } angle = newArr; newArr = new float[size_pnt]; for (int o = 0; o < resp.length; o++) { resp[o] = resp[o]; } resp = newArr; } row[num_pnt] = cont[k].row[j]; col[num_pnt] = cont[k].col[j]; beta = cont[k].angle[j]; if (beta >= Math.PI) beta -= Math.PI; diff1 = Math.abs(beta - last_beta); if (diff1 >= Math.PI) diff1 = 2.0 * Math.PI - diff1; diff2 = Math.abs(beta + Math.PI - last_beta); if (diff2 >= Math.PI) diff2 = 2.0 * Math.PI - diff2; if (diff1 < diff2) angle[num_pnt] = (float) beta; else angle[num_pnt] = (float) (beta + Math.PI); resp[num_pnt] = cont[k].response[j]; num_pnt++; } } /* * Add the junction point only if it is not one of the * other line's endpoints. */ if (j > 0 && j < cont[k].num - 1) { /* Determine contour class. */ if (it == 1) cls = LinesUtil.contour_class.cont_start_junc; else if (cls == LinesUtil.contour_class.cont_start_junc) cls = LinesUtil.contour_class.cont_both_junc; else cls = LinesUtil.contour_class.cont_end_junc; /* Add the new junction. */ junc[num_junc].cont1 = k; junc[num_junc].cont2 = num_cont; junc[num_junc].pos = j; junc[num_junc].x = row[(num_pnt - 1)]; junc[num_junc].y = col[(num_pnt - 1)]; num_junc++; } break; } label[pos] = (num_cont + 1); if (!(indx[pos] == 0)) cross[(indx[pos] - 1)].done = true; } } if (num_pnt > 1) { /* Only add lines with at least two points. */ if (num_cont >= size_cont) { size_cont = (int) Math.floor((double) (size_cont * LinesUtil.REALLOC_FACTOR)); Line[] conth = new Line[size_cont]; for (int o = 0; o < conth.length; o++) { // true ? (conth[o] = cont[0]) : (conth[o] = new // contour()); if (o < cont.length) conth[o] = cont[o]; else conth[o] = new Line(); } cont = conth; } cont[num_cont] = new Line(); cont[num_cont].row = java.util.Arrays.copyOf(row, num_pnt); cont[num_cont].col = java.util.Arrays.copyOf(col, num_pnt); cont[num_cont].angle = java.util.Arrays.copyOf(angle, num_pnt); cont[num_cont].response = java.util.Arrays.copyOf(resp, num_pnt); cont[num_cont].width_r = null; cont[num_cont].width_l = null; cont[num_cont].asymmetry = null; cont[num_cont].intensity = null; cont[num_cont].num = num_pnt; cont[num_cont].setContourClass(cls); num_cont++; } else { /* * Delete the point from the label image; we can use maxx and * maxy as the coordinates in the label image in this case. */ for (i = -1; i <= 1; i++) { for (j = -1; j <= 1; j++) { pos = LinesUtil.LINCOOR(LinesUtil.BR(maxx + i, height), LinesUtil.BC(maxy + j, width), width); if (label[pos] == num_cont + 1) label[pos] = 0; } } } } /* * Now try to extend the lines at their ends to find additional * junctions. */ if (extend_lines) { /* Sign by which the gradient has to be multiplied below. */ if (mode == LinesUtil.MODE_LIGHT) s = 1; else s = -1; double MAX_LINE_EXTENSION = 2.5 * sigma; length = MAX_LINE_EXTENSION; max_line = (int) Math.ceil(length * 3); line = new Offset[max_line]; for (int o = 0; o < line.length; o++) { line[o] = new Offset(); } extx = new float[max_line]; exty = new float[max_line]; for (i = 0; i < num_cont; i++) { tmp_cont = cont[i]; num_pnt = tmp_cont.num; if (num_pnt == 1) continue; if (tmp_cont.getContourClass() == LinesUtil.contour_class.cont_closed) continue; trow = tmp_cont.row; tcol = tmp_cont.col; tangle = tmp_cont.angle; tresp = tmp_cont.response; /* Check both ends of the line (it==-1: start, it==1: end). */ for (it = -1; it <= 1; it += 2) { /* * Determine the direction of the search line. This is done * by using the normal to the line (angle). Since this * normal may point to the left of the line (see below) we * have to check for this case by comparing the normal to * the direction of the line at its respective end point. */ if (it == -1) { /* Start point of the line. */ if (tmp_cont.getContourClass() == LinesUtil.contour_class.cont_start_junc || tmp_cont.getContourClass() == LinesUtil.contour_class.cont_both_junc) continue; dx = trow[1] - trow[0]; dy = tcol[1] - tcol[0]; alpha = tangle[0]; nx = Math.cos(alpha); ny = Math.sin(alpha); if (nx * dy - ny * dx < 0) { /* Turn the normal by +90 degrees. */ mx = -ny; my = nx; } else { /* Turn the normal by -90 degrees. */ mx = ny; my = -nx; } px = trow[0]; py = tcol[0]; response = tresp[0]; } else { /* End point of the line. */ if (tmp_cont.getContourClass() == LinesUtil.contour_class.cont_end_junc || tmp_cont.getContourClass() == LinesUtil.contour_class.cont_both_junc) continue; dx = trow[(num_pnt - 1)] - trow[(num_pnt - 2)]; dy = tcol[(num_pnt - 1)] - tcol[(num_pnt - 2)]; alpha = tangle[(num_pnt - 1)]; nx = Math.cos(alpha); ny = Math.sin(alpha); if (nx * dy - ny * dx < 0) { /* Turn the normal by -90 degrees. */ mx = ny; my = -nx; } else { /* Turn the normal by +90 degrees. */ mx = -ny; my = nx; } px = trow[(num_pnt - 1)]; py = tcol[(num_pnt - 1)]; response = tresp[(num_pnt - 1)]; } /* * Determine the current pixel and calculate the pixels on * the search line. */ x = (int) Math.floor(px + 0.5); y = (int) Math.floor(py + 0.5); dx = px - x; dy = py - y; w.bresenham(mx, my, dx, dy, length, line, num_line); /* * Now determine whether we can go only uphill (bright * lines) or downhill (dark lines) until we hit another * line. */ num_add = 0; add_ext = false; for (k = 0; k < num_line.intValue(); k++) { nextx = x + line[k].x; nexty = y + line[k].y; MutableDouble hnextpx = new MutableDouble(nextpx); MutableDouble hnextpy = new MutableDouble(nextpy); closest_point(px, py, mx, my, (double) nextx, (double) nexty, hnextpx, hnextpy, t); nextpx = hnextpx.getValue(); nextpy = hnextpy.getValue(); /* * Ignore points before or less than half a pixel away * from the true end point of the line. */ if (t.getValue() <= 0.5) continue; /* * Stop if the gradient can't be interpolated any more * or if the next point lies outside the image. */ if (nextpx < 0 || nextpy < 0 || nextpx >= height - 1 || nextpy >= width - 1 || nextx < 0 || nexty < 0 || nextx >= height || nexty >= width) break; MutableDouble hgx = new MutableDouble(); MutableDouble hgy = new MutableDouble(); interpolate_gradient(gradx, grady, nextpx, nextpy, width, hgx, hgy); gx = hgx.getValue(); gy = hgy.getValue(); /* * Stop if we can't go uphill anymore. This is * determined by the dot product of the line direction * and the gradient. If it is smaller than 0 we go * downhill (reverse for dark lines). */ nextpos = LinesUtil.LINCOOR(nextx, nexty, width); if (s * (mx * gx + my * gy) < 0 && label[nextpos] == 0) break; /* Have we hit another line? */ if (label[nextpos] > 0) { m = label[nextpos] - 1; /* Search for the junction point on the other line. */ mindist = Double.MAX_VALUE; j = -1; for (l = 0; l < cont[m].num; l++) { dx = nextpx - cont[m].row[l]; dy = nextpy - cont[m].col[l]; dist = Math.sqrt(dx * dx + dy * dy); if (dist < mindist) { mindist = dist; j = l; } } /* * This should not happen... But better safe than * sorry... */ if (mindist > 3.0) { break; } extx[num_add] = cont[m].row[j]; exty[num_add] = cont[m].col[j]; end_resp = cont[m].response[j]; end_angle = cont[m].angle[j]; beta = end_angle; if (beta >= Math.PI) beta -= Math.PI; diff1 = Math.abs(beta - alpha); if (diff1 >= Math.PI) diff1 = 2.0 * Math.PI - diff1; diff2 = Math.abs(beta + Math.PI - alpha); if (diff2 >= Math.PI) diff2 = 2.0 * Math.PI - diff2; if (diff1 < diff2) end_angle = beta; else end_angle = beta + Math.PI; num_add++; add_ext = true; break; } else { extx[num_add] = (float) nextpx; exty[num_add] = (float) nextpy; num_add++; } } if (add_ext) { /* Make room for the new points. */ num_pnt += num_add; float[] newArr = new float[num_pnt]; for (int o = 0; o < trow.length; o++) { newArr[o] = trow[o]; } trow = newArr; newArr = new float[num_pnt]; for (int o = 0; o < tcol.length; o++) { newArr[o] = tcol[o]; } tcol = newArr; newArr = new float[num_pnt]; for (int o = 0; o < tangle.length; o++) { newArr[o] = tangle[o]; } tangle = newArr; newArr = new float[num_pnt]; for (int o = 0; o < tresp.length; o++) { newArr[o] = tresp[o]; } tresp = newArr; tmp_cont.row = trow; tmp_cont.col = tcol; tmp_cont.angle = tangle; tmp_cont.response = tresp; tmp_cont.num = num_pnt; if (it == -1) { /* Move points on the line up num_add places. */ for (k = num_pnt - 1 - num_add; k >= 0; k--) { trow[(k + num_add)] = trow[k]; tcol[(k + num_add)] = tcol[k]; tangle[(k + num_add)] = tangle[k]; tresp[(k + num_add)] = tresp[k]; } /* Insert points at the beginning of the line. */ for (k = 0; k < num_add; k++) { trow[k] = extx[(num_add - 1 - k)]; tcol[k] = exty[(num_add - 1 - k)]; tangle[k] = (float) alpha; tresp[k] = (float) response; } tangle[0] = (float) end_angle; tresp[0] = (float) end_resp; /* Adapt indices of the previously found junctions. */ for (k = 0; k < num_junc; k++) { if (junc[k].cont1 == i) junc[k].pos += num_add; } } else { /* Insert points at the end of the line. */ for (k = 0; k < num_add; k++) { trow[(num_pnt - num_add + k)] = extx[k]; tcol[(num_pnt - num_add + k)] = exty[k]; tangle[(num_pnt - num_add + k)] = (float) alpha; tresp[(num_pnt - num_add + k)] = (float) response; } tangle[(num_pnt - 1)] = (float) end_angle; tresp[(num_pnt - 1)] = (float) end_resp; } /* If necessary, make room for the new junction. */ if (num_junc >= size_junc) { size_junc = (int) Math.floor((double) (size_junc * LinesUtil.REALLOC_FACTOR)); Junction[] junch = new Junction[size_junc]; for (int o = 0; o < junch.length; o++) { if (o < junc.length) junch[o] = junc[o]; else junch[o] = new Junction(); } junc = junch; } /* * Add the junction point only if it is not one of the * other line's endpoints. */ if (j > 0 && j < cont[m].num - 1) { if (it == -1) { if (tmp_cont.getContourClass() == LinesUtil.contour_class.cont_end_junc) tmp_cont.setContourClass(LinesUtil.contour_class.cont_both_junc); else tmp_cont.setContourClass(LinesUtil.contour_class.cont_start_junc); } else { if (tmp_cont.getContourClass() == LinesUtil.contour_class.cont_start_junc) tmp_cont.setContourClass(LinesUtil.contour_class.cont_both_junc); else tmp_cont.setContourClass(LinesUtil.contour_class.cont_end_junc); } junc[num_junc].cont1 = m; junc[num_junc].cont2 = i; junc[num_junc].pos = j; if (it == -1) { junc[num_junc].x = trow[0]; junc[num_junc].y = tcol[0]; } else { junc[num_junc].x = trow[(num_pnt - 1)]; junc[num_junc].y = tcol[(num_pnt - 1)]; } num_junc++; } } } } } /* Done with linking. Now split the lines at the junction points. */ java.util.Arrays.sort(junc); for (i = 0; i < num_junc; i += k) { j = junc[i].cont1; tmp_cont = cont[j]; num_pnt = tmp_cont.num; /* Count how often line j needs to be split. */ for (k = 0; junc[(i + k)].cont1 == j && i + k < num_junc; k++) ; if (k == 1 && tmp_cont.row[0] == tmp_cont.row[(num_pnt - 1)] && tmp_cont.col[0] == tmp_cont.col[(num_pnt - 1)]) { /* * If only one junction point is found and the line is closed it * only needs to be rearranged cyclically, but not split. */ begin = junc[i].pos; trow = tmp_cont.row; tcol = tmp_cont.col; tangle = tmp_cont.angle; tresp = tmp_cont.response; tmp_cont.row = new float[num_pnt]; tmp_cont.col = new float[num_pnt]; tmp_cont.angle = new float[num_pnt]; tmp_cont.response = new float[num_pnt]; for (l = 0; l < num_pnt; l++) { pos = begin + l; /* Skip starting point so that it is not added twice. */ if (pos >= num_pnt) pos = begin + l - num_pnt + 1; tmp_cont.row[l] = trow[pos]; tmp_cont.col[l] = tcol[pos]; tmp_cont.angle[l] = tangle[pos]; tmp_cont.response[l] = tresp[pos]; } /* Modify contour class. */ tmp_cont.setContourClass(LinesUtil.contour_class.cont_both_junc); } else { /* Otherwise the line has to be split. */ for (l = 0; l <= k; l++) { if (l == 0) begin = 0; else begin = junc[(i + l - 1)].pos; if (l == k) end = tmp_cont.num - 1; else end = junc[(i + l)].pos; num_pnt = end - begin + 1; if (num_pnt == 1 && k > 1) { /* Do not add one point segments. */ continue; } if (num_cont >= size_cont) { size_cont = (int) Math.floor((double) (size_cont * LinesUtil.REALLOC_FACTOR)); Line[] conth = new Line[size_cont]; for (int o = 0; o < cont.length; o++) { conth[o] = cont[o]; } cont = conth; } cont[num_cont] = new Line(); cont[num_cont].row = new float[num_pnt]; cont[num_cont].col = new float[num_pnt]; cont[num_cont].angle = new float[num_pnt]; cont[num_cont].response = new float[num_pnt]; System.arraycopy(tmp_cont.row, begin, cont[num_cont].row, 0, num_pnt); System.arraycopy(tmp_cont.col, begin, cont[num_cont].col, 0, num_pnt); System.arraycopy(tmp_cont.angle, begin, cont[num_cont].angle, 0, num_pnt); System.arraycopy(tmp_cont.response, begin, cont[num_cont].response, 0, num_pnt); cont[num_cont].width_r = null; cont[num_cont].width_l = null; cont[num_cont].asymmetry = null; cont[num_cont].intensity = null; cont[num_cont].num = num_pnt; /* Modify contour class. */ if (l == 0) { if (tmp_cont.getContourClass() == LinesUtil.contour_class.cont_start_junc || tmp_cont.getContourClass() == LinesUtil.contour_class.cont_both_junc) cont[num_cont].setContourClass(LinesUtil.contour_class.cont_both_junc); else cont[num_cont].setContourClass(LinesUtil.contour_class.cont_end_junc); } else if (l == k) { if (tmp_cont.getContourClass() == LinesUtil.contour_class.cont_end_junc || tmp_cont.getContourClass() == LinesUtil.contour_class.cont_both_junc) cont[num_cont].setContourClass(LinesUtil.contour_class.cont_both_junc); else cont[num_cont].setContourClass(LinesUtil.contour_class.cont_start_junc); } else { cont[num_cont].setContourClass(LinesUtil.contour_class.cont_both_junc); } num_cont++; } cont[j] = cont[--num_cont]; } } /* Finally, check whether all angles point to the right of the line. */ for (i = 0; i < num_cont; i++) { tmp_cont = cont[i]; num_pnt = tmp_cont.num; if (num_pnt > 1) { trow = tmp_cont.row; tcol = tmp_cont.col; tangle = tmp_cont.angle; /* * One point of the contour is enough to determine the * orientation. */ k = (num_pnt - 1) / 2; /* * The next few lines are ok because lines have at least two * points. */ dx = trow[(k + 1)] - trow[k]; dy = tcol[(k + 1)] - tcol[k]; nx = Math.cos(tangle[k]); ny = Math.sin(tangle[k]); /* * If the angles point to the left of the line they have to be * adapted. The orientation is determined by looking at the * z-component of the cross-product of (dx,dy,0) and (nx,ny,0). */ if (nx * dy - ny * dx < 0) { for (j = 0; j < num_pnt; j++) { tangle[j] += Math.PI; if (tangle[j] >= 2 * Math.PI) tangle[j] -= 2 * Math.PI; } } } } for (Line c : cont) { if (c != null && c.getContourClass() != null) { c.setFrame(contours.getFrame()); contours.add(c); } } for (Junction jun : junc) { if (jun != null && !(jun.cont1 == 0 && jun.cont2 == 0)) { junctions.add(jun); } } num_result.setValue(num_cont); }
From source file:pl.edu.icm.visnow.lib.utils.ImageUtilities.java
public static BufferedImage cylindricalMapping(BufferedImage img, double f) { if (img == null) { return null; }/*from w w w .ja v a 2 s . c om*/ int w = img.getWidth(); int h = img.getHeight(); BufferedImage out = new BufferedImage(w, h, img.getType()); //System.out.println("w:"+w+", h:"+h); int x0 = (int) Math.floor(w / 2) + 1; int y0 = (int) Math.floor(h / 2) + 1; double tmax = Math.atan2((double) (w - x0), f); double tmin = Math.atan2(-((double) x0), f); double tstep = (tmax - tmin) / ((double) w); double vmax = ((double) (h - y0)) / f; double vmin = (-(double) y0) / f; double vstep = (vmax - vmin) / ((double) h); double theta, tan, cos; int x, y; for (int t = 0; t < w; t++) { theta = tmin + (double) t * tstep; tan = Math.tan(theta); cos = Math.cos(theta); x = (int) Math.round(f * tan) + x0; for (int v = 0; v < h; v++) { //nearest neighbour--------------------------------------- //x = (int)Math.round(f*tan) + x0; y = (int) Math.round((vmin + (double) v * vstep) * f / cos) + y0; if (x >= 0 && y >= 0 && x < w && y < h) { //piksel nowy x,y = piksel stary xd,yd out.setRGB(t, v, img.getRGB(x, y)); } //--------------------------------------------------------- } } return out; }
From source file:chat.client.gui.ChatActivity.java
private static double getDistance(double lat1, double lon1, double lat2, double lon2) { final double Radius = 6371 * 1E3; // Earth's mean radius double dLat = Math.toRadians(lat2 - lat1); double dLon = Math.toRadians(lon2 - lon1); double a = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.cos(Math.toRadians(lat1)) * Math.cos(Math.toRadians(lat2)) * Math.sin(dLon / 2) * Math.sin(dLon / 2); double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a)); return Radius * c; }
From source file:org.eclipse.smarthome.binding.astro.internal.calc.SunCalc.java
private double getAzimuth(double th, double a, double phi, double d) { double h = th - a; return Math.atan2(Math.sin(h), Math.cos(h) * Math.sin(phi) - Math.tan(d) * Math.cos(phi)); }
From source file:edu.uci.ics.jung.visualization.picking.ShapePickSupport.java
/** * Retrieves the shape template for <code>e</code> and * transforms it according to the positions of its endpoints * in <code>layout</code>./* w ww .j a va 2s. c om*/ * @param layout the <code>Layout</code> which specifies * <code>e</code>'s endpoints' positions * @param e the edge whose shape is to be returned * @return */ private Shape getTransformedEdgeShape(Layout<V, E> layout, E e) { Pair<V> pair = layout.getGraph().getEndpoints(e); V v1 = pair.getFirst(); V v2 = pair.getSecond(); boolean isLoop = v1.equals(v2); Point2D p1 = vv.getRenderContext().getMultiLayerTransformer().transform(Layer.LAYOUT, layout.transform(v1)); Point2D p2 = vv.getRenderContext().getMultiLayerTransformer().transform(Layer.LAYOUT, layout.transform(v2)); if (p1 == null || p2 == null) return null; float x1 = (float) p1.getX(); float y1 = (float) p1.getY(); float x2 = (float) p2.getX(); float y2 = (float) p2.getY(); // translate the edge to the starting vertex AffineTransform xform = AffineTransform.getTranslateInstance(x1, y1); Shape edgeShape = vv.getRenderContext().getEdgeShapeTransformer() .transform(Context.<Graph<V, E>, E>getInstance(vv.getGraphLayout().getGraph(), e)); if (isLoop) { // make the loops proportional to the size of the vertex Shape s2 = vv.getRenderContext().getVertexShapeTransformer().transform(v2); Rectangle2D s2Bounds = s2.getBounds2D(); xform.scale(s2Bounds.getWidth(), s2Bounds.getHeight()); // move the loop so that the nadir is centered in the vertex xform.translate(0, -edgeShape.getBounds2D().getHeight() / 2); } else { float dx = x2 - x1; float dy = y2 - y1; // rotate the edge to the angle between the vertices double theta = Math.atan2(dy, dx); xform.rotate(theta); // stretch the edge to span the distance between the vertices float dist = (float) Math.sqrt(dx * dx + dy * dy); xform.scale(dist, 1.0f); } // transform the edge to its location and dimensions edgeShape = xform.createTransformedShape(edgeShape); return edgeShape; }
From source file:org.gearvrf.keyboard.util.Util.java
public static float getZRotationAngle(GVRSceneObject rotatingObject, GVRTransform targetObject) { float angle = (float) Math .toDegrees(Math.atan2(targetObject.getPositionY() - rotatingObject.getTransform().getPositionY(), targetObject.getPositionX() - rotatingObject.getTransform().getPositionX())); return angle; }
From source file:org.freecine.filmscan.ScanStrip.java
/** Calculate the affine transform from scanStrip to a single frame (frame rotated to straight position, top left corner translated to (0,0) @param frame//from ww w .ja v a2s . c om @return */ AffineTransform getFrameXform(int frame) { /** Estimate film rotation from max 5 perforations */ int f1 = frame - 1; int f2 = frame + 1; int x1 = (f1 >= 0) ? perforations.get(f1).x : perforations.get(0).x; int x2 = (f2 < perforations.size()) ? perforations.get(f2).x : perforations.get(perforations.size() - 1).x; int y1 = (f1 >= 0) ? perforations.get(f1).y : perforations.get(0).y; int y2 = (f2 < perforations.size()) ? perforations.get(f2).y : perforations.get(perforations.size() - 1).y; double rot = Math.atan2((double) x2 - x1, (double) (y2 - y1)); // Translate the center of perforation to origin AffineTransform xform = new AffineTransform(); xform.translate(0, FRAME_HEIGHT / 2); xform.rotate(rot); xform.translate(-perforations.get(frame).x, -perforations.get(frame).y); // System.out.println( String.format( "frame %d: (%d, %d), rot %f", // frame,perforations.get(frame).x, -perforations.get(frame).y, rot )); return xform; }
From source file:uk.ac.diamond.scisoft.analysis.dataset.function.MapToPolarAndIntegrate.java
/** * This method uses applies weighting factors to every sampled data point to calculate integration profiles * //from w w w . j ava2s. c o m * @param datasets * input 2D dataset * @return 4 1D datasets for integral over radius, integral over azimuth (for given input and a uniform input) */ private List<Dataset> simple_value(IDataset... datasets) { if (datasets.length == 0) { return null; } List<Dataset> result = new ArrayList<Dataset>(); for (IDataset ids : datasets) { if (ids.getRank() != 2) { throw new IllegalArgumentException("operating on 2d arrays only"); } Dataset ds = DatasetUtils.convertToDataset(ids); int npts = (int) (erad - srad + 1); int apts = (int) (erad * (ephi - sphi) + 1); double dphi = (ephi - sphi) / apts; // Final intensity is 1D data float[] intensity = new float[npts]; float[] azimuth = new float[apts]; // Calculate bounding rectangle around the sector int nxstart = (int) Math.max(0, cx - erad); int nx = (int) Math.min(ds.getShapeRef()[1], cx + erad); int nystart = (int) Math.max(0, cy - erad); int ny = (int) Math.min(ds.getShapeRef()[0], cy + erad); for (int j = nystart; j < ny; j++) { for (int i = nxstart; i < nx; i++) { double Theta = Math.atan2((j - cy), (i - cx)); Theta = MathUtils.normalizeAngle(Theta, sphi + Math.PI); if ((Theta >= sphi) && (Theta <= ephi)) { double xR = i - cx; double yR = j - cy; double R = Math.sqrt(xR * xR + yR * yR); if ((R >= srad) && (R < erad)) { int k = (int) (R - srad); if (mask != null && !mask.getBoolean(j, i)) { continue; } double val = ds.getDouble(j, i); // Each point participating to the sector integration is weighted depending on // how far/close it is from the following point i+1 double fFac = (k + 1) - (R - srad); if (doRadial) { // Evaluates the intensity and the frequency intensity[k] += fFac * val; if (k < (npts - 1)) { intensity[k + 1] += (1 - fFac) * val; } } if (doAzimuthal) { double dk = 1. / R; int ak1 = Math.max(0, (int) ((Theta - dk / 2.0 - sphi) / dphi)); int ak2 = Math.min(apts - 1, (int) ((Theta + dk / 2.0 - sphi) / dphi)); for (int n = ak1; n <= ak2; n++) { fFac = ak2 - ak1 + 1.0; azimuth[n] += val / fFac; } } } } } } result.add(new FloatDataset(azimuth, apts)); result.add(new FloatDataset(intensity, npts)); } return result; }
From source file:org.orekit.tle.TLEPropagator.java
/** Retrieves the position and velocity. * @return the computed PVCoordinates.// w w w .j a v a2 s. c om * @exception OrekitException if current orbit is out of supported range * (too large eccentricity, too low perigee ...) */ private PVCoordinates computePVCoordinates() throws OrekitException { // Long period periodics final double axn = e * Math.cos(omega); double temp = 1.0 / (a * (1.0 - e * e)); final double xlcof = 0.125 * TLEConstants.A3OVK2 * sini0 * (3.0 + 5.0 * cosi0) / (1.0 + cosi0); final double aycof = 0.25 * TLEConstants.A3OVK2 * sini0; final double xll = temp * xlcof * axn; final double aynl = temp * aycof; final double xlt = xl + xll; final double ayn = e * Math.sin(omega) + aynl; final double elsq = axn * axn + ayn * ayn; final double capu = MathUtils.normalizeAngle(xlt - xnode, Math.PI); double epw = capu; double ecosE = 0; double esinE = 0; double sinEPW = 0; double cosEPW = 0; // Dundee changes: items dependent on cosio get recomputed: final double cosi0Sq = cosi0 * cosi0; final double x3thm1 = 3.0 * cosi0Sq - 1.0; final double x1mth2 = 1.0 - cosi0Sq; final double x7thm1 = 7.0 * cosi0Sq - 1.0; if (e > (1 - 1e-6)) { throw new OrekitException( "eccentricity becomes too large for TLE propagation " + "(e: {0}, satellite number: {1})", e, tle.getSatelliteNumber()); } if ((a * (1.0 - e) < 1.0) || (a * (1.0 + e) < 1.0)) { throw new OrekitException( "too small perigee radius for TLE propagation " + "(r: {0}, satellite number: {1})", a * (1. - e), tle.getSatelliteNumber()); } // Solve Kepler's' Equation. final double newtonRaphsonEpsilon = 1e-12; for (int j = 0; j < 10; j++) { boolean doSecondOrderNewtonRaphson = true; sinEPW = Math.sin(epw); cosEPW = Math.cos(epw); ecosE = axn * cosEPW + ayn * sinEPW; esinE = axn * sinEPW - ayn * cosEPW; final double f = capu - epw + esinE; if (Math.abs(f) < newtonRaphsonEpsilon) { break; } final double fdot = 1.0 - ecosE; double delta_epw = f / fdot; if (j == 0) { final double maxNewtonRaphson = 1.25 * Math.abs(e); doSecondOrderNewtonRaphson = false; if (delta_epw > maxNewtonRaphson) { delta_epw = maxNewtonRaphson; } else if (delta_epw < -maxNewtonRaphson) { delta_epw = -maxNewtonRaphson; } else { doSecondOrderNewtonRaphson = true; } } if (doSecondOrderNewtonRaphson) { delta_epw = f / (fdot + 0.5 * esinE * delta_epw); } epw += delta_epw; } // Short period preliminary quantities temp = 1.0 - elsq; final double pl = a * temp; final double r = a * (1.0 - ecosE); double temp2 = a / r; final double betal = Math.sqrt(temp); temp = esinE / (1.0 + betal); final double cosu = temp2 * (cosEPW - axn + ayn * temp); final double sinu = temp2 * (sinEPW - ayn - axn * temp); final double u = Math.atan2(sinu, cosu); final double sin2u = 2.0 * sinu * cosu; final double cos2u = 2.0 * cosu * cosu - 1.0; final double temp1 = TLEConstants.CK2 / pl; temp2 = temp1 / pl; // Update for short periodics final double rk = r * (1.0 - 1.5 * temp2 * betal * x3thm1) + 0.5 * temp1 * x1mth2 * cos2u; final double uk = u - 0.25 * temp2 * x7thm1 * sin2u; final double xnodek = xnode + 1.5 * temp2 * cosi0 * sin2u; final double xinck = i + 1.5 * temp2 * cosi0 * sini0 * cos2u; // Orientation vectors final double sinuk = Math.sin(uk); final double cosuk = Math.cos(uk); final double sinik = Math.sin(xinck); final double cosik = Math.cos(xinck); final double sinnok = Math.sin(xnodek); final double cosnok = Math.cos(xnodek); final double xmx = -sinnok * cosik; final double xmy = cosnok * cosik; final double ux = xmx * sinuk + cosnok * cosuk; final double uy = xmy * sinuk + sinnok * cosuk; final double uz = sinik * sinuk; // Position and velocity final double cr = 1000 * rk * TLEConstants.EARTH_RADIUS; final Vector3D pos = new Vector3D(cr * ux, cr * uy, cr * uz); final double rdot = TLEConstants.XKE * Math.sqrt(a) * esinE / r; final double rfdot = TLEConstants.XKE * Math.sqrt(pl) / r; final double xn = TLEConstants.XKE / (a * Math.sqrt(a)); final double rdotk = rdot - xn * temp1 * x1mth2 * sin2u; final double rfdotk = rfdot + xn * temp1 * (x1mth2 * cos2u + 1.5 * x3thm1); final double vx = xmx * cosuk - cosnok * sinuk; final double vy = xmy * cosuk - sinnok * sinuk; final double vz = sinik * cosuk; final double cv = 1000.0 * TLEConstants.EARTH_RADIUS / 60.0; final Vector3D vel = new Vector3D(cv * (rdotk * ux + rfdotk * vx), cv * (rdotk * uy + rfdotk * vy), cv * (rdotk * uz + rfdotk * vz)); return new PVCoordinates(pos, vel); }
From source file:pl.edu.icm.visnow.lib.utils.ImageUtilities.java
public static BufferedImage sphericalMapping(BufferedImage img, double f) { if (img == null) { return null; }//ww w. j av a 2 s . c o m int w = img.getWidth(); int h = img.getHeight(); BufferedImage out = new BufferedImage(w, h, img.getType()); //System.out.println("w:"+w+", h:"+h); int x0 = (int) Math.floor(w / 2) + 1; int y0 = (int) Math.floor(h / 2) + 1; double tmax = Math.atan2((double) (w - x0), f); double tmin = Math.atan2(-((double) x0), f); double tstep = (tmax - tmin) / ((double) w); double fimax = Math.atan2((double) (h - y0), Math.sqrt(f * f)); double fimin = Math.atan2(-(double) y0, Math.sqrt(f * f)); double fistep = (fimax - fimin) / ((double) h); //System.out.println("fimax:"+fimax+", fimin:"+fimin); double theta, tantheta, costheta, tanfi, phi; int x, y; for (int t = 0; t < w; t++) { theta = tmin + (double) t * tstep; tantheta = Math.tan(theta); x = (int) Math.round(f * tantheta) + x0; for (int fi = 0; fi < h; fi++) { //nearest neighbour--------------------------------------- phi = fimin + (double) fi * fistep; tanfi = Math.tan(phi); //x = (int)Math.round(f*tantheta) + x0; y = (int) Math.round(Math.sqrt((x - x0) * (x - x0) + f * f) * tanfi) + y0; if (x >= 0 && y >= 0 && x < w && y < h) { //piksel nowy x,y = piksel stary xd,yd out.setRGB(t, fi, img.getRGB(x, y)); } //--------------------------------------------------------- } } return out; }