| blob | b01aaf4fe51533f4944509dde69e3ffc40836fcd |
1 /**
3 TrakEM2 plugin for ImageJ(C).
4 Copyright (C) 2006, 2007 Albert Cardona.
6 This program is free software; you can redistribute it and/or
7 modify it under the terms of the GNU General Public License
8 as published by the Free Software Foundation (http://www.gnu.org/licenses/gpl.txt )
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19 You may contact Albert Cardona at acardona at ini.phys.ethz.ch
20 Institute of Neuroinformatics, University of Zurich / ETH, Switzerland.
21 **/
96 }
101 }
104 Utils.addEnablerListener((Checkbox)gd.getCheckboxes().get(0), new Component[]{(Component)gd.getCheckboxes().get(1)}, null);
123 }
124 }
126 // check that the calibration units are the same
130 if (!matchUnits(unit1, ref[i].getRootLayerSet().getCalibration().getUnit(), ref[i].getTitle())) {
132 }
133 }
134 }
136 return findSimilar(pipe, ref, ignore_orientation, ignore_calibration, mirror, chain_branches, true);
137 }
144 }
156 final String[] preset_names = new String[]{"X - 'medial lobe', Y - 'dorsal lobe', Z - 'peduncle'"};
160 final ArrayList<ZDisplayable> pipes_ref = all[iother].getRootLayerSet().getZDisplayables(Line3D.class, true);
164 // automatically find for the first preset
169 // check if none found
173 }
184 }
191 // refresh reference project choices
205 }
206 }
218 }
221 }
222 }
223 });
228 Utils.addEnablerListener((Checkbox)gd.getCheckboxes().get(0), new Component[]{(Component)gd.getNumericFields().get(0), (Component)gd.getNumericFields().get(1)}, null);
241 //////
246 // ok, ready
253 };
261 };
270 }
281 }
284 // check that the Calibration units are the same
285 if (!matchUnits(pipe.getLayerSet().getCalibration().getUnit(), all[iproject].getRootLayerSet().getCalibration().getUnit(), all[iproject].getTitle())) {
287 }
289 return findSimilarWithAxes(pipe, axes, axes_ref, pipes_ref, skip_ends, max_mut, min_chunk, transform_type, chain_branches, true, score_mut, substring_matching, direct, delta, distance_type);
290 }
292 /** Finds the first three X,Y,Z axes as specified by the names in preset. */
293 static private int[] findFirstXYZAxes(final String[] preset, final ArrayList<ZDisplayable> pipes, final String[] pipe_names) {
299 // Already all found, just filling names
300 next++;
302 }
307 next++;
308 }
310 }
312 /** Find the indices of the axes pipes in the pipes array list, and fills in the pipe_names array. */
313 static private int[] findXYZAxes(final Line3D[] axes, final ArrayList<ZDisplayable> pipes, final String[] pipe_names) {
322 next++;
323 }
325 }
327 static private void trySetAsAxis(final Line3D[] axes, final String[] preset, final ProjectThing pt) {
331 final String title = pipe.getProject().getShortMeaningfulTitle(pt, (ZDisplayable)pipe).toLowerCase();
332 //Utils.log2("title is " + title);
336 }
338 /** Finds the 3 pipes named according to the presets (such as "medial lobe", "dorsal lobe" and "peduncle"), that are structurally closest to the query_pipe in the Project Tree. In this fashion, if there are more than one hemisphere, the correct set of axes will be found for the query pipe.*/
343 }
345 // Step up level by level looking for pipes with the given preset names
358 }
360 }
363 }
365 /** Generate calibrated origin of coordinates. */
366 static public Object[] obtainOrigin(final Line3D[] axes, final int transform_type, final Vector3d[] o_ref) {
367 // pipe's axes
371 final Calibration cal = (null != axes[0].getLayerSet() ? axes[0].getLayerSet().getCalibration() : null);
372 // 1 - calibrate
375 }
376 // 2 - resample (although it's done before transforming, it's only for aesthetic purposes: it doesn't matter, won't ever go into dynamic programming machinery)
382 // return origin vectors for pipe's project
383 final Vector3d[] o = VectorString3D.createOrigin(vs[0], vs[1], vs[2], transform_type, o_ref); // requires resampled vs
386 }
388 /** Compare pipe to all pipes in pipes_ref, by first transforming to match both sets of axes. */
389 static public final Bureaucrat findSimilarWithAxes(final Line3D pipe, final Line3D[] axes, final Line3D[] axes_ref, final ArrayList<ZDisplayable> pipes_ref, final boolean skip_ends, final int max_mut, final float min_chunk, final int transform_type, final boolean chain_branches, final boolean show_gui, final boolean score_mut, final boolean substring_matching, final boolean direct, final double delta, final int sort_by_distance_type) {
399 }
404 }
406 // NEW style: bring query brain to reference brain space.
408 final Calibration cal1 = (null != pipe.getLayerSet() ? pipe.getLayerSet().getCalibration() : null);
411 // obtain origin vectors for reference project
418 // obtain axes origin vectors for pipe's project
425 // The transfer transform: from query axes to reference axes: M_ref * M_query^{-1}
429 // Transform the query axes only (the reference ones stay the same)
431 // axes are already calibrated
432 //Utils.log2("Axis " + i + " Length before transforming: " + vs_axes[i].computeLength());
434 //Utils.log2("Axis " + i + " Length after transforming: " + vs_axes[i].computeLength());
435 }
437 // If chain_branches, the query must also be split into as many branches as it generates.
438 // Should generate a tab for each potential branch? Or better, a tab but not with a table but with a list of labels, one per potential branch, and underneath as many tables in more tabs?
440 // the queries to do, according to how many different chains the query pipe is part of
446 // create all chained ref branches
447 chains_ref = createPipeChains(pipes_ref.get(0).getProject().getRootProjectThing(), pipes_ref.get(0).getLayerSet()); // TODO WARNING: if the parent has more than one pipe, this will query them all!
449 // add all possible query chains, starting at the parent of the chosen pipe
450 for (Chain chain : createPipeChains((ProjectThing)pipe.getProject().findProjectThing(pipe).getParent(), pipe.getLayerSet())) {
452 }
454 // no branching: single query of one single-pipe chain
457 // just add a single-pipe chain for each ref pipe
461 }
462 }
464 // set and calibrate them all
468 // each thread handles a ref pipe, which is to be matched against all queries
480 ////
483 // obtain a calibrated ref chain, to be uniquely processed by this thread
485 // match it against all queries
491 final Object[] ob = findBestMatch(vs1, vs2, delta1, skip_ends, max_mut, min_chunk, 1, direct, substring_matching);
493 //qh.addMatch(query, ref, ed, seq_sim, ed.getPhysicalDistance(skip_ends, max_mut, min_chunk));
496 1.0f
500 qm[next++].cm[k] = new ChainMatch(query, ref, ed, stats, prop_len, score(ed.getSimilarity(), ed.getDistance(), ed.getStatistics(skip_ends, max_mut, min_chunk, false)[3], Compare.W));
501 }
504 }
505 }
506 ////
507 }
508 });
509 }
513 // done!
514 // put result into the Worker
518 // add to the GUI (will sort them by phys_dist)
522 }
529 }
530 }
531 };
537 }
539 /** For a given pipe, create a VectorString3D for each possible path, as determined by the Project Tree structure and the parent/child relationships.
540 * A pipe is considered a potential branch when it is contained in an abstract child at the same tree level that the pipe itself in the tree. So:
541 * - lineage
542 * - branch 1
543 * -pipe 1
544 * - branch 2
545 * - pipe 2
546 * - branch 3
547 * - pipe 3
548 *
549 * Results in an ArrayList with:
550 * VS 1
551 * VS 1 + 2
552 * VS 1 + 3
553 *
554 * An so on, recursively from the given pipe's parent.
555 */
556 static public ArrayList<Chain> createPipeChains(final ProjectThing root_pt, final LayerSet ls) throws Exception {
560 }
562 /** Recursive. */
563 static private void appendAndFork(final ProjectThing parent, Chain chain, HashSet<ProjectThing> hs_c_done, final ArrayList<Chain> chains, final LayerSet ls) throws Exception {
576 if (!pipe.getLayerSet().equals(ls) || pipe.length() < 2) continue; // not from the same LayerSet, maybe from a nested one.
582 }
583 // find other children in the parent who contain children with child pipes
589 // c is at the same tree level as child (which contains a pipe directly)
592 Chain ca;
594 // just append
598 // otherwise make a copy to branch out
599 ca = base.duplicate(); // can't duplicate from chain itself, because it would have the previous child added to it.
601 }
603 }
604 }
605 // pipe wrapping ProjectThing objects cannot have any children
607 }
609 // if it does not have direct pipe children, cut chain - but keep inspecting
612 }
614 // inspect others down the unvisited tree nodes
616 }
617 }
619 /** Represents a list of concatenated pipes, where each pipe is parent of the next within the ProjectTree. */
627 }
629 //if (pipes.contains(p)) throw new Exception("Already contains pipe #" + p.getId());
632 }
638 }
644 }
650 }
655 //else, chain the ids of the rest
659 }
660 /** Returns max 10 chars, solely the name of the parent's parent node of the root pipe (aka the [lineage] containing the [branch]) or the id if too long. Intended for the 10-digit limitation in the problem in .dis files for Phylip. */
665 // investigate the [branch] title
669 // investigate the lineage title
674 }
675 // check length
678 }
679 // else fall back to unique id
683 }
686 }
688 }
689 /** Returns the color of the root pipe. */
692 }
695 }
696 /** Show centered, set visible and select. */
706 }
709 }
710 }
713 /** Contains all que query chains created from the single pipe selected for matching against another reference project, their matches with the chains made from the reference project, and some general data such as the transforms and the axes. */
718 Calibration cal1;
719 Calibration cal2;
720 // The transfer transform: from query space to reference space
721 Transform3D T;
723 final Hashtable<Chain,ArrayList<ChainMatch>> matches = new Hashtable<Chain,ArrayList<ChainMatch>>();
730 // these chains are kept only calibrated, NOT transformed. Because each query will resample it to its own delta, and then transform it.
737 }
739 /** Will calibrate and transform the chain's VectorString3D. */
742 // Order is important:
743 // 1 - calibrate: bring to user-space microns, whatever
745 // 2 - transform into reference
747 // 3 - resample, within reference space
750 // Store all
752 }
759 }
761 }
763 /** Will calibrate them all to cal2. */
769 }
770 }
772 /** Returns a resampled and transformed copy of the pipe's VectorString3D. */
775 }
776 /** Returns a resampled and transformed copy of the pipe's VectorString3D. */
779 }
785 }
787 /** Bring Chain q (query) to reference space. */
790 }
791 /** Bring Chain q (query) to reference space. Will resample after transforming it, to its own average delta. */
794 }
796 /** Returns a resampled and transformed copy of the pipe's VectorString3D. */
799 }
803 //Utils.log2("VS1: Length before transforming: " + vs.computeLength());
805 //Utils.log2("VS1: Length after transforming: " + vs.computeLength());
806 // Resample after transforming, of course!
809 }
811 /** For each list of matches corresponding to each query chain, sort the matches by physical distance. */
815 }
816 }
817 /** Returns all pipes involved in the query chains. */
822 }
824 }
835 }
836 }
838 }
840 // add all matches
844 }
845 }
846 }
847 // One table entry per query chain
853 QueryHolderTableModel qt = new QueryHolderTableModel(this, new Visualizer(this, vs_axes, vs_axes_ref), query, matches.get(query));
861 }
862 }
865 // from the queries (and thus the tabs as well)
873 }
876 }
877 }
878 // from the matches list of each query chain
883 }
884 }
885 }
887 }
889 /** Remove all queries and refs that belong to the given project. */
891 // from the queries (and thus the tabs as well)
899 }
902 }
903 }
904 // from the matches list of each query chain
909 }
910 }
911 }
913 }
914 }
920 }
922 }
924 /** Represents the scored match between any two Chain objects. */
926 Chain query;
927 Chain ref;
928 Editions ed;
929 double score; // combined score, made from several of the parameters below (S, L and M as of 20080823)
930 double seq_sim; // similarity measure made of num 1 - ((num insertions + num deletions) / max (len1, len2)).
935 double prop_mut; // the proportion of mutation pairs relative to the length of the queried sequence
937 double proximity; // unitless value: cummulative distance of pairs relative to query sequence length
938 double proximity_mut; // unitless value: cummulative distance of only mutation pairs relative to query sequence length
939 ChainMatch(final Chain query, final Chain ref, final Editions ed, final double[] stats, final float prop_len, final double score) {
953 }
954 }
957 Chain query;
962 }
963 }
966 /** Sort by the given distance type. */
970 }
974 // select for smallest physical distance of the center of mass
975 // double val = cm1.phys_dist - cm2.phys_dist;
976 /*
977 final double val = cm1.median - cm2.median;
978 if (val < 0) return -1; // m1 is closer
979 if (val > 0) return 1; // m1 is further away
980 return 0; // same distance
981 */
983 // Select the largest score
984 //
995 val = cm1.seq_sim - cm2.seq_sim; // INVERTED ORDER because each would need an inversion: 1 - cmX.seq_sim
1009 case PROXIMITY: // cummulative distance relative to largest physical length of the two sequences
1012 case PROXIMITY_MUT: // cummulative distance of mutation pairs relative to largest physical length of the two sequences
1015 }
1020 }
1021 }
1026 // select for largest score
1031 }
1032 }
1035 QueryHolder qh;
1036 LayerSet common; // calibrated to queried pipe space, which is now also the space of all others.
1044 // create common LayerSet space
1048 this.common = new LayerSet(pipe.getProject(), pipe.getProject().getLoader().getNextId(), "Common", 10, 10, 0, 0, 0, ls.getLayerWidth() * cal1.pixelWidth, ls.getLayerHeight() * cal1.pixelHeight, false, 2, new AffineTransform());
1052 }
1053 /** Shows the matched chain in 3D. */
1058 // The LayerSet is that of the Line3D being queried, not the given pipe to show which belongs to the reference project (and thus not to the queried pipe project)
1062 }
1068 }
1070 }
1081 }
1082 }
1083 }
1086 ProjectThing pt = (ProjectThing)root.getProject().findProjectThing((ZDisplayable)root).getParent();
1092 VectorString3D vs;
1097 }
1099 }
1100 }
1108 // if already there, ignore request
1114 }
1116 }
1124 }
1129 }
1131 }
1136 // if already there, ignore request
1140 }
1143 }
1144 }
1146 static protected final Object[] findBestMatch(final VectorString3D vs1, final VectorString3D vs2, double delta, boolean skip_ends, int max_mut, float min_chunk) {
1148 }
1150 /** Since comparing two sequences starting from one end or starting from the other
1151 * is not the same at all, this method performs the match starting first from one
1152 * end and then from the other.
1153 * Then it performs a match starting from the middle of the longest stretches of
1154 * pure mutations detected in the best of the matches above.
1155 * Also, since strings may be reversed, the test against the reversed one is done as well.
1156 *
1157 * vs1 vs2
1158 * vs1.reversed() vs2.reversed()
1159 *
1160 * vs1 vs2.reversed()
1161 * vs1.reversed() vs2
1162 *
1163 * ASSUMES both VectorString3D are open.
1164 *
1165 * @param direct Whether to test vs1 against vs2 only, or to try all 4 possible combinations of reversed versus non-reversed and pick the best.
1166 * */
1167 static protected final Object[] findBestMatch(final VectorString3D vs1, final VectorString3D vs2, double delta, boolean skip_ends, int max_mut, float min_chunk, final int distance_type, final boolean direct, final boolean substring_matching) {
1170 // identify shorter chain
1174 // iterate matching of shorter string inside longer string:
1175 // (so that the match is always between two equally long strings)
1176 // aaaaaaaa : 8 elements
1177 // bbbbb : 5 elements
1178 // bbbbb --- total 4 matches to try
1179 // bbbbb
1180 // bbbbb
1181 //
1186 //Utils.log2("substring_matching lengths: shorter, longer : " + shorter.length() + ", " + longer_sub.length());
1192 }
1199 }
1200 }
1201 }
1214 static private final String[] distance_types = {"Levenshtein", "Dissimilarity", "Average physical distance", "Median physical distance", "Cummulative physical distance", "Standard deviation", "Combined SLM", "Proximity", "Proximity of mutation pairs"};
1216 // Weights as empirically approximated with some lineages, with S. Preibisch ( see Test_Scoring.java )
1217 static public final double[] W = new double[]{1.3345290383577453, -0.0012626693452889859, -0.012764729437173508, -0.13344076489951817};
1219 static public final double score(final double seq_sim, final double levenshtein, final double median_phys_dist, final double[] w) {
1220 // S L M
1222 }
1224 /** Zero is best; gets bad towards positive infinite. */
1225 static private final double getScore(Editions ed, boolean skip_ends, int max_mut, float min_chunk, int distance_type) {
1242 return 1 / score(ed.getSimilarity(), ed.getDistance(), ed.getStatistics(skip_ends, max_mut, min_chunk, false)[3], Compare.W);
1243 case PROXIMITY: // cummulative distance relative to largest physical length of the two sequences
1245 case PROXIMITY_MUT: // cummulative distance of mutation pairs relative to largest physical length of the two sequences
1247 }
1249 }
1251 static private final Object[] matchDirect(final VectorString3D vs1, final VectorString3D vs2, double delta, boolean skip_ends, int max_mut, float min_chunk, int distance_type) {
1252 // Levenshtein is unfortunately not commutative: must try both
1260 }
1262 // Match in all possible ways
1263 static private final Object[] matchFwdRev(final VectorString3D vs1, final VectorString3D vs2, double delta, boolean skip_ends, int max_mut, float min_chunk, int distance_type) {
1270 // vs1 vs2
1272 // vs1rev vs2rev
1274 // vs1 vs2rev
1276 // vs1rev vs2
1279 //double best_score1 = 0;
1288 //best_score1 = score1;
1289 }
1290 }
1291 //Utils.log2("score, score1: " + best_score + ", " + best_score1);
1293 // now test also starting from the middle of the longest mutation chunk of the best matching
1296 // is null if no chunks were found
1302 }
1303 }
1306 }
1309 }
1313 // else, record problem
1314 Utils.log("WARNING: the calibration units of the queried pipe (" + unit1 + ") does not match with that of the reference project '" + project_title + "' (" + unit2 + ").");
1316 }
1318 /** Compare the given pipe with other pipes in the given standard project(s). WARNING: the calibrations will work ONLY IF all pipes found to compare with come from LayerSets which have the same units of measurement! For example, all in microns. */
1319 static public final Bureaucrat findSimilar(final Line3D pipe, final Project[] ref, final boolean ignore_orientation, final boolean ignore_calibration, final boolean mirror, final boolean chain_branches, final boolean show_gui) {
1328 }
1333 if (!ignore_calibration) cal = (null != pipe.getLayerSet() ? pipe.getLayerSet().getCalibration() : null);
1341 // collect the full set of ref pipes from each project
1346 }
1347 }
1348 // add all possible query chains, starting at the parent of the chosen pipe
1349 for (Chain chain : createPipeChains((ProjectThing)pipe.getProject().findProjectThing(pipe).getParent(), pipe.getLayerSet())) {
1354 }
1356 }
1358 // no branching: single query of one single-pipe chain
1366 }
1367 }
1368 }
1373 // each thread handles a ref pipe, which is to be matched against all queries and reversed queries
1385 ////
1389 // obtain a ref chain, to be uniquely processed by this thread
1398 //
1403 }
1405 Editions ed;
1412 // the higher the better
1417 }
1422 }
1423 //qh.addMatch(query, ref, ed, score, ed.getPhysicalDistance(false, 0, 1));
1426 qm[q].cm[k] = new ChainMatch(query, ref, ed, stats, prop_len, score(ed.getSimilarity(), ed.getDistance(), ed.getStatistics(false, 0, 0, false)[3], Compare.W));
1427 }
1431 }
1432 }
1433 ////
1434 }
1435 });
1436 }
1439 // done!
1440 // put result into the Worker
1444 // Now, sort matches by physical distance.
1447 qh.sortMatches(ignore_orientation ? new ChainMatchComparatorSim() : new ChainMatchComparator(MEDIAN_PHYS_DIST));
1449 }
1456 }
1457 }
1458 };
1464 }
1467 Displayable displ;
1468 Editions ed;
1475 }
1479 }
1480 }
1486 // select for largest score
1488 //double val = m1.ed.getDistance() - m2.ed.getDistance();
1492 }
1493 }
1499 // select for smallest physical distance of the center of mass
1504 }
1505 }
1508 /* // WOULD HAVE to create my own sorter: can't sort numerically even (it's bit sort), plus the getValueAt(int row, int col) gets messed up - no proper backend data update. But to create my own sorter, I need a TableSorter class which is new 1.6.0, i..e would have tobe very convolutedly generated.
1509 if (ij.IJ.isJava16()) {
1510 try {
1511 java.lang.reflect.Method msort = JTable.class.getMethod("setAutoCreateRowSorter", new Class[]{Boolean.TYPE});
1512 msort.invoke(table, new Object[]{Boolean.TRUE});
1513 } catch (Exception e) {
1514 e.printStackTrace();
1515 }
1516 }
1517 */
1528 }
1536 }
1537 }
1545 }
1546 });
1550 // a listener to change the label text when the tab is selected
1553 if (null == frame || null == ht_tabs || null == tabs || null == label) return; // the event fires during instantiation ... pffff!
1559 }
1560 };
1564 }
1565 };
1573 }
1574 });
1586 tabs.addChangeListener(tabs_listener); // to avoid firing it during instantiation, must be added last
1588 }
1591 }
1600 }
1601 }
1613 }
1615 }
1619 QueryHolder qh;
1620 Chain query;
1622 Visualizer vis;
1629 }
1647 }
1648 }
1656 case 3: return Utils.cutNumber(cm.get(row).seq_sim * 100, 2) + " %"; // 1 - seq_sim, to convert from dissimilarity to similarity
1662 case 9: return Utils.cutNumber(cm.get(row).prop_mut * 100, 2) + " %"; // the proportion of mutations, relative to the query sequence length
1663 case 10: return Utils.cutNumber(cm.get(row).prop_len * 100, 2) + " %"; // the proportion of lengths = len(query) / len(ref)
1667 }
1668 }
1671 }
1674 }
1691 }
1693 }
1696 // select row under mouse
1707 // TODO: need better
1708 // - to show the matched chains alone
1709 // - to show the full matched nodes
1715 // for now, use the first selected only
1720 }
1732 }
1733 }
1734 };
1735 JMenuItem item;
1746 }
1747 }
1748 }
1750 /** Compare all to all parameters. */
1760 public final String[] preset_names = new String[]{"X - 'medial lobe', Y - 'dorsal lobe', Z - 'peduncle'"};
1776 public boolean setup(final boolean to_file, final String regex, final boolean plot, final boolean condense) {
1783 Utils.addEnablerListener((Checkbox)gd.getCheckboxes().get(0), new Component[]{(Component)gd.getNumericFields().get(0), (Component)gd.getNumericFields().get(1)}, null);
1798 }
1808 }
1813 }
1815 //////
1828 }
1850 }
1854 }
1857 }
1858 }
1860 /** Gather chains for all projects considering the cp.regex. */
1861 static public final Object[] gatherChains(final Project[] p, final CATAParameters cp) throws Exception {
1862 // gather all chains
1863 final ArrayList[] p_chains = new ArrayList[p.length]; // to keep track of each project's chains
1869 // Search (shallow) for cp.regex matches
1874 }
1876 }
1878 // calibrate
1881 }
1884 // register all, or relative
1886 // '3' means relative
1887 // compute global average delta
1891 }
1892 }
1898 }
1901 // '1', '2' and '3' involve a 3D affine
1902 // no need //VectorString3D[][] vs_axes = new VectorString3D[p.length][];
1905 // 1 - find pipes to work as axes for each project
1910 }
1913 // if axes are -1, forget it: not found
1918 }
1920 // obtain axes and origin
1925 o[0]); // will be null for the first, which will then be non-null and act as the reference for the others.
1927 // no need //vs_axes[i] = (VectorString3D[])pack[0];
1929 }
1930 /* // OLD WAY
1931 // match the scales to make the largest be 1.0
1932 final double scaling_factor = VectorString3D.matchOrigins(o, transform_type);
1933 Utils.log2("matchOrigins scaling factor: " + scaling_factor + " for transform_type " + transform_type);
1934 */
1935 // transform all except the first (which acts as reference)
1940 // The transfer T transform: from query space to reference space.
1945 }
1946 }
1947 }
1948 // else, direct
1950 // compute global average delta, after correcting calibration and transformation
1954 }
1955 }
1958 // After calibration and transformation, resample all to the same delta
1960 }
1963 }
1965 /** Gets pipes for all open projects, and generates a matrix of dissimilarities, which gets passed on to the Worker thread and also to a file, if desired.
1966 *
1967 * @param to_file Whether to save the results to a file and popup a save dialog for it or not. In any case the results are stored in the worker's load, which you can retrieve like:
1968 * Bureaucrat bu = Compare.compareAllToAll(true);
1969 * Object result = bu.getWorker().getResult();
1970 * float[][] scores = (float[][])result[0];
1971 * ArrayList<Compare.Chain> chains = (ArrayList<Compare.Chain>)result[1];
1972 * @param normalize Whether to normalize the score values so that the maximum value is 1 and the minimum is 0, or not.
1973 */
1976 // gather all open projects
1988 }
2000 }
2003 }
2012 }
2016 // compare all to all
2019 final float[][] scores = Compare.scoreAllToAll(vs, cp.distance_type, cp.delta, cp.skip_ends, cp.max_mut, cp.min_chunk, cp.direct, cp.substring_matching, this);
2024 }
2026 // store matrix and chains into the worker
2029 // write to file
2033 }
2036 final OutputStreamWriter dos = new OutputStreamWriter(new BufferedOutputStream(new FileOutputStream(f)), "8859_1"); // encoding in Latin 1 (for macosx not to mess around
2038 // Normalize matrix to largest value of 1.0
2044 }
2045 }
2049 }
2050 }
2051 }
2053 // write chain titles, with project prefix
2055 // as csv:
2064 titles[next] = new StringBuffer().append('\"').append(prefix).append(' ').append(chain.getCellTitle()).append('\"');
2066 next++;
2067 }
2068 }
2076 }
2080 }
2082 // as XML:
2084 StringBuffer sb = new StringBuffer("<?xml version=\"1.0\"?>\n<!DOCTYPE ggobidata SYSTEM \"ggobi.dtd\">\n");
2089 sb.append("<variables count=\"0\">\n</variables>\n"); // ggobi: what a crappy XML parser it has
2090 sb.append("<records count=\"").append(chains.size()).append("\" glyph=\"fr 1\" color=\"3\">\n");
2096 sb.append("<record id=\"").append(next+1).append("\" label=\"").append(prefix).append(' ').append(chain.getCellTitle()).append("\" ").append(color).append("></record>\n");
2097 next++;
2098 }
2099 }
2105 sb.append("<records count=\"").append(n_chains*(n_chains-1)).append("\" glyph=\"fr 1\" color=\"0\">\n");
2109 sb.append("<record source=\"").append(i+1).append("\" destination=\"").append(j+1).append("\">").append(scores[i][j]).append("</record>\n");
2110 }
2111 }
2121 }
2123 // as Phylip .dis
2125 // collect different projects
2130 }
2137 String title = chains.get(i).getShortCellTitle().replace(' ', '_').replace('\t', '_').replace('[', '-').replace(']', '-');
2138 // Crop title to 6 chars
2143 }
2146 // Prepend a project char identifier
2149 project_name = Utils.getCharacter(projects.indexOf(chains.get(i).getRoot().getProject()) + 1).toLowerCase();
2151 }
2152 // Append a char index when name is used multiple times (mostly because of branching, and other reasons)
2160 k++;
2162 }
2164 //
2171 count++;
2178 }
2179 }
2182 }
2186 }
2187 }
2196 }
2197 }
2198 };
2200 }
2202 /** Returns the half matrix of scores, with values copied from one half matrix to the other, and a diagonal of zeros.
2203 * @param distance_type ranges from 0 to 5, and includes: 0=Levenshtein, 1=Dissimilarity, 2=Average physical distance, 3=Median physical distance, 4=Cummulative physical distance and 5=Standard deviation. */
2204 static public float[][] scoreAllToAll(final VectorString3D[] vs, final int distance_type, final double delta, final boolean skip_ends, final int max_mut, final float min_chunk, final boolean direct, final boolean substring_matching, final Worker worker) {
2213 ////
2219 final Object[] ob = findBestMatch(vs[i], vs[j], delta, skip_ends, max_mut, min_chunk, distance_type, direct, substring_matching); // TODO should add 'distance_type' as well for the selection of the best match when not direct.
2220 /*
2221 switch (distance_type) {
2222 case 0: // Levenshtein
2223 scores[i][j] = (float)((Editions)ob[0]).getDistance();
2224 break;
2225 case 1: // dissimilarity
2226 scores[i][j] = (float)((Double)ob[1]).doubleValue();
2227 break;
2228 case 2: // average physical distance between mutation pairs
2229 scores[i][j] = (float)((Editions)ob[0]).getPhysicalDistance(skip_ends, max_mut, min_chunk, true);
2230 break;
2231 case 3: // median physical distance between mutation pairs
2232 scores[i][j] = (float)((Editions)ob[0]).getStatistics(skip_ends, max_mut, min_chunk, false)[3]; // 3 is median
2233 break;
2234 case 4: // cummulative physical distance between mutation pairs
2235 scores[i][j] = (float)((Editions)ob[0]).getPhysicalDistance(skip_ends, max_mut, min_chunk, false);
2236 break;
2237 case 5: // stdDev of distances between mutation pairs
2238 scores[i][j] = (float)((Editions)ob[0]).getStdDev(skip_ends, max_mut, min_chunk);
2239 break;
2240 }
2241 */
2247 // mirror value
2249 }
2250 }
2252 ////
2253 }};
2254 }
2260 }
2262 /** Creates a transform with the 4 given vectors: X, Y, Z and translation of origin. */
2265 // X Y Z Trans
2270 }
2274 }
2276 /** If @param reference_project is null, then the first one found in the Project.getProjects() lists is used.
2277 * If regex is not null, then only ProjectThing nodes with the matching regex are analyzed (shallow: none of their children are questioned, but pipes will be built from them all).
2278 * if show_plots, plots are shown -- otherwise, a directory is asked for and they are saved as png files.
2279 * */
2280 static public Bureaucrat variabilityAnalysis(final Project reference_project, final String regex, final boolean show_plots, final boolean show_3D, final Map<String,VectorString3D> map_condensed) {
2281 // gather all open projects
2283 // make the reference_project be the first in the array
2290 }
2291 }
2292 }
2305 }
2306 cp.with_source = true; // so source points are stored in VectorString3D for each resampled and interpolated point
2310 // Save plots
2316 }
2319 }
2323 Object[] ob = gatherChains(p, cp); // will transform them as well to the reference found in the first project in the p array
2330 }
2334 // Sort out into groups by unique names of lineage bundles
2339 // ignore:
2347 // DEBUG:
2348 //if (! (title.startsWith("DPLd") || title.startsWith("BAmv1")) ) continue;
2358 }
2360 }
2370 // Condense each into a single VectorString3D
2376 }
2381 }
2385 // Gather source for each, compute stdDev at each point and make a plot with it
2386 // X axis: from first to last point
2387 // Y axis: the stdDev at each point, computed from the group of points that contribute to each
2392 //FAILS//plot.addLabel(10, cp.plot_height-5, name); // must be added after setting size
2394 else if (null != plot_dir) new FileSaver(plot.getImagePlus()).saveAsPng(plot_dir + name.replace('/', '-') + ".png");
2395 }
2397 // Show all condensed in 3D
2399 final LayerSet common_ls = new LayerSet(p[0], -1, "Common", 10, 10, 0, 0, 0, 512, 512, false, 2, new AffineTransform());
2406 }
2407 }
2410 }
2418 }
2419 }
2420 };
2422 }
2424 static public Plot makePlot(final CATAParameters cp, final String name, final VectorString3D c) {
2428 final Plot plot = new Plot(name, name + " -- Point index (delta: " + Utils.cutNumber(c.getDelta(), 2) + " " + c.getCalibrationCopy().getUnits() + ")", "Std Dev", index, stdDev);
2433 }
2435 /** From a condensed VectorString3D, create the radius at each point. */
2438 // one std dev
2440 }
2445 // two or three std dev
2450 }
2452 // average distance from condensed to all sources
2458 i++;
2459 }
2461 // max distance from condensed to all sources
2467 i++;
2468 }
2469 }
2472 }
2479 }
2485 }
2486 }
2488 /** Do an all-to-all distance matrix of the given vs, then do a neighbor joining, do a weighted merge of the two VectorString3D being merged, and then finally output the resulting condensed unique VectorString3D with its source array full with all points that make each point in it. Expects VectorString3D which are already calibrated and transformed. */
2489 static public VectorString3D condense(final CATAParameters cp, final VectorString3D[] vs, final Worker worker) throws Exception {
2490 // Trivial case 1:
2493 // Estimate delta
2497 }
2499 }
2500 // Resample all:
2503 // Trivial case 2:
2505 if (2 == vs.length) VectorString3D.createInterpolatedPoints(new Editions(vs[0], vs[1], cp.delta, false), 0.5f);
2509 }
2511 // Else, do neighbor joining
2512 final float[][] scores = Compare.scoreAllToAll(vs, cp.distance_type, cp.delta, cp.skip_ends, cp.max_mut, cp.min_chunk, cp.direct, cp.substring_matching, worker);
2513 final HashMap<Compare.Cell<VectorString3D>,Float> table = new HashMap<Compare.Cell<VectorString3D>,Float>();
2514 // Input the half matrix only into the table, since it's mirrored. And without the diagonal of zeros:
2518 }
2519 }
2527 }
2528 // find smallest value
2536 }
2537 }
2538 // pop cells from table
2539 // done below//table.remove(cell);
2545 }
2546 }
2547 // pop the two merged VectorString3D
2551 // merge, weighted by number of sources of each
2552 final double alpha = (double)(cell.t1.getNSources()) / (double)(cell.t1.getNSources() + cell.t2.getNSources()); // in createInterpolated, the alpha is the opposite of what one would think: a 0.2 alpha means 0.8 for the first and 0.2 for the second. So alpha should be 1-alpha
2556 // crop ends to eliminate strings of insertions or deletions sparsed by strings of max cp.max_mut mutations inside
2557 // (This reduces or eliminates variability noise caused by unequal sequence length)
2564 n_mut++;
2568 }
2569 }
2570 }
2574 n_mut++;
2578 }
2579 }
2580 }
2584 }
2587 // add a new cell for each possible comparison with all other unique vs
2589 Object[] ob = findBestMatch(vs_merged, v, cp.delta, cp.skip_ends, cp.max_mut, cp.min_chunk, cp.distance_type, cp.direct, cp.substring_matching);
2593 }
2595 // add the new VectorString3D
2597 }
2599 // Only the last vs_merged should remain:
2601 // test:
2604 }
2607 }
2609 /** Iterates all tabs, removing those that contain queries from the given project. */
2612 // avoid concurrent modifications
2620 // uff!
2626 }
2627 }
2628 }
2630 /** Transform all points of all VectorString3D in vs using a Moving Least Squares Transform defined by the pairing of points in source to those in target.
2631 * In short, bring source into target. */
2632 static public List<VectorString3D> transferVectorStrings(final List<VectorString3D> vs, final List<Tuple3d> source, final List<Tuple3d> target) throws Exception {
2634 Utils.log2("Could not generate a MovingLeastSquaresTransform: different number of source and target points.");
2636 }
2640 }
2642 // 1 - Create the MovingLeastSquaresTransform from the point matches
2648 pm.add(new PointMatch(new mpicbg.models.Point(new float[]{(float)sp.x, (float)sp.y, (float)sp.z}), new mpicbg.models.Point(new float[]{(float)tp.x, (float)tp.y, (float)tp.z}), 1));
2649 }
2657 // 1.1 - Test: transfer source points
2670 }
2672 // 2 - Transfer each VectorString3D in vs with mls
2675 // The points of the VectorString3D:
2679 // Empty arrays to fill with the points to transfer:
2683 // Transfer point by point:
2692 }
2696 }
2699 }
2701 /** Transfer vs via a moving least squares transform by matching source named points into equally named target named points.
2702 * If no points in common, returns null. */
2703 static public List<VectorString3D> transferVectorStrings(final List<VectorString3D> vs, final Map<String,Tuple3d> source, final Map<String,Tuple3d> target) throws Exception {
2712 }
2713 }
2717 }
2720 }
2722 static public List<VectorString3D> transferVectorStrings(final List<VectorString3D> vs, final ProjectThing source_fiduciary, final ProjectThing target_fiduciary) throws Exception {
2723 return transferVectorStrings(vs, extractPoints(source_fiduciary), extractPoints(target_fiduciary));
2724 }
2726 /** Extracts the list of fiduciary points from the fiducial parent and, if their name is different than "ball", adds their title as key and their first ball as a fiduciary point value of the returned map. The map is never null but could be empty.
2727 * The values are calibrated. */
2732 }
2737 }
2744 }
2745 // pick the first one only
2750 // get the first one only
2753 }
2754 }
2756 }
2757 }