-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathCellMigration.ipf
2458 lines (2241 loc) · 81.6 KB
/
CellMigration.ipf
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3 // Use modern global access method and strict wave access
#pragma DefaultTab={3,20,4} // Set default tab width in Igor Pro 9 and later
#pragma version = 1.15 // version number of Migrate()
#pragma ModuleName = CellMigration
#include <Waves Average>
#include <ColorWaveEditor>
// CellMigration will analyse 2D cell migration in IgorPro
// Use ImageJ to track the cells. Outputs from tracking are saved either
// 1) as sheets in an Excel Workbook, 1 per condition, or
// 2) as direct outputs from Manual Tracking, in csv format
//
// Select CellMigr > Cell Migration...
//
// Tell the dialog how many conditions you want to load and the magnification and time resolution
// Next, give each condition a label (name) and then tell Igor where to find the data
// Either an Excel workbook (per condition) or a directory of CSVs (per condition)
// There is also the ability to define "offsetting data", in the case of XY drift during the experiment.
//
// For Excel workbooks:
// NOTE no headers in Excel file. Keep data to columns A-H, max of 2000 rows
// columns are
// A - 0 - ImageJ row
// B - 1 - Track No
// C - 2 - Slice No
// D - 3 - x (in px)
// E - 4 - y (in px)
// F - 5 - distance
// G - 6 - velocity (actually speed)
// H - 7 - pixel value
////////////////////////////////////////////////////////////////////////
// Menu items
////////////////////////////////////////////////////////////////////////
Menu "CellMigr"
"Cell Migration...", /Q, SetUpMigration(0)
"Superplot...", /Q, SetUpMigration(1)
"Save Reports...", /Q, SaveAllReports()
"Recolor Everything", /Q, RecolorAllPlotsWrapper()
"Rerun Analysis", /Q, RerunAnalysis()
Submenu "Manual tracking conversion"
"Excel to Converted CSV", /Q, Excel2CSV()
"CSV to Converted CSV", /Q, CSV2CSV()
End
"About CellMigration", /Q, AboutCellMigr()
End
////////////////////////////////////////////////////////////////////////
// Master functions and wrappers
////////////////////////////////////////////////////////////////////////
Function SetUpMigration(optVar)
Variable optVar
SetDataFolder root:
// kill all windows and waves before we start
CleanSlate()
SetUp_Panel(optVar)
End
Function ProceedToMigrate()
SetDataFolder root:
WAVE/Z paramWave
MakeColorWave(paramWave[0],"colorWave")
if(paramWave[4] > 0)
Superplot_Panel(paramWave[0],paramWave[4])
else
myIO_Panel(paramWave[0])
endif
End
////////////////////////////////////////////////////////////////////////
// Main functions
////////////////////////////////////////////////////////////////////////
// Loads the data and performs migration analysis
Function Migrate()
WAVE/Z paramWave = root:paramWave
if(!WaveExists(paramWave))
DoAlert 0, "Setup has failed. Missing paramWave."
return -1
endif
// pick up global values needed
Variable cond = paramWave[0]
Variable tStep = paramWave[1]
Variable pxSize = paramWave[2]
Variable superPlot = 0, reps
WAVE/Z colorWave = root:colorWave
WAVE/T/Z condWave = root:condWave
WAVE/T/Z condSplitWave = root:condSplitWave
if(WaveExists(condSplitWave) == 1)
if(cond < numpnts(condSplitWave))
superPlot = 1
cond = numpnts(condSplitWave)
reps = cond / paramWave[0]
endif
endif
// because the user may have used illegal characters in condWave, we make a clean version
// for use in Igor and a copy of the original called labelWave to use in plots and layouts
if(superPlot == 1)
WAVE/T labelWave = CleanUpCondWave(condSplitWave)
WAVE/T labelWave = CleanUpCondWave(condWave)
else
WAVE/T labelWave = CleanUpCondWave(condWave)
endif
// make summary plot windows
String fullList = "cdPlot;ivPlot;ivHPlot;dDPlot;MSDPlot;DAPlot;angleHPlot"
Variable nPlots = ItemsInList(fullList)
String name
Variable i,j
for(i = 0; i < nPlots; i += 1)
name = StringFromList(i, fullList)
KillWindow/Z $name
Display/N=$name/HIDE=1
endfor
String dataFolderName = "root:data"
NewDataFolder/O $dataFolderName // make root:data: but don't put anything in it yet
String condName, pref
Variable moviemax1, moviemax2
for(i = 0; i < cond; i += 1)
if(superPlot == 1)
condName = condSplitWave[i]
else
condName = condWave[i]
endif
// make data folder for each condition
dataFolderName = "root:data:" + condName
NewDataFolder/O/S $dataFolderName
// run other procedures
moviemax1 = LoadMigration(i)
moviemax2 = CorrectMigration(i)
if(moviemax1 != moviemax2)
if(moviemax2 == -1)
print "No correction applied to", condName
else
print "Caution: different number of stationary tracks compared with real tracks."
endif
endif
if(superPlot == 0)
// for each condition go and make tracks and plot everything out
MakeTracks(i)
endif
SetDataFolder root:
endfor
// if we are making a superplot, copy the data into folders named by condWave
if(superPlot == 1)
String sourceDF = ""
cond = numpnts(condWave)
for(i = 0; i < cond; i += 1)
condName = condWave[i]
// make data folder for each master condition
dataFolderName = "root:data:" + condName
NewDataFolder/O $dataFolderName
for(j = 0; j < reps; j += 1)
sourceDF = "root:data:" + condSplitWave[i * reps + j]
// merge data into master condition folder
DuplicateDataFolder/O=2 $sourceDF, $dataFolderName
endfor
SetDataFolder $dataFolderName
// as in the main program: for each condition go and make tracks and plot everything out
MakeTracks(i)
SetDataFolder root:
endfor
endif
// make the image quilt, spakline and joint histogram and then sort out the layouts
Variable optDur = MakeImageQuilt(10) // this aims for a quilt of 100 = 10^2 tracks
MakeJointHistogram(optDur)
TidyCondSpecificLayouts()
KillWindow/Z summaryLayout
NewLayout/N=summaryLayout
TidyUpSummaryLayout()
// when we get to the end, print (pragma) version number
Print "*** Executed Migrate v", GetProcedureVersion("CellMigration.ipf")
KillWindow/Z FilePicker
End
// This function will load the tracking data
/// @param ii variable containing row number from condWave
Function LoadMigration(ii)
Variable ii
WAVE/T/Z condSplitWave = root:condSplitWave
WAVE/T/Z condWave = root:condWave
String condName
if(WaveExists(condSplitWave) == 1)
condName = condSplitWave[ii]
else
condName = condWave[ii]
endif
WAVE/T PathWave1 = root:PathWave1
String pathString = PathWave1[ii]
String sheet, prefix, matName, wList
String fileList
Variable moviemax,csvOrNot
Variable i
if(StringMatch(pathString, "*.xls*") == 1)
// set variable to indicate Excel Workbook
csvOrNot = 0
// Works out what sheets are in Excel Workbook and then loads each.
XLLoadWave/J=1 PathWave1[ii]
fileList = S_value
else
// set variable to indicate csv file
csvOrNot = 1
// Work out what files are in directory
NewPath/O/Q ExpDiskFolder, pathString
fileList = IndexedFile(expDiskFolder,-1,".csv")
endif
fileList = SortList(fileList, ";", 16)
moviemax = ItemsInList(fileList)
for(i = 0; i < moviemax; i += 1)
sheet = StringFromList(i, fileList)
prefix = condName + "_c_" + num2str(i)
matName = condName + "_" + num2str(i)
if(csvOrNot == 0)
XLLoadWave/S=sheet/R=(A1,H2000)/O/K=0/N=$prefix/Q PathWave1[ii]
else
LoadWave/A=$prefix/J/K=1/L={0,1,0,0,0}/O/P=expDiskFolder/Q sheet
endif
wList = wavelist(prefix + "*",";","") // make matrix for each sheet
Concatenate/O/KILL wList, $matName
// now we need to check that the matrix is OK
Wave matTrax = $matName
// check we have a valid matrix with 8 columns
CheckColumnsOfMatrix(matTrax)
// make sure 1st point is -1
matTrax[0][5,6] = -1
// check distances and speeds are correct
CheckDistancesAndSpeeds(matTrax)
endfor
Print "*** Condition", condName, "was loaded from", pathString
// return moviemax back to calling function for checking
return moviemax
End
// This was added in v1.12 Dec 2018. At some point in the last year, the output
// of Manual Tracking changed so that there was no longer a first (0) column containing
// the row numbers. Add this column if that is the case.
STATIC Function CheckColumnsOfMatrix(matTrax)
WAVE matTrax
Variable numCols = DimSize(matTrax,1)
if(numCols == 8)
return 1
elseif(numCols == 7)
// insert new 0th column
InsertPoints/M=1 0,1, MatTrax
MatTrax[][0] = p+1
return 0
else
Print NameOfWave(matTrax), "does not have 8 columns of data."
return -1
endif
End
// The purpose of this function is to work out whether the distances (and speeds) in the
// original data are correct. Currently it just corrects them rather than testing and correcting if needed.
STATIC Function CheckDistancesAndSpeeds(matTrax)
WAVE matTrax
WAVE/Z paramWave = root:paramWave
Variable tStep = paramWave[1]
Variable pxSize = paramWave[2]
// make new distance column
Duplicate/O/RMD=[][3,4]/FREE matTrax,tempDist // take offset coords
Differentiate/METH=2/DIM=0 tempDist
tempDist[][] = (matTrax[p][5] == -1) ? 0 : tempDist[p][q]
MatrixOp/O/FREE tempNorm = sqrt(sumRows(tempDist * tempDist))
tempNorm[] *= pxSize // convert to real distance
// MatrixOp/O/FREE tempReal = sumcols(tempNorm - col(matTrax,5)) // for checking
matTrax[][5] = (matTrax[p][5] == -1) ? -1 : tempNorm[p] // going to leave first point as -1
// correct speed column
matTrax[][6] = (matTrax[p][6] == -1) ? -1 : tempNorm[p] / tStep
// make sure 1st point is -1
matTrax[0][5,6] = -1
End
// This function will load the tracking data from an Excel Workbook
/// @param ii variable containing row number from condWave
Function CorrectMigration(ii)
Variable ii
WAVE/T/Z condSplitWave = root:condSplitWave
WAVE/T/Z condWave = root:condWave
String condName
if(WaveExists(condSplitWave) == 1)
condName = condSplitWave[ii]
else
condName = condWave[ii]
endif
WAVE/T PathWave2 = root:PathWave2
String pathString = PathWave2[ii]
Variable len = strlen(pathString)
if(len == 0)
return -1
elseif(numtype(len) == 2)
return -1
endif
String sheet, prefix, matName, wList, mName
String fileList
Variable moviemax,csvOrNot
Variable i
if(StringMatch(pathString, "*.xls*") == 1)
// set variable to indicate Excel Workbook
csvOrNot = 0
// Works out what sheets are in Excel Workbook and then loads each.
XLLoadWave/J=1 PathWave2[ii]
fileList = S_value
else
// set variable to indicate csv file
csvOrNot = 1
// Work out what files are in directory
NewPath/O/Q ExpDiskFolder, pathString
fileList = IndexedFile(expDiskFolder,-1,".csv")
endif
fileList = SortList(fileList, ";", 16)
moviemax = ItemsInList(fileList)
for(i = 0; i < moviemax; i += 1)
sheet = StringFromList(i,fileList)
prefix = "stat_" + "c_" + num2str(i) // use stat prefix
matName = "stat_" + num2str(i)
if(csvOrNot == 0)
XLLoadWave/S=sheet/R=(A1,H2000)/O/K=0/N=$prefix/Q PathWave2[ii]
else
LoadWave/A=$prefix/J/K=1/L={0,1,0,0,0}/O/P=expDiskFolder/Q sheet
endif
wList = wavelist(prefix + "*",";","") // make matrix for each sheet
Concatenate/O/KILL wList, $matName
Wave matStat = $matName
// Find corresponding movie matrix
mName = ReplaceString("stat_",matname,condName + "_")
Wave matTrax = $mName
OffsetAndRecalc(matStat,matTrax)
endfor
Print "*** Offset data for condition", condName, "was loaded from", pathString
// return moviemax back to calling function for checking
return moviemax
End
// This function uses matStat to offset matTrax
Function OffsetAndRecalc(matStat,matTrax)
Wave matStat,matTrax
// Work out offset for the stat_* waves
Variable x0 = matStat[0][3]
Variable y0 = matStat[0][4]
matStat[][3] -= x0
matStat[][4] -= y0
MatrixOp/O/FREE mStat2 = col(matStat,2)
Variable maxFrame = WaveMax(mStat2)
Variable j // because i refers to rows
// offsetting loop
for(j = 1; j < maxFrame + 1; j += 1)
FindValue/V=(j) mStat2
if(V_Value == -1)
x0 = 0
y0 = 0
else
x0 = matStat[V_Value][3]
y0 = matStat[V_Value][4]
endif
matTrax[][3] = (matTrax[p][2] == j) ? matTrax[p][3] - x0 : matTrax[p][3]
matTrax[][4] = (matTrax[p][2] == j) ? matTrax[p][4] - y0 : matTrax[p][4]
endfor
WAVE/Z paramWave = root:paramWave
Variable tStep = paramWave[1]
Variable pxSize = paramWave[2]
// make new distance column
Duplicate/O/RMD=[][3,4]/FREE matTrax,tempDist // take offset coords
Differentiate/METH=2 tempDist
tempDist[][] = (matTrax[p][5] == -1) ? 0 : tempDist[p][q]
MatrixOp/O/FREE tempNorm = sqrt(sumRows(tempDist * tempDist))
tempNorm[] *= pxSize // convert to real distance
matTrax[][5] = (matTrax[p][5] == -1) ? -1 : tempNorm[p] // going to leave first point as -1
// correct speed column
matTrax[][6] = (matTrax[p][6] == -1) ? -1 : tempNorm[p] / tStep
// put 1st point as -1
matTrax[0][5,6] = -1
End
// This function will make cumulative distance waves for each cell. They are called cd_*
/// @param ii variable containing row number from condWave
Function MakeTracks(ii)
Variable ii
WAVE/T condWave = root:condWave
String condName = condWave[ii]
WAVE/Z paramWave = root:paramWave
Variable tStep = paramWave[1]
Variable pxSize = paramWave[2]
Variable minPoints = paramWave[5]
WAVE/Z colorWave = root:colorWave
WAVE/Z/T unitWave = root:unitWave
String wList0 = WaveList(condName + "_*",";","") // find all matrices
Variable nWaves = ItemsInList(wList0)
Variable nTrack,nTrace=0
String mName0, newName, plotName, avList, avName, errName
Variable i, j
String layoutName = condName + "_layout"
KillWindow/Z $layoutName // Kill the layout if it exists
NewLayout/HIDE=1/N=$layoutName
// cumulative distance and plot over time
plotName = condName + "_cdplot"
KillWindow/Z $plotName // set up plot
Display/N=$plotName/HIDE=1
for(i = 0; i < nWaves; i += 1)
mName0 = StringFromList(i,wList0)
WAVE m0 = $mName0
Duplicate/O/RMD=[][5,5] m0, $"tDistW" // distance
Duplicate/O/RMD=[][1,1] m0, $"tCellW" // cell number
WAVE tDistW,tCellW
Redimension/N=-1 tDistW, tCellW // make 1D
nTrack = WaveMax(tCellW) // find maximum track number
for(j = 1; j < (nTrack+1); j += 1) // index is 1-based
newName = "cd_" + mName0 + "_" + num2str(j)
Duplicate/O tDistW, $newName
WAVE w2 = $newName
w2 = (tCellW[p] == j) ? tDistW[p] : NaN
WaveTransform zapnans w2
if(numpnts(w2) < minPoints)
KillWaves/Z w2 // delete short tracks and any tracks that didn't exist
else
w2[0] = 0 // first point in distance trace is -1 so correct this
Integrate/METH=0 w2 // make cumulative distance
SetScale/P x 0,tStep, w2
AppendtoGraph/W=$plotName $newName
nTrace += 1
endif
endfor
KillWaves/Z tDistW
endfor
Variable alphaLevel = DecideOpacity(nTrace)
ModifyGraph/W=$plotName rgb=(colorWave[ii][0],colorWave[ii][1],colorWave[ii][2],alphaLevel)
avList = Wavelist("cd*",";","WIN:"+ plotName)
avName = "W_Ave_cd_" + condName
errName = ReplaceString("Ave", avName, "Err")
fWaveAverage(avList, "", 3, 1, avName, errName)
AppendToGraph/W=$plotName $avName
Label/W=$plotName left "Cumulative distance ("+unitWave[1]+")"
Label/W=$plotName bottom "Time ("+unitWave[0]+")"
ErrorBars/W=$plotName $avName SHADE= {0,0,(0,0,0,0),(0,0,0,0)},wave=($ErrName,$ErrName)
ModifyGraph/W=$plotName lsize($avName)=2,rgb($avName)=(0,0,0)
SetAxis/W=$plotName/A/N=1 left
AppendLayoutObject/W=$layoutName graph $plotName
// instantaneous speed over time
plotName = condName + "_ivplot"
KillWindow/Z $plotName // set up plot
Display/N=$plotName/HIDE=1
for(i = 0; i < nWaves; i += 1)
mName0 = StringFromList(i,wList0)
WAVE m0 = $mName0
Duplicate/O/RMD=[][5,5] m0, $"tDistW" // distance
Duplicate/O/RMD=[][1,1] m0, $"tCellW" // cell number
WAVE tDistW,tCellW
Redimension/N=-1 tDistW, tCellW // make 1D
nTrack = WaveMax(tCellW) // find maximum track number
for(j = 1; j < (nTrack+1); j += 1) // index is 1-based
newName = "iv_" + mName0 + "_" + num2str(j)
Duplicate/O tDistW, $newName
WAVE w2 = $newName
w2 = (tCellW[p] == j) ? tDistW[p] : NaN
WaveTransform zapnans w2
if(numpnts(w2) <= minPoints)
KillWaves w2
else
w2[0] = 0 // first point in distance trace is -1, so correct this
w2 /= tStep // make instantaneous speed
SetScale/P x 0,tStep, w2
AppendtoGraph/W=$plotName $newName
endif
endfor
KillWaves/Z tDistW
endfor
ModifyGraph/W=$plotName rgb=(colorWave[ii][0],colorWave[ii][1],colorWave[ii][2],alphaLevel)
avList = Wavelist("iv*",";","WIN:"+ plotName)
avName = "W_Ave_iv_" + condName
errName = ReplaceString("Ave", avName, "Err")
fWaveAverage(avList, "", 3, 1, avName, errName)
AppendToGraph/W=$plotName $avName
Label/W=$plotName left "Instantaneous Speed ("+unitWave[2]+")"
Label/W=$plotName bottom "Time ("+unitWave[0]+")"
ErrorBars/W=$plotName $avName SHADE= {0,0,(0,0,0,0),(0,0,0,0)},wave=($ErrName,$ErrName)
ModifyGraph/W=$plotName lsize($avName)=2,rgb($avName)=(0,0,0)
SetAxis/W=$plotName/A/N=1 left
AppendLayoutObject/W=$layoutName graph $plotName
// print a message to say how many valid tracks we have in this condition
Print ItemsInList(avList), "valid tracks plotted for", condName
plotName = condName + "_ivHist"
KillWindow/Z $plotName //set up plot
Display/N=$plotName/HIDE=1
Concatenate/O/NP avList, tempwave
newName = "h_iv_" + condName // note that this makes a name like h_iv_Ctrl
// before v 1.15 we used sqrt((3*pxsize)^2)/tStep
print wavemax(tempwave)
Variable nBins = ceil(1 + log(numpnts(tempwave))/log(2))
Variable binWidth = wavemax(tempwave)/nBins
Make/O/N=(nBins) $newName
Histogram/B={0,binWidth,nBins} tempwave,$newName
AppendToGraph/W=$plotName $newName
ModifyGraph/W=$plotName rgb=(colorWave[ii][0],colorWave[ii][1],colorWave[ii][2])
ModifyGraph/W=$plotName mode=5,hbFill=4
SetAxis/W=$plotName/A/N=1/E=1 left
SetAxis/W=$plotName/A/N=1/E=1 bottom
Label/W=$plotName left "Frequency"
Label/W=$plotName bottom "Instantaneous Speed ("+unitWave[2]+")"
KillWaves/Z tempwave
AppendLayoutObject/W=$layoutName graph $plotName
// plot out tracks
plotName = condName + "_tkplot"
KillWindow/Z $plotName //set up plot
Display/N=$plotName/HIDE=1
Variable off
for(i = 0; i < nWaves; i += 1)
mName0 = StringFromList(i,wList0)
WAVE m0 = $mName0
Duplicate/O/RMD=[][3,3] m0, $"tXW" //x pos
Duplicate/O/RMD=[][4,4] m0, $"tYW" //y pos
Duplicate/O/RMD=[][1,1] m0, $"tCellW" //track number
WAVE tXW,tYW,tCellW
Redimension/N=-1 tXW,tYW,tCellW
nTrack = WaveMax(tCellW) //find maximum track number
for(j = 1; j < (nTrack+1); j += 1) //index is 1-based
newName = "tk_" + mName0 + "_" + num2str(j)
Duplicate/O tXW, $"w3"
WAVE w3
w3 = (tCellW[p] == j) ? w3[p] : NaN
WaveTransform zapnans w3
if(numpnts(w3) <= minPoints)
KillWaves w3
else
off = w3[0]
w3 -= off //set to origin
w3 *= pxSize
// do the y wave
Duplicate/O tYW, $"w4"
WAVE w4
w4 = (tCellW[p] == j) ? w4[p] : NaN
WaveTransform zapnans w4
off = w4[0]
w4 -= off
w4 *= pxSize
Concatenate/O/KILL {w3,w4}, $newName
Wave w5 = $newName
AppendtoGraph/W=$plotName w5[][1] vs w5[][0]
endif
endfor
Killwaves/Z tXW,tYW,tCellW //tidy up
endfor
ModifyGraph/W=$plotName rgb=(colorWave[ii][0],colorWave[ii][1],colorWave[ii][2],alphaLevel)
// set these graph limits temporarily
SetAxis/W=$plotName left -250,250
SetAxis/W=$plotName bottom -250,250
ModifyGraph/W=$plotName width={Plan,1,bottom,left}
ModifyGraph/W=$plotName mirror=1
ModifyGraph/W=$plotName grid=1
ModifyGraph/W=$plotName gridRGB=(32767,32767,32767)
AppendLayoutObject/W=$layoutName graph $plotName
// calculate d/D directionality ratio
plotName = condName + "_dDplot"
KillWindow/Z $plotName // setup plot
Display/N=$plotName/HIDE=1
String wName0, wName1
Variable len
wList0 = WaveList("tk_" + condName + "_*", ";","")
nWaves = ItemsInList(wList0)
for(i = 0; i < nWaves; i += 1)
wName0 = StringFromList(i,wList0) // tk wave
wName1 = ReplaceString("tk",wName0,"cd") // cd wave
WAVE w0 = $wName0
WAVE w1 = $wName1
newName = ReplaceString("tk",wName0,"dD")
Duplicate/O w1 $newName
WAVE w2 = $newName
len = numpnts(w2)
w2[] = (w1[p] == 0) ? 1 : sqrt(w0[p][0]^2 + w0[p][1]^2) / w1[p]
w2[0] = NaN // d/D at point 0 is not a number
AppendtoGraph/W=$plotName w2
endfor
ModifyGraph/W=$plotName rgb=(colorWave[ii][0],colorWave[ii][1],colorWave[ii][2],alphaLevel)
avList = Wavelist("dD*",";","WIN:"+ plotName)
avName = "W_Ave_dD_" + condName
errName = ReplaceString("Ave", avName, "Err")
fWaveAverage(avList, "", 3, 1, avName, errName)
AppendToGraph/W=$plotName $avName
Label/W=$plotName left "Directionality ratio (d/D)"
Label/W=$plotName bottom "Time ("+unitWave[0]+")"
ErrorBars/W=$plotName $avName SHADE= {0,0,(0,0,0,0),(0,0,0,0)},wave=($ErrName,$ErrName)
ModifyGraph/W=$plotName lsize($avName)=2,rgb($avName)=(0,0,0)
AppendLayoutObject/W=$layoutName graph $plotName
// calculate MSD (overlapping method)
plotName = condName + "_MSDplot"
KillWindow/Z $plotName //setup plot
Display/N=$plotName/HIDE=1
wList0 = WaveList("tk_" + condName + "_*", ";","")
nWaves = ItemsInList(wList0)
Variable k
for(i = 0; i < nWaves; i += 1)
wName0 = StringFromList(i,wList0) // tk wave
WAVE w0 = $wName0
len = DimSize(w0,0)
newName = ReplaceString("tk",wName0,"MSD") // for results of MSD per cell
Make/O/N=(len-1,len-1,2)/FREE tempMat0,tempMat1
// make 2 3D waves. 0 is end point to measure MSD, 1 is start point
// layers are x and y
tempMat0[][][] = (p >= q) ? w0[p+1][r] : 0
tempMat1[][][] = (p >= q) ? w0[p-q][r] : 0
MatrixOp/O/FREE tempMat2 = (tempMat0 - tempMat1) * (tempMat0 - tempMat1))
Make/O/N=(len-1)/FREE countOfMSDPnts = (len-1)-p
MatrixOp/O $newName = sumcols(sumbeams(tempMat2))^t / countOfMSDPnts
SetScale/P x 0,tStep, $newName
AppendtoGraph/W=$plotName $newName
endfor
ModifyGraph/W=$plotName rgb=(colorWave[ii][0],colorWave[ii][1],colorWave[ii][2],alphaLevel)
avList = Wavelist("MSD*",";","WIN:"+ plotName)
avName = "W_Ave_MSD_" + condName
errName = ReplaceString("Ave", avName, "Err")
fWaveAverage(avList, "", 3, 1, avName, errName)
AppendToGraph/W=$plotName $avName
ModifyGraph/W=$plotName log=1
SetAxis/W=$plotName/A/N=1 left
len = numpnts($avName)*tStep
SetAxis/W=$plotName bottom tStep,(len/2)
Label/W=$plotName left "MSD"
Label/W=$plotName bottom "Time ("+unitWave[0]+")"
ErrorBars/W=$plotName $avName SHADE= {0,0,(0,0,0,0),(0,0,0,0)},wave=($errName,$errName)
ModifyGraph/W=$plotName lsize($avName)=2,rgb($avName)=(0,0,0)
AppendLayoutObject /W=$layoutName graph $plotName
// calculate direction autocorrelation
plotName = condName + "_DAplot"
KillWindow/Z $plotName // setup plot
Display/N=$plotName/HIDE=1
for(i = 0; i < nWaves; i += 1)
wName0 = StringFromList(i,wList0) // tk wave
WAVE w0 = $wName0
len = DimSize(w0,0) // len is number of frames
Differentiate/METH=2/DIM=0/EP=1 w0 /D=vWave // make vector wave (vWave). nVectors is len-1
MatrixOp/O/FREE magWave = sqrt(sumrows(vWave * vWave)) // calculate magnitude of each vector
vWave[][] /= magWave[p] // normalise vectors
newName = ReplaceString("tk",wName0,"DA") // for results of DA per cell
Make/O/N=(len-2,len-2,2)/FREE tempDAMat0,tempDAMat1
tempDAMat0[][][] = (p >= q) ? vWave[p-q][r] : 0
tempDAMat1[][][] = (p >= q) ? vWave[p+1][r] : 0
MatrixOp/O/FREE dotWave = (tempDAMat0 * tempDAMat1)
MatrixOp/O/FREE alphaWave = sumBeams(dotWave)
// Make average. Previously we did this:
// Make/O/N=(len-2)/FREE countOfDAPnts = (len-2)-p
// MatrixOp/O $newName = sumcols(alphaWave)^t / countOfDAPnts
// Now we need to get rid of NaNs in the alphaWave and count the non-NaN points
MatrixOp/O/FREE alphaWave = replaceNans(alphaWave,0)
Make/O/FREE/N=(dimsize(alphaWave,0),dimsize(alphaWave,1)) countOfDAPnts
countOfDAPnts[][] = (abs(alphaWave[p][q]) > 0) ? 1 : 0
MatrixOp/O $newName = sumcols(alphaWave)^t / sumCols(countOfDAPnts)^t
SetScale/P x 0,tStep, $newName
AppendtoGraph/W=$plotName $newName
endfor
Killwaves/Z vWave
ModifyGraph/W=$plotName rgb=(colorWave[ii][0],colorWave[ii][1],colorWave[ii][2],alphaLevel)
avList = Wavelist("DA*",";","WIN:"+ plotName)
avName = "W_Ave_DA_" + condName
errName = ReplaceString("Ave", avName, "Err")
fWaveAverage(avList, "", 3, 1, avName, errName)
AppendToGraph/W=$plotName $avName
SetAxis/W=$plotName left -1,1
Label/W=$plotName left "Direction autocorrelation"
Label/W=$plotName bottom "Time ("+unitWave[0]+")"
ErrorBars/W=$plotName $avName SHADE= {0,0,(0,0,0,0),(0,0,0,0)},wave=($errName,$errName)
ModifyGraph/W=$plotName lsize($avName)=2,rgb($avName)=(0,0,0)
AppendLayoutObject/W=$layoutName graph $plotName
// calculate the angle distribution
plotName = condName + "_anglePlot"
KillWindow/Z $plotName // setup plot
Display/N=$plotName/HIDE=1
wList0 = WaveList(condName + "_*",";","") // find all matrices
nWaves = ItemsInList(wList0)
Concatenate/O/NP=0 wList0, allTempW // make long matrix of all tracks
Variable nSteps = dimsize(allTempW,0)
Make/O/N=(nSteps)/FREE tempDistThreshW
Make/O/D/N=(nSteps) angleWave = NaN
// quality filter so that minimal steps are not analysed
tempDistThreshW[] = (allTempW[p][5] > 4 * pxSize) ? 1 : 0
Variable successVar = 0
// this will find angles for all tracks even short tracks
// valid track length has a lower bound defined by user, minimum is 6 points (v. 1.15)
// at tStep = 10 this will find max 4 angles in a track that is not found elsewhere
for(i = 0; i < nSteps; i += 1)
// find a proper step
if(tempDistThreshW[i] == 1)
Make/O/N=(2,2)/FREE matAA,matBB
// set first vector offset to origin (need to transpose later)
matAA[][] = allTempW[p + (i-1)][q+3] - allTempW[i-1][q+3]
for(j = i+1; j < nSteps; j += 1)
if(allTempW[j][1] != allTempW[i][1])
// check there is another proper step from the same cell
successVar = 0
break
elseif(tempDistThreshW[j] == 1)
// find a proper step from same cell and set second vector as this, offset to origin
successVar = 1
matBB[][] = allTempW[p + (j-1)][q+3] - allTempW[j-1][q+3]
break
else
successVar = 0
endif
endfor
if(successVar == 1)
// matrices need transposing
MatrixTranspose matAA
MatrixTranspose matBB
// find cos(theta)
MatrixOp/O/FREE matCC = matAA . matBB / (sqrt(sum(matAA * matAA)) * sqrt(sum(matBB * matBB)))
// angle in radians
AngleWave[i] = acos(matCC[0])
endif
endif
endfor
KillWaves/Z allTempW
// zapnans on AngleWave so I can count valid angles
WaveTransform zapnans AngleWave
newName = "h_angle_" + condName
Make/N=41/O $newName
Histogram/B={0,pi/40,41} angleWave,$newName
AppendToGraph/W=$plotName $newName
ModifyGraph/W=$plotName rgb=(colorWave[ii][0],colorWave[ii][1],colorWave[ii][2])
ModifyGraph/W=$plotName mode=5, hbFill=4
SetAxis/W=$plotName/A/N=1/E=1 left
SetAxis/W=$plotName bottom 0,pi
Make/O/N=5 axisW = {0,pi/4,pi/2,3*pi/4,pi}
// vulgar fractions and pi as Unicode
Make/O/N=5/T axisTW = {"0","\u00BC\u03C0","\u00BD\u03C0","\u00BE\u03C0","\u03C0"}
ModifyGraph/W=$plotName userticks(bottom)={axisW,axisTW}
Label/W=$plotName left "Density"
Label/W=$plotName bottom "Cell turning"
AppendLayoutObject/W=$layoutName graph $plotName
// print message about number of angles
Print numpnts(AngleWave), "valid angles found from all tracks for", condName
// Plot these averages to summary windows at the end
avName = "W_Ave_cd_" + condName
errName = ReplaceString("Ave", avName, "Err")
AppendToGraph/W=cdPlot $avName
ErrorBars/W=cdPlot $avName SHADE= {0,0,(0,0,0,0),(0,0,0,0)},wave=($errName,$errName)
ModifyGraph/W=cdPlot lsize($avName)=2,rgb($avName)=(colorWave[ii][0],colorWave[ii][1],colorWave[ii][2])
avName = "W_Ave_iv_" + condName
errName = ReplaceString("Ave", avName, "Err")
AppendToGraph/W=ivPlot $avName
ErrorBars/W=ivPlot $avName SHADE= {0,0,(0,0,0,0),(0,0,0,0)},wave=($errName,$errName)
ModifyGraph/W=ivPlot lsize($avName)=2,rgb($avName)=(colorWave[ii][0],colorWave[ii][1],colorWave[ii][2])
newName = "h_iv_" + condName
AppendToGraph/W=ivHPlot $newName
ModifyGraph/W=ivHPlot rgb($newName)=(colorWave[ii][0],colorWave[ii][1],colorWave[ii][2])
avName = "W_Ave_dD_" + condName
errName = ReplaceString("Ave", avName, "Err")
AppendToGraph/W=dDPlot $avName
ErrorBars/W=dDPlot $avName SHADE= {0,0,(0,0,0,0),(0,0,0,0)},wave=($errName,$errName)
ModifyGraph/W=dDPlot lsize($avName)=2,rgb($avName)=(colorWave[ii][0],colorWave[ii][1],colorWave[ii][2])
avName = "W_Ave_MSD_" + condName
errName = ReplaceString("Ave", avName, "Err")
AppendToGraph/W=MSDPlot $avName
ErrorBars/W=MSDPlot $avName SHADE= {0,0,(0,0,0,0),(0,0,0,0)},wave=($errName,$errName)
ModifyGraph/W=MSDPlot lsize($avName)=2,rgb($avName)=(colorWave[ii][0],colorWave[ii][1],colorWave[ii][2])
avName = "W_Ave_DA_" + condName
errName = ReplaceString("Ave", avName, "Err")
AppendToGraph/W=DAPlot $avName
ErrorBars/W=DAPlot $avName SHADE= {0,0,(0,0,0,0),(0,0,0,0)},wave=($errName,$errName)
ModifyGraph/W=DAPlot lsize($avName)=2,rgb($avName)=(colorWave[ii][0],colorWave[ii][1],colorWave[ii][2])
newName = "h_angle_" + condName
AppendToGraph/W=angleHPlot $newName
ModifyGraph/W=angleHPlot rgb($newName)=(colorWave[ii][0],colorWave[ii][1],colorWave[ii][2])
End
/// @param plotString String such as "*_tkPlot" to determine which plot get rescaled
/// @param limitVar Variable containing the furthest point
STATIC Function RescalePlotsToLimits(plotString,limitVar)
String plotString
Variable limitVar
// previously a temporary limit of -250,250 was set to all tkPlots and then expanded here if necessary
// as of v 1.15, we just use the furthest point
if(limitVar > 1000)
limitVar = ceil(limitVar / 100) * 100
elseif(limitVar >= 10 && limitVar < 100)
limitVar = ceil(limitVar / 10) * 10
elseif(limitVar >= 1 && limitVar < 10)
limitVar = ceil(limitVar / 1) * 1
elseif(limitVar >= 0.1 && limitVar < 1)
limitVar = ceil(limitVar / 0.1) * 0.1
elseif(limitVar < 0.1)
limitVar = limitVar
else
limitVar = ceil(limitVar / 50) * 50 // we expect values of 150-300
endif
String plotList = WinList(plotString,";","WIN:1")
String plotName
Variable i
for(i = 0; i < ItemsInList(plotList); i += 1)
plotName = StringFromList(i,plotList)
SetAxis/W=$plotName left -limitVar,limitVar
SetAxis/W=$plotName/Z bottom -limitVar,limitVar
SetAxis/W=$plotName/Z top -limitVar,limitVar
endfor
return 1
End
// This function will make ImageQuilts and sparklines of 2D tracks for all conditions
/// @param qSize Variable to indicate desired size of image quilt (qSize^2 tracks)
Function MakeImageQuilt(qSize)
Variable qSize
WAVE/T condWave = root:condWave
Variable cond = numpnts(condWave)
Wave paramWave = root:paramWave
Variable tStep = paramWave[1]
Variable pxSize = paramWave[2]
Wave colorWave = root:colorWave
String condName, dataFolderName, wName
Variable longestCond = 0 , mostFrames = 0
Variable i
for(i = 0; i < cond; i += 1)
condName = condWave[i]
dataFolderName = "root:data:" + condName
SetDataFolder datafolderName
mostFrames = FindSolution()
longestCond = max(longestCond,mostFrames)
endfor
SetDataFolder root:
// Now they're all done, cycle again to find optimum quilt size
Make/O/N=(longestCond,qSize+1,cond)/FREE optiMat
for(i = 0; i < cond; i += 1)
condName = condWave[i]
wName = "root:data:" + condName + ":solutionWave"
Wave w0 = $wName
optiMat[][][i] = (w0[p][0] >= q^2) ? 1 : 0
endfor
optiMat /= cond
// make a 1D wave where row = qSize and value = frames that can be plotted for all cond
MatrixOp/O/FREE quiltSizeMat = sumcols(floor(sumBeams(optiMat)))^t
// find optimum
quiltSizeMat *= p^2
WaveStats/Q quiltSizeMat
Variable optQSize = V_maxRowLoc
Variable optDur = (V_max / V_maxRowLoc^2) - 1 // because 0-based
Print qSize, "x", qSize, "quilt requested.", optQSize, "x", optQSize, "quilt with", optDur, "frames shown (", optDur * tStep, "min)."
String plotName,sampleWName
Variable startVar,endVar,xShift,yShift
Variable spaceVar = WorkOutSpacing(pxSize)
Variable j,k
for(i = 0; i < cond; i += 1)
condName = condWave[i]
dataFolderName = "root:data:" + condName
SetDataFolder datafolderName
// make image quilt for each condition
plotName = condName + "_quilt"
KillWindow/Z $plotName
Display/N=$plotName/HIDE=1
WAVE segValid, trackDurations
WAVE/T trackNames
segValid[] = (trackDurations[p] > optDur * tStep) ? p : NaN
WaveTransform zapnans segValid
StatsSample/N=(optQSize^2) segValid
WAVE/Z W_Sampled
sampleWName = "segSelected_" + condName
Duplicate/O W_Sampled, $sampleWName
Wave segSelected = $sampleWName
Make/O/N=(optQSize^2)/T segNames
Make/O/N=(optQSize^2) segLengths
for(j = 0; j < optQSize^2; j += 1)
segNames[j] = trackNames[segSelected[j]]
wName = ReplaceString("tk_",segNames[j],"cd_") // get cum dist wave name
Wave cdW0 = $wName
segLengths[j] = cdW0[optDur] // store cum dist at point optDur
endfor
Sort segLengths, segLengths, segNames
// plot segNamed waves out
Make/O/N=(optQSize^2*(optDur+1),2) quiltBigMat = NaN
for(j = 0; j < optQSize^2; j += 1)
wName = segNames[j]
Wave tkW0 = $wName
// put each track into the big quilt wave leaving a NaN between each
startVar = j * optDur + (j * 1)
endVar = startVar + optDur - 1
quiltBigMat[startVar,endVar][] = tkW0[p-startVar][q]
xShift = mod(j,optQSize) * spaceVar
yShift = floor(j/optQSize) * spaceVar
quiltBigMat[startVar,endVar][0] += xShift
quiltBigMat[startVar,endVar][1] += yShift
endfor
// Add to plot and then format
AppendToGraph/W=$plotName quiltBigMat[][1] vs quiltBigMat[][0]
SetAxis/W=$plotName left (optQsize+1.5) * spaceVar,-2.5*spaceVar
SetAxis/W=$plotName bottom -2.5*spaceVar,(optQsize+1.5) * spaceVar
ModifyGraph/W=$plotName width={Aspect,1}
ModifyGraph/W=$plotName manTick={0,spaceVar,0,0},manMinor={0,0}
ModifyGraph/W=$plotName noLabel=2,mirror=1,standoff=0,tick=3
ModifyGraph/W=$plotName grid=1,gridRGB=(34952,34952,34952)
ModifyGraph axRGB=(65535,65535,65535)
ModifyGraph/W=$plotName rgb=(colorWave[i][0],colorWave[i][1],colorWave[i][2])
ModifyGraph/W=$plotName margin=14
// Append to appropriate layout (page 2)
String layoutName = condName + "_layout"
LayoutPageAction/W=$layoutName appendPage
AppendLayoutObject/W=$layoutName/PAGE=(2) graph $plotName
ModifyLayout/W=$layoutName/PAGE=(2) left($plotName)=21,top($plotName)=291,width($plotName)=261,height($plotName)=261
// make sparkline graphic
plotName = condName + "_sprkln"
KillWindow/Z $plotName
Display/N=$plotName/HIDE=1
// plot the diagonal of segNamed waves out laterally
Make/O/N=(optQSize*(optDur+1),2) sprklnBigMat = NaN
Variable theta
k = 0
for(j = 0; j < optQSize^2; j += optQsize+1)
wName = segNames[j]
Wave tkW0 = $wName
Duplicate/O/FREE tkW0, sprkW0
theta = (1.5 * pi) - atan2(sprkW0[optDur-1][1],sprkW0[optDur-1][0])
Make/O/N=(2,2)/FREE rotMat = {{cos(theta),-sin(theta)},{sin(theta),cos(theta)}}
MatrixMultiply sprkW0, rotMat
WAVE/Z M_product
// put each track into the big quilt wave leaving a NaN between each
startVar = k * optDur + (k * 1)
endVar = startVar + optDur - 1
sprklnBigMat[startVar,endVar][] = M_Product[p-startVar][q]
xShift = mod(k,optQSize) * spaceVar
sprklnBigMat[startVar,endVar][0] += xShift
k += 1
endfor
KillWaves/Z M_product
// Add to plot and then format
AppendToGraph/W=$plotName sprklnBigMat[][1] vs sprklnBigMat[][0]
SetAxis/W=$plotName left 0.75 * spaceVar,-2.5*spaceVar
SetAxis/W=$plotName bottom -0.5*spaceVar,(optQsize+0.5) * spaceVar
ModifyGraph/W=$plotName width={Plan,1,bottom,left}
ModifyGraph/W=$plotName manTick={0,spaceVar,0,0},manMinor={0,0}
ModifyGraph/W=$plotName noLabel=2,mirror=1,standoff=0,tick=3
ModifyGraph/W=$plotName grid=1,gridRGB=(34952,34952,34952)
ModifyGraph axRGB=(65535,65535,65535)
ModifyGraph/W=$plotName rgb=(colorWave[i][0],colorWave[i][1],colorWave[i][2])
ModifyGraph/W=$plotName margin=14
SetDrawLayer/W=$plotName UserBack