/
rgcAnalyzer.m
1918 lines (1756 loc) · 80.9 KB
/
rgcAnalyzer.m
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
% Various MATLAB functions used in the paper
% 'A genetic and computational approach to structural classification of neuronal types'
% Software developed by: Uygar Sümbül <uygar@mit.edu>
% THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE.
% IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DAMAGES WHATSOEVER.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rgcAnalyzer(arborFileName,OnSACFilename,OffSACFilename,voxelRes,conformalJump)
% read the arbor trace file - add 1 to node positions because FIJI format for arbor tracing starts from 0
[nodes,edges,radii,nodeTypes,abort] = readArborTrace(arborFileName,[-1 0 1 2 3 4 5]); nodes = nodes + 1;
arborBoundaries(1) = min(nodes(:,1)); arborBoundaries(2) = max(nodes(:,1));
arborBoundaries(3) = min(nodes(:,2)); arborBoundaries(4) = max(nodes(:,2));
%generate the SAC surfaces from annotations
thisVZminmesh = fitSurfaceToSACAnnotation(OnSACFilename);
thisVZmaxmesh = fitSurfaceToSACAnnotation(OffSACFilename);
% find conformal maps of the ChAT surfaces onto the median plane
surfaceMapping = calcWarpedSACsurfaces(thisVZminmesh,thisVZmaxmesh,arborBoundaries,conformalJump);
warpedArbor = calcWarpedArbor(nodes,edges,radii,surfaceMapping,voxelRes,conformalJump);
% hard-coded representation and filtering parameters used in arbor density representation (ADR) calculation
% filter parameters were optimized based on representation size choices to minimize discrepancy with genetic
% labels as explained in the manuscript
zLen=120; % number of pixels in z of the ADR
repLenX=20; % number of pixels in x (and y) of the ADR
zExt=5; % z-extent of the ADR in terms of the distance between SAC surfaces
xExt=420; % x-extent (and y-extent) of the ADR in um (x resolution of ADR=xExtent/replenX (in um))
xLen=720; % x-extent (and y-extent) of the initial nearest neighbor gridding (in um)
fLR=0.6; % ratio of the extent of the post-registration in-plane low-pass filter to representation length (repLenX)
betaX=11;betaY=29; % beta parameters for the Kaiser-Bessel filters used in post-registration in-plane low-pass filtering
% calculate the anisotropic post-registration in-plane low-pass filter for individual dimensions (x and y)
postFilterWx = ceil(fLR*repLenX); x=-(repLenX/2):(repLenX/2)-1;
kbX=sin(sqrt((pi*postFilterWx*x/repLenX).^2-betaX^2))./sqrt((pi*postFilterWx*x/repLenX).^2-betaX^2);
kbY=sin(sqrt((pi*postFilterWx*x/repLenX).^2-betaY^2))./sqrt((pi*postFilterWx*x/repLenX).^2-betaY^2);
kbX(isnan(kbX))=0; kbX=kbX/max(kbX); kbY(isnan(kbY))=0; kbY=kbY/max(kbY);
% calculate the separable 2-d (xy) filter and rotate it 45 degrees because the longer extent of the
% arbor trace is aligned with the main diagonal in returnDist3d for storage/computation efficiency
xyFilter = imrotate(kbX'*kbY,45,'bicubic','crop');lpfilt=permute(repmat(xyFilter,[1 1 zLen]),[3 1 2]);
% calculate the ADR and use nodes reassigned to edge centers weighted by edge masses
[dist3d,nodes,density] = calc3dDist(warpedArbor.nodes,warpedArbor.edges,warpedArbor.medVZmin,warpedArbor.medVZmax,repLenX,zLen,xLen,zExt,xExt,lpfilt);
% calculate depth profile directly from the arbor and save as png image
gridPointCount = 120; % number of points on the z-profile plot
[token,remain] = strtok(datestr(clock),' '); commonSaveFileName = strcat(arborFileName,token,'_',remain(2:end));
saveZProfile(warpedArbor.medVZmin,warpedArbor.medVZmax,density,nodes,voxelRes(3),zExt,gridPointCount,commonSaveFileName);
%% load allDist here -- all arbor distributions for classification
%% to find how the new cell will be clustered
% myLinkage = elinkage(pdist([allDist dist3d(:)]','euclidean')); myLinkage(:,3) = myLinkage(:,3)/myLinkage(end,3);
% labelDendrogram(myLinkage,readLabels,commonSaveFileName);
% [AR,RI,MI,HI,rowSplits,columnSplits,typeConfusions] = reportConfusionsAndRI(c1,c2);
% hardClusterStats = cutDendrogram(myLinkage,clusterIDsForKnownCells,knownCellPositions,maxCutLevel);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [nodes,edges,radii,nodeTypes,abort] = readArborTrace(fileName,validNodeTypes)
abort = false; nodes = []; edges = []; nodeTypes = [];
validNodeTypes = setdiff(validNodeTypes,1); % 1 is for soma
% read the SWC file
[nodeID, nodeType, xPos, yPos, zPos, radii, parentNodeID] = textread(fileName, '%u%d%f%f%f%f%d','commentstyle', 'shell');
% every tree should start from a node of type 1 (soma)
nodeType(find(parentNodeID==-1))=1;
% find the first soma node in the list (more than one node can be labeled as soma)
firstSomaNode = find(nodeType == 1 & parentNodeID == -1, 1);
% find the average position of all the soma nodes, and assign it as THE soma node position
somaNodes = find(nodeType == 1);
somaX = mean(xPos(somaNodes)); somaY = mean(yPos(somaNodes)); somaZ = mean(zPos(somaNodes));
somaRadius = mean(radii(somaNodes));
xPos(firstSomaNode) = somaX; yPos(firstSomaNode) = somaY; zPos(firstSomaNode) = somaZ;
radii(firstSomaNode) = somaRadius;
% change parenthood so that there is a single soma parent
parentNodeID(ismember(parentNodeID,somaNodes)) = firstSomaNode;
% delete all the soma nodes except for the firstSomaNode
nodesToDelete = setdiff(somaNodes,firstSomaNode);
nodeID(nodesToDelete)=[]; nodeType(nodesToDelete)=[];
xPos(nodesToDelete)=[]; yPos(nodesToDelete)=[]; zPos(nodesToDelete)=[];
radii(nodesToDelete)=[]; parentNodeID(nodesToDelete)=[];
% reassign node IDs due to deletions
for kk = 1:numel(nodeID)
while ~any(nodeID==kk)
nodeID(nodeID>kk) = nodeID(nodeID>kk)-1;
parentNodeID(parentNodeID>kk) = parentNodeID(parentNodeID>kk)-1;
end
end
% of all the nodes, retain the ones indicated in validNodeTypes
% ensure connectedness of the tree if a child is marked as valid but not some of its ancestors
validNodes = nodeID(ismember(nodeType,validNodeTypes));
additionalValidNodes = [];
for kk = 1:numel(validNodes)
thisParentNodeID = parentNodeID(validNodes(kk)); thisParentNodeType = nodeType(thisParentNodeID);
while ~ismember(thisParentNodeType,validNodeTypes)
if thisParentNodeType == 1
break;
end
additionalValidNodes = union(additionalValidNodes, thisParentNodeID); nodeType(thisParentNodeID) = validNodeTypes(1);
thisParentNodeID = parentNodeID(thisParentNodeID); thisParentNodeType = nodeType(thisParentNodeID);
end
end
% retain the valid nodes only
% the soma node is always a valid node to ensure connectedness of the tree
validNodes = [firstSomaNode; validNodes; additionalValidNodes']; validNodes = unique(validNodes);
nodeID = nodeID(validNodes); nodeType = nodeType(validNodes); parentNodeID = parentNodeID(validNodes);
xPos = xPos(validNodes); yPos = yPos(validNodes); zPos = zPos(validNodes); radii = radii(validNodes);
% reassign node IDs after deletions
for kk = 1:numel(nodeID)
while ~any(nodeID==kk)
nodeID(nodeID>kk) = nodeID(nodeID>kk)-1;
parentNodeID(parentNodeID>kk) = parentNodeID(parentNodeID>kk)-1;
end
end
% return the resulting tree data
nodes = [xPos yPos zPos];
edges = [nodeID parentNodeID];
edges(any(edges==-1,2),:) = [];
nodeTypes = unique(nodeType)';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function vzmesh = fitSurfaceToSACAnnotation(annotationFilename)
% read the annotation file
[t1,t2,t3,t4,t5,x,z,y] = textread(annotationFilename, '%u%u%d%d%d%d%d%d','headerlines',1);
% add 1 to x and z, butnot to y because FIJI point tool starts from 0 for pixels of an image
% but from 1 for slice of a stack
x=x+1; z=z+1;
% find the maximum boundaries
xMax = max(x); yMax = max(y);
% smoothened fit (with extrapolation) to a grid - to save time, make the grid coarser (every 3 pixels)
[zgrid,xgrid,ygrid] = gridfit(x,y,z,[[1:3:xMax-1] xMax],[[1:3:yMax-1] yMax],'smoothness',1);
% linearly (fast) interpolate to fine grid
[xi,yi]=meshgrid(1:xMax,1:yMax); xi = xi'; yi = yi';
vzmesh=interp2(xgrid,ygrid,zgrid,xi,yi,'*linear',mean(zgrid(:)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function surfaceMapping = calcWarpedSACsurfaces(thisVZminmesh,thisVZmaxmesh,arborBoundaries,conformalJump)
minXpos = arborBoundaries(1); maxXpos = arborBoundaries(2); minYpos = arborBoundaries(3); maxYpos = arborBoundaries(4);
% retain the minimum grid of SAC points, where grid resolution is determined by conformalJump
thisx = [max(minXpos-1,1):conformalJump:min(maxXpos+1,size(thisVZmaxmesh,1))];
thisy = [max(minYpos-1,1):conformalJump:min(maxYpos+1,size(thisVZmaxmesh,2))];
thisminmesh=thisVZminmesh(thisx,thisy); thismaxmesh=thisVZmaxmesh(thisx,thisy);
% calculate the traveling distances on the diagonals of the two SAC surfaces - this must be changed to Dijkstra's algorithm for exact results
[mainDiagDistMin, skewDiagDistMin] = calculateDiagLength(thisx,thisy,thisminmesh);
[mainDiagDistMax, skewDiagDistMax] = calculateDiagLength(thisx,thisy,thismaxmesh);
% average the diagonal distances on both surfaces for more stability against band tracing errors - not ideal
mainDiagDist = (mainDiagDistMin+mainDiagDistMax)/2; skewDiagDist = (skewDiagDistMin+skewDiagDistMax)/2;
% quasi-conformally map individual SAC surfaces to planes
mappedMinPositions = conformalMap_indepFixedDiagonals(mainDiagDist,skewDiagDist,thisx,thisy,thisminmesh);
mappedMaxPositions = conformalMap_indepFixedDiagonals(mainDiagDist,skewDiagDist,thisx,thisy,thismaxmesh);
% align the two independently mapped surfaces so their flattest regions are registered to each other
xborders = [thisx(1) thisx(end)]; yborders = [thisy(1) thisy(end)];
mappedMaxPositions = alignMappedSurfaces(thisVZminmesh,thisVZmaxmesh,mappedMinPositions,mappedMaxPositions,xborders,yborders,conformalJump);
% return original and mapped surfaces with grid information
surfaceMapping.mappedMinPositions=mappedMinPositions; surfaceMapping.mappedMaxPositions=mappedMaxPositions;;
surfaceMapping.thisVZminmesh=thisVZminmesh; surfaceMapping.thisVZmaxmesh=thisVZmaxmesh; surfaceMapping.thisx=thisx; surfaceMapping.thisy=thisy;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [mainDiagDist, skewDiagDist] = calculateDiagLength(xpos,ypos,VZmesh)
M = size(VZmesh,1); N = size(VZmesh,2);
[ymesh,xmesh] = meshgrid(ypos,xpos);
mainDiagDist = 0; skewDiagDist = 0;
% travel on the diagonals and not necessarily the grid points) and accumulate the 3d distance traveled
if N >= M
xKnots = interp2(ymesh,xmesh,xmesh, ypos, [xpos(1):(xpos(end)-xpos(1))/(N-1):xpos(end)]');
yKnots = interp2(ymesh,xmesh,ymesh, ypos, [xpos(1):(xpos(end)-xpos(1))/(N-1):xpos(end)]');
zKnotsMainDiag = griddata(xmesh(:),ymesh(:),VZmesh(:), [xpos(1):(xpos(end)-xpos(1))/(N-1):xpos(end)]', ypos');
zKnotsSkewDiag = griddata(xmesh(:),ymesh(:),VZmesh(:), [xpos(1):(xpos(end)-xpos(1))/(N-1):xpos(end)]', ypos(end:-1:1)');
dxKnots = diag(xKnots); dyKnots = diag(yKnots); mainDiagDist = sum(sqrt(diff(dxKnots).^2 + diff(dyKnots).^2 + diff(zKnotsMainDiag).^2));
sdxKnots = xKnots(N:N-1:end-1)'; sdyKnots = yKnots(N:N-1:end-1)'; skewDiagDist = sum(sqrt(diff(sdxKnots).^2 + diff(sdyKnots).^2 + diff(zKnotsSkewDiag).^2));
else
xKnots = interp2(ymesh,xmesh,xmesh, [ypos(1):(ypos(end)-ypos(1))/(M-1):ypos(end)], xpos');
yKnots = interp2(ymesh,xmesh,ymesh, [ypos(1):(ypos(end)-ypos(1))/(M-1):ypos(end)], xpos');
zKnotsMainDiag = griddata(xmesh(:),ymesh(:),VZmesh(:), xpos', [ypos(1):(ypos(end)-ypos(1))/(M-1):ypos(end)]');
zKnotsSkewDiag = griddata(xmesh(:),ymesh(:),VZmesh(:), xpos', [ypos(end):-(ypos(end)-ypos(1))/(M-1):ypos(1)]');
dxKnots = diag(xKnots); dyKnots = diag(yKnots); mainDiagDist = sum(sqrt(diff(dxKnots).^2 + diff(dyKnots).^2 + diff(zKnotsMainDiag).^2));
sdxKnots = xKnots(M:M-1:end-1)'; sdyKnots = yKnots(M:M-1:end-1)'; skewDiagDist = sum(sqrt(diff(sdxKnots).^2 + diff(sdyKnots).^2 + diff(zKnotsSkewDiag).^2));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mappedMaxPositions = alignMappedSurfaces(VZminmesh,VZmaxmesh,mappedMinPositions,mappedMaxPositions,validXborders,validYborders,conformalJump,patchSize)
if nargin < 8
patchSize = 21;
if nargin < 7
conformalJump = 1;
end
end
patchSize = ceil(patchSize/conformalJump);
% pad the surfaces by one pixel so that the size remains the same after the difference operation
VZminmesh = [[VZminmesh 10*max(VZminmesh(:))*ones(size(VZminmesh,1),1)]; 10*max(VZminmesh(:))*ones(1,size(VZminmesh,2)+1)];
VZmaxmesh = [[VZmaxmesh 10*max(VZmaxmesh(:))*ones(size(VZmaxmesh,1),1)]; 10*max(VZmaxmesh(:))*ones(1,size(VZmaxmesh,2)+1)];
% calculate the absolute differences in xy between neighboring pixels
tmp1 = diff(VZminmesh,1,1); tmp1 = tmp1(:,1:end-1); tmp2 = diff(VZminmesh,1,2); tmp2 = tmp2(1:end-1,:); dMinSurface = abs(tmp1+i*tmp2);
tmp1 = diff(VZmaxmesh,1,1); tmp1 = tmp1(:,1:end-1); tmp2 = diff(VZmaxmesh,1,2); tmp2 = tmp2(1:end-1,:); dMaxSurface = abs(tmp1+i*tmp2);
% retain the region of interest, with a resolution specified by conformalJump
dMinSurface = dMinSurface(validXborders(1):conformalJump:validXborders(2), validYborders(1):conformalJump:validYborders(2));
dMaxSurface = dMaxSurface(validXborders(1):conformalJump:validXborders(2), validYborders(1):conformalJump:validYborders(2));
% calculate the cost as the sum of absolute slopes on both surfaces
patchCosts = conv2(dMinSurface+dMaxSurface, ones(patchSize),'valid');
% find the minimum cost
[row,col] = find(patchCosts == min(min(patchCosts)),1);
row = row+(patchSize-1)/2; col = col+(patchSize-1)/2;
% calculate and apply the corresponding shift
linearInd = sub2ind(size(dMinSurface),row,col);
shift1 = mappedMaxPositions(linearInd,1)-mappedMinPositions(linearInd,1);
mappedMaxPositions(:,1) = mappedMaxPositions(:,1) - shift1;
shift2 = mappedMaxPositions(linearInd,2)-mappedMinPositions(linearInd,2);
mappedMaxPositions(:,2) = mappedMaxPositions(:,2) - shift2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mappedPositions = conformalMap_indepFixedDiagonals(mainDiagDist,skewDiagDist,xpos,ypos,VZmesh)
% implements quasi-conformal mapping suggested in
% Levy et al., 'Least squares conformal maps for automatic texture atlas generation', 2002, ACM Transactions on Graphics
M = size(VZmesh,1); N = size(VZmesh,2);
col1 = kron([1;1],[1:M-1]');
temp1 = kron([1;M+1],ones(M-1,1));
temp2 = kron([M+1;M],ones(M-1,1));
oneColumn = [col1 col1+temp1 col1+temp2];
% every pixel is divided into 2 triangles
triangleCount = (2*M-2)*(N-1);
vertexCount = M*N;
triangulation = zeros(triangleCount, 3);
% store the positions of the vertices for each triangle - triangles are oriented consistently
for kk = 1:N-1
triangulation((kk-1)*(2*M-2)+1:kk*(2*M-2),:) = oneColumn + (kk-1)*M;
end
Mreal = sparse([],[],[],triangleCount,vertexCount,triangleCount*3);
Mimag = sparse([],[],[],triangleCount,vertexCount,triangleCount*3);
% calculate the conformality condition (Riemann's theorem)
for triangle = 1:triangleCount
for vertex = 1:3
nodeNumber = triangulation(triangle,vertex);
xind = rem(nodeNumber-1,M)+1;
yind = floor((nodeNumber-1)/M)+1;
trianglePos(vertex,:) = [xpos(xind) ypos(yind) VZmesh(xind,yind)];
end
[w1,w2,w3,zeta] = assignLocalCoordinates(trianglePos);
denominator = sqrt(zeta/2);
Mreal(triangle,triangulation(triangle,1)) = real(w1)/denominator;
Mreal(triangle,triangulation(triangle,2)) = real(w2)/denominator;
Mreal(triangle,triangulation(triangle,3)) = real(w3)/denominator;
Mimag(triangle,triangulation(triangle,1)) = imag(w1)/denominator;
Mimag(triangle,triangulation(triangle,2)) = imag(w2)/denominator;
Mimag(triangle,triangulation(triangle,3)) = imag(w3)/denominator;
end
% minimize the LS error due to conformality condition of mapping triangles into triangles
% take the two fixed points required to solve the system as the corners of the main diagonal
mainDiagXdist = mainDiagDist*M/sqrt(M^2+N^2); mainDiagYdist = mainDiagXdist*N/M;
A = [Mreal(:,2:end-1) -Mimag(:,2:end-1); Mimag(:,2:end-1) Mreal(:,2:end-1)];
b = -[Mreal(:,[1 end]) -Mimag(:,[1 end]); Mimag(:,[1 end]) Mreal(:,[1 end])]*[[xpos(1);xpos(1)+mainDiagXdist];[ypos(1);ypos(1)+mainDiagYdist]];
mappedPositions1 = A\b;
mappedPositions1 = [[xpos(1);mappedPositions1(1:end/2);xpos(1)+mainDiagXdist] [ypos(1);mappedPositions1(1+end/2:end);ypos(1)+mainDiagYdist]];
% take the two fixed points required to solve the system as the corners of the skew diagonal
skewDiagXdist = skewDiagDist*M/sqrt(M^2+N^2); skewDiagYdist = skewDiagXdist*N/M; freeVar = [[1:M-1] [M+1:M*N-M] [M*N-M+2:M*N]];
A = [Mreal(:,freeVar) -Mimag(:,freeVar); Mimag(:,freeVar) Mreal(:,freeVar)];
b = -[Mreal(:,[M M*N-M+1]) -Mimag(:,[M M*N-M+1]); Mimag(:,[M M*N-M+1]) Mreal(:,[M M*N-M+1])]*[[xpos(1)+skewDiagXdist;xpos(1)];[ypos(1);ypos(1)+skewDiagYdist]];
mappedPositions2 = A\b;
mappedPositions2 = [[mappedPositions2(1:M-1);xpos(1)+skewDiagXdist;mappedPositions2(M:M*N-M-1);xpos(1);mappedPositions2(M*N-M:end/2)] ...
[mappedPositions2(1+end/2:M-1+end/2);ypos(1);mappedPositions2(end/2+M:end/2+M*N-M-1);ypos(1)+skewDiagYdist;mappedPositions2(end/2+M*N-M:end)]];
% take the mean of the two independent solutions and return it as the solution
mappedPositions = (mappedPositions1+mappedPositions2)/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [w1,w2,w3,zeta] = assignLocalCoordinates(triangle)
% triangle is a 3x3 matrix where the rows represent the x,y,z coordinates of the vertices.
% The local triangle is defined on a plane, where the first vertex is at the origin(0,0) and the second vertex is at (0,-d12).
d12 = norm(triangle(1,:)-triangle(2,:));
d13 = norm(triangle(1,:)-triangle(3,:));
d23 = norm(triangle(2,:)-triangle(3,:));
y3 = ((-d12)^2+d13^2-d23^2)/(2*(-d12));
x3 = sqrt(d13^2-y3^2); % provisional
w2 = -x3-i*y3; w1 = x3+i*(y3-(-d12));
zeta = i*(conj(w2)*w1-conj(w1)*w2); % orientation indicator
w3 = i*(-d12);
zeta = abs(real(zeta));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function warpedArbor = calcWarpedArbor(nodes,edges,radii,surfaceMapping,voxelDim,conformalJump)
% voxelDim: physical size of voxels in um, 1x3
mappedMinPositions=surfaceMapping.mappedMinPositions; mappedMaxPositions=surfaceMapping.mappedMaxPositions;
thisVZminmesh=surfaceMapping.thisVZminmesh; thisVZmaxmesh=surfaceMapping.thisVZmaxmesh; thisx=surfaceMapping.thisx; thisy=surfaceMapping.thisy;
% generate correspondence points for the points on the surfaces
[tmpymesh,tmpxmesh] = meshgrid([thisy(1):conformalJump:thisy(end)],[thisx(1):conformalJump:thisx(end)]);
tmpminmesh = thisVZminmesh(thisx(1):conformalJump:thisx(end),thisy(1):conformalJump:thisy(end));
tmpmaxmesh = thisVZmaxmesh(thisx(1):conformalJump:thisx(end),thisy(1):conformalJump:thisy(end));
topInputPos = [tmpxmesh(:) tmpymesh(:) tmpminmesh(:)]; botInputPos = [tmpxmesh(:) tmpymesh(:) tmpmaxmesh(:)];
topOutputPos = [mappedMinPositions(:,1) mappedMinPositions(:,2) median(tmpminmesh(:))*ones(size(mappedMinPositions,1),1)];
botOutputPos = [mappedMaxPositions(:,1) mappedMaxPositions(:,2) median(tmpmaxmesh(:))*ones(size(mappedMaxPositions,1),1)];
% use the correspondence points to calculate local transforms and use those local transforms to map points on the arbor
nodes = localLSregistration(nodes,topInputPos,botInputPos,topOutputPos,botOutputPos,conformalJump);
% switch to physical dimensions (in um)
nodes(:,1) = nodes(:,1)*voxelDim(1); nodes(:,2) = nodes(:,2)*voxelDim(2); nodes(:,3) = nodes(:,3)*voxelDim(3);
% calculate median band positions in z
medVZminmesh = median(tmpminmesh(:)); medVZmaxmesh = median(tmpmaxmesh(:));
% return the warped arbor and the corresponding median SAC surface values
warpedArbor.nodes=nodes; warpedArbor.edges=edges; warpedArbor.radii=radii;
warpedArbor.medVZmin=medVZminmesh; warpedArbor.medVZmax=medVZmaxmesh;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function nodes = localLSregistration(nodes,topInputPos,botInputPos,topOutputPos,botOutputPos,conformalJump)
% aX^2+bXY+cY^2+dX+eY+f+gZX^2+hXYZ+iZY^2+jXZ+kYZ+lZ -> at least 12 equations needed
maxOrder=2; % maximum multinomial order in xy
for kk = 1:size(nodes,1)
window=conformalJump; % neighborhood extent for the LS registration
% gradually increase 'window' to avoid an underdetermined system
while(true)
% fetch the band points around an xy neighborhood for each point on the arbor
xpos = nodes(kk,1); ypos = nodes(kk,2); zpos = nodes(kk,3);
lx = round(xpos-window); ux = round(xpos+window); ly = round(ypos-window); uy = round(ypos+window);
thisInT = topInputPos(topInputPos(:,1)>=lx & topInputPos(:,1)<=ux & topInputPos(:,2)>=ly & topInputPos(:,2)<=uy,:);
thisInB = botInputPos(botInputPos(:,1)>=lx & botInputPos(:,1)<=ux & botInputPos(:,2)>=ly & botInputPos(:,2)<=uy,:);
thisIn = [thisInT; thisInB];
thisOutT = topOutputPos(topInputPos(:,1)>=lx & topInputPos(:,1)<=ux & topInputPos(:,2)>=ly & topInputPos(:,2)<=uy,:);
thisOutB = botOutputPos(botInputPos(:,1)>=lx & botInputPos(:,1)<=ux & botInputPos(:,2)>=ly & botInputPos(:,2)<=uy,:);
thisOut = [thisOutT; thisOutB];
% convert band correspondence data into local coordinates
xShift = mean(thisIn(:,1)); yShift = mean(thisIn(:,2));
thisIn(:,1) = thisIn(:,1)-xShift; thisOut(:,1) = thisOut(:,1)-xShift;
thisIn(:,2) = thisIn(:,2)-yShift; thisOut(:,2) = thisOut(:,2)-yShift;
% calculate the transformation that maps the band points to their correspondences
quadData = [thisIn(:,1:2) ones(size(thisIn,1),1)];
for totalOrd = 2:maxOrder; for ord = 0:totalOrd; % slightly more general - it can handle higher lateral polynomial orders
quadData = [(thisIn(:,1).^ord).*(thisIn(:,2).^(totalOrd-ord)) quadData];
end; end;
quadData = [quadData kron(thisIn(:,3),ones(1,size(quadData,2))).*quadData];
if rank(quadData)<12
window=window+1;
else
break;
end
end
transformMat = lscov(quadData,thisOut);
shiftedNode = [xpos-xShift ypos-yShift];
quadData = [shiftedNode(1:2) 1];
for totalOrd = 2:maxOrder; for ord = 0:totalOrd; % slightly more general - it can handle higher lateral polynomial orders
quadData = [(shiftedNode(1).^ord).*(shiftedNode(2).^(totalOrd-ord)) quadData];
end; end;
quadData = [quadData kron(zpos,ones(1,size(quadData,2))).*quadData];
% transform the arbor points using the calculated transform
nodes(kk,:) = quadData*transformMat; nodes(kk,1) = nodes(kk,1) + xShift; nodes(kk,2) = nodes(kk,2) + yShift;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dist3d,nodes,density] = calc3dDist(nodes,edges,medVZmin,medVZmax,repLenX,zLen,xLen,zExtent,xExtent,lpfilt)
% calculate the mass of each edge and assign it to nodes located at the center of each edge
% this is because edges can be 1, sqrt(2), or sqrt(3) units long in 3d
[density,nodes] = segmentLengths(nodes,edges);
% shift the nodes so that the xy center of mass is at the origin
allXYpos = [nodes(:,1)-density'*nodes(:,1)/sum(density) nodes(:,2)-density'*nodes(:,2)/sum(density)];
% rotate the nodes so that the xy component of the longest principal axis coincides with the main diagonal - to save space and computation
allXYpos = align2dDataWithMainDiagonal([allXYpos nodes(:,3)-density'*nodes(:,3)/sum(density)],density);
% normalize nodes by the SAC surface unit
allZpos = (2*nodes(:,3) - medVZmin)/(medVZmax - medVZmin);
% prepare for gridding by scaling the data so that data of interest lies within [-0.5,0.5]
allZpos = allZpos/zExtent; allXYpos = allXYpos/xExtent;
% grid the data and perform filtering
dist3d = gridder_KBinZ_NNinXY(allZpos,allXYpos(:,1),allXYpos(:,2),density,zLen,xLen,repLenX,lpfilt);
% normalize the arbor density representation so that the norm is equal to the total length of the arbor
dist3d = sqrt(sum(density))*dist3d/sqrt(dist3d(:)'*dist3d(:));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xy = align2dDataWithMainDiagonal(xyz, weights)
if nargin < 2; weights = ones(size(xyz,1),1); end;
% find the inertia tensor
inertiaTensor = zeros(3);
inertiaTensor(1,1) = sum(weights .* (xyz(:,2).^2 + xyz(:,3).^2)); inertiaTensor(2,2) = sum(weights .* (xyz(:,1).^2 + xyz(:,3).^2));
inertiaTensor(3,3) = sum(weights .* (xyz(:,1).^2 + xyz(:,2).^2)); inertiaTensor(1,2) = -sum(weights .* xyz(:,1) .* xyz(:,2));
inertiaTensor(1,3) = -sum(weights .* xyz(:,1) .* xyz(:,3)); inertiaTensor(2,3) = -sum(weights .* xyz(:,2) .* xyz(:,3));
inertiaTensor(2,1) = inertiaTensor(1,2); inertiaTensor(3,1) = inertiaTensor(1,3); inertiaTensor(3,2) = inertiaTensor(2,3);
% find the principal axes of the inertia tensor
[principalAxes, evMatrix] = eig(inertiaTensor);
% take the projection of the 1st principle axis onto the xy plane
pA = principalAxes(1:2,1); pA = pA/norm(pA);
% find the rotation to align pA with the x-axis
pA = pA * sign(xyz(1,1:2)*pA);
%pA = pA * sign(weights'*xyz(:,1:2)*pA);
rotAngle1 = -sign(pA(2))*acos(pA(1)); % (negative of) angle with x-axis
rotMatrix1 = cos(rotAngle1)*eye(2) + sin(rotAngle1)*[0 -1;1 0]; rotMatrix1 = sqrt(1/2)*[1 1;-1 1]*rotMatrix1;
xy = xyz(:,1:2)*rotMatrix1';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [localMass,newNodes] = segmentLengths(nodes,edges)
% assign new nodes at the center of mass of each edge and calculate the mass (length) of each edge
localMass = zeros(size(nodes,1),1); newNodes = zeros(size(nodes,1),3);
for kk=1:size(nodes,1);
parent = edges(find(edges(:,1)==kk),2);
if ~isempty(parent)
localMass(kk) = norm(nodes(parent,:)-nodes(kk,:)); newNodes(kk,:) = (nodes(parent,:)+nodes(kk,:))/2;
else
localMass(kk) = 0; newNodes(kk,:) = nodes(kk,:);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function interpolated = gridder_KBinZ_NNinXY(zSamples,xSamples,ySamples,density,n,m,repLenX,lpfilt)
% data must be in [-0.5, 0.5]
alphaZ=2; Wz=3; error=1e-3; Sz=ceil(0.91/error/alphaZ); beta=pi*sqrt((Wz/alphaZ*(alphaZ-1/2))^2-0.8);
Gz=alphaZ*n; F_kbZ=besseli(0,beta*sqrt(1-([-1:2/(Sz*Wz):1]).^2)); z=-(alphaZ*n/2):(alphaZ*n/2)-1; F_kbZ=F_kbZ/max(F_kbZ(:));
kbZ=sqrt(alphaZ*n)*sin(sqrt((pi*Wz*z/Gz).^2-beta^2))./sqrt((pi*Wz*z/Gz).^2-beta^2); % generate Fourier transform of 1d interpolating kernels
% zero out output array in the alpha grid - use a vector representation to be able to use sparse matrix structure
n = alphaZ*n; interpolated = sparse(n*m*m,1);
% convert samples to matrix indices
nz = (n/2+1) + n*zSamples;
% nearest neighbor interpolation in XY results in some frequency aliasing.
% Low-pass filtering to obtain an arbor density estimate will remove aliased frequencies anyway
nxt = min(m,max(1,round((m/2+1)+m*xSamples))); nyt = min(m,max(1,round((m/2+1)+m*ySamples)));
% loop over samples in kernel
for lz = -(Wz-1)/2:(Wz-1)/2,
nzt = round(nz+lz); zpos=Sz*((nz-nzt)-(-Wz/2))+1; kwz=F_kbZ(round(zpos))'; nzt = min(max(nzt,1),n);
linearIndices = sub2ind([n m m],nzt,nxt,nyt); %nzt+(nxt-1)*n+(nyt-1)*n*m; % compute linear indices
interpolated = interpolated+sparse(linearIndices,1,density.*kwz,n*m*m,1); % use sparse matrix to turn k-space trajectory into 2D matrix
end
interpolated=reshape(full(interpolated),n,m,m); interpolated([1 n],:,:)=0; % edges may be due to samples outside the matrix
n = n/alphaZ; interpolated = myifft3(interpolated,n,repLenX,repLenX); % pre-filter for decimation, low-pass filtering for arbor density
kbZtmp = kbZ((alphaZ-1)*n/2+1:(alphaZ+1)*n/2); deapod = repmat(kbZtmp',[1 repLenX repLenX]); interpolated = interpolated./deapod; % deapodize
% further low-pass filtering. also removes ringing
interpolated=abs(myfft3(interpolated.*lpfilt));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F=myfft3(f,u,v,w)
% does the centering(shifts) and normalization
% needs a slight modification to actually support rectangular images
if nargin<4
F=fftshift(fftn(fftshift(f)))/sqrt(prod(size(f)));
else
F=zeros(u,v,w);
zoffset=(u-size(f,1))/2; xoffset=(v-size(f,2))/2; yoffset=(w-size(f,3))/2;
F(zoffset+1:zoffset+size(f,1),xoffset+1:xoffset+size(f,2),yoffset+1:yoffset+size(f,3))=f;
F=fftshift(fftn(fftshift(F)))/sqrt(prod(size(f)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f=myifft3(F,u,v,w)
% does the centering(shifts) and normalization
% needs a slight modification to actually support rectangular images
if nargin<4
f=ifftshift(ifftn(ifftshift(F)))*sqrt(prod(size(F)));
else
f=ifftshift(ifftn(ifftshift(F)))*sqrt(u*v*w);
zoffset=ceil((size(f,1)-u)/2); xoffset=ceil((size(f,2)-v)/2); yoffset=ceil((size(f,3)-w)/2);
f=f(zoffset+1:zoffset+u,xoffset+1:xoffset+v,yoffset+1:yoffset+w);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function saveZProfile(medVzmin,medVzmax,density,nodes,zRes,zExt,gridPointCount,saveName)
% express z-positions in the universal z-coordinate
allZpos = (nodes(:,3)/zRes - medVzmin)/(medVzmax - medVzmin);
% divide z-positions by the total z-extent so the region of interest is within [-0.5, 0.5]
allZpos = allZpos/zExt;
% grid the 1d z-position data onto a 120-point grid (grid Resolution= zExt*SACdistance/gridPointCount)
zDist = gridder1d(allZpos,density,gridPointCount);
% normalize the profile so that the sum of the profile equals total dendritic length
zDist = zDist * (sum(density)/sum(zDist)/zRes); % assuming each z-bin is 0.5um, sum of profile equals total dendritic length
bins = [-zExt/2:zExt/gridPointCount:zExt/2-zExt/gridPointCount];
% plot the profile
figure;
Xup = gridPointCount*zRes/2; Xlow = -Xup;
h=plot(bins*(gridPointCount*zRes/zExt),zDist); set(h,'LineWidth',2); xlim([[Xlow Xup]]);
ylabel('Dendritic length distribution','FontSize',14); xlabel(strcat('depth in IPL (\mu','m)'),'FontSize',14);
h = gcf; pixelsPerInch = 70; set(h,'Color','w','PaperUnits', 'inches', 'PaperPosition', [0 0 700 400]/pixelsPerInch); set(gca,'box','off');
% save the profile and close the plot;
print(gcf,'-dpng',sprintf('-r%d',pixelsPerInch), strcat(saveName,'_zProfile.png')); close;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function interpolated = gridder1d(zSamples,density,n,Wz,betaZ)
% data must be in [-0.5, 0.5]
alpha=2; W=5; error=1e-3; S=ceil(0.91/error/alpha); beta=pi*sqrt((W/alpha*(alpha-1/2))^2-0.8);
% generate interpolating low-pass filter
Gz=alpha*n; F_kbZ=besseli(0,beta*sqrt(1-([-1:2/(S*W):1]).^2)); z=-(alpha*n/2):(alpha*n/2)-1; %F_kbZ=(Gz/W)*besseli(0,beta*sqrt(1-([-1:2/(S*W):1]).^2));
F_kbZ=F_kbZ/max(F_kbZ(:));
% generate Fourier transform of 1d interpolating kernels
kbZ=sqrt(alpha*n)*sin(sqrt((pi*W*z/Gz).^2-beta^2))./sqrt((pi*W*z/Gz).^2-beta^2);
% zero out output array in the alphaX grid - use a vector representation to be able to use sparse matrix structure
n = alpha*n; interpolated = zeros(n,1);
% convert samples to matrix indices
nz = (n/2+1) + n*zSamples;
% loop over samples in kernel
for lz = -(W-1)/2:(W-1)/2,
% find nearest samples
nzt = round(nz+lz);
% compute weighting for KB kernel using linear interpolation
zpos=S*((nz-nzt)-(-W/2))+1; kwz = F_kbZ(round(zpos))';
% map samples outside the matrix to the edges
nzt = max(nzt,1); nzt = min(nzt,n);
% use sparse matrix to turn k-space trajectory into 2D matrix
interpolated = interpolated+sparse(nzt,1,density.*kwz,n,1);
end;
interpolated = full(interpolated);
% zero out edge samples, since these may be due to samples outside the matrix
interpolated(1) = 0; interpolated(n) = 0;
n = n/alpha;
interpolated = myifft(interpolated,n);
% generate Fourier domain deapodization function and deapodize
deapod = kbZ(n/2+1:3*n/2)'; interpolated = interpolated./deapod;
if nargin>3
z=-(n/2):(n/2)-1;kbZ=sin(sqrt((pi*Wz*z/n).^2-betaZ^2))./sqrt((pi*Wz*z/n).^2-betaZ^2); kbZ(isnan(kbZ))=0; kbZ=kbZ/max(kbZ);
interpolated = abs(myfft3(interpolated.*kbZ'));
else
interpolated = abs(myfft3(interpolated));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f=myifft(F,u)
% does the centering(shifts) and normalization
% needs a slight modification to actually support rectangular images
if nargin<2
f=ifftshift(ifft(ifftshift(F)))*sqrt(numel(F));
else
f=ifftshift(ifft(ifftshift(F)))*sqrt(u);
zoffset=ceil((size(f,1)-u)/2);
f=f(zoffset+1:zoffset+u);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function labelDendrogram(thisLinkage,labels,colorLabels,saveName)
figure;h = dendrogram(thisLinkage,0);
% cell IDs in the dendrogram leaves
curLabels=get(gca,'XTickLabel'); for kk = 1:size(curLabels,1); numCurLabels(kk) = str2num(curLabels(kk,:)); end;
% label leaves
for kk = 1:size(curLabels,1);
x = kk; y = -0.01; t=text(x,y,labels{numCurLabels(kk)},'Color',colorLabels{numCurLabels(kk)});
set(t,'HorizontalAlignment','right','VerticalAlignment','top'); %,'Rotation',90);
end
%Prettier
set(gca,'XLim',[0,size(thisLinkage,1)+2],'YLim',[0,1.05], 'XTickLabel', {}, 'XTick',[], 'YTickLabel', {'0','0.2','0.4','0.6','0.8','1'}, 'YTick',[0:0.2:1]);
for kk = 1:numel(h); set(h(kk),'Color','black','LineWidth',2); end;
pixelsPerInch = 70; set(gcf,'Color','w','PaperUnits', 'inches', 'PaperPosition',[100 100 4000 800]/pixelsPerInch); set(gca,'box','off');
print(gcf,'-dpng',sprintf('-r%d',pixelsPerInch), strcat(saveName,'_dendrogram.png')); close;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function labels=readLabels
% cell types of the 363 cells of the dataset as presented in the paper
labels=cell(364,1); labels{364}='N';
colorLabels=cell(364,1); colorLabels{364}='r';
AA=[84:93 127 304 355]; BB=[94:99 142 153 170 186 201 209 240 246 259 264 285 290 316 318 321 328 354];
CC=[100:105 125 126 130 136 137 149 150 151 172 173 195 205 206 214 215 228 235 237 238 293 299 305 315 319 342 344 345 346 347 348 356];
DD=[65:83 112 119 131 133 145 154 156 157 162 163 164 168 169 171 177 178 181 185 190 191 194 197 198 202 210 218];
DD=[DD 220 221 239 241 245 247 250 251 253 255 268 272 274 286 287 296 306 307 309 310 327 329 343 350 352];
EE=[42:57 108:111 216]; FF=[33:41 114 115 159 160 176 203 225 242 265 276 300 311 330];
GG=[58:64 263 337 363];
HH=[1:32 200 323];
II=[106:107 117 121 122 123 128 141 158 166 179 219 222 226 231 243 244 266 269 288 292 298];
VV=[193 224 283 336 358 359 362];
UU=[129 146 161 167 175 187 189 192 196 204 212 213 229 232 254 281 282 308 312 314 320 322 331 334 349 360];
YY=[139 140 147 148 174 217 227 230 234 252 267 273 277 289 294 295 297 340];
WW=[258 270 280 303 325 341 353 357];
XX=[113 116 132 134 138 152 155 165 180 182 183 184 188 207 233 236 248 256 257 271 275 278 284 302 313 324 326 333 338 339];
ZZ=[118 120 124 135 143 144 199 208 211 223 249 260 261 262 279 291 301 317 332 335 351 361];
labels(AA)={'a'}; labels(BB)={'b'}; labels(CC)={'c'}; labels(DD)={'d'}; labels(EE)={'e'}; labels(FF)={'f'}; labels(GG)={'g'}; labels(HH)={'h'};
labels(II)={'i'}; labels(UU)={'u'}; labels(VV)={'v'}; labels(WW)={'w'}; labels(XX)={'x'}; labels(YY)={'y'}; labels(ZZ)={'z'};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function hardClusterStats = cutDendrogram(myLinkage,clusterIDsForKnownCells,knownCellPositions,maxCutLevel)
% clusterIDsForKnownCells: array of integers of length numel(knowmCellPositions), where cells with same IDs belong to the same cluster
if nargin < 4; maxCutLevel = 0.1; end;
% normalize linkage values
myLinkage(:,3) = myLinkage(:,3)/myLinkage(end,3);
% cutting levels for the linkage
cutStep = 0.001; threshold = cutStep:cutStep:maxCutLevel;
% results of all cuts
randIndices = zeros(size(threshold)); adjustedRandIndices = zeros(size(threshold)); typeConfs = zeros(size(threshold)); widths = zeros(size(threshold));
rowSplits = zeros(size(threshold)); columnSplits = zeros(size(threshold));
% initialize hard-clustering variable for all cutting levels (allT)
allT=ones(size(myLinkage,1)+1,numel(threshold));
for kk = 1:numel(threshold)
% cut the dendrogram at threshold(kk)
T = cluster(myLinkage,'cutoff',threshold(kk),'criterion','distance');
% calculate clustering accuracy metrics
allT(:,kk) = T; [AR,RI,MI,HI,rS,cS,typeConfusions]=reportConfusionsAndRI(T(knownCellPositions),clusterIDsForKnownCells);
randIndices(kk) = RI; adjustedRandIndices(kk) = AR; rowSplits(kk) = rS; columnSplits(kk) = cS; typeConfs(kk) = typeConfusions;
end
% calculate the width of the range of cutting levels yielding the exact same clustering, for each cutting level
for kk = 1:numel(threshold)
flag=true;width=1;while flag&&(kk-width>0);[AR,RI1,MI,HI,~,~,~]=reportConfusionsAndRI(allT(:,kk-width),allT(:,kk));
if RI1==1;width=width+1;else;flag=false;end;end;
flag=true;tmpwidth=0;while flag&&(kk+tmpwidth+1<=numel(threshold));[AR,RI2,MI,HI,~,~,~]=reportConfusionsAndRI(allT(:,kk+tmpwidth+1),allT(:,kk));
if RI2==1;tmpwidth=tmpwidth+1;else;flag=false;end;end;
widths(kk) = width+tmpwidth;
end
[mini0, pos0] = min(typeConfs); allMin = find(typeConfs==mini0); % among minimum type-confusions ...
[~, postemp] = max(randIndices(allMin)); pos0=allMin(postemp); th0=threshold(pos0); width=widths(pos0);
myRowSplits = rowSplits(pos0); myColumnSplits = columnSplits(pos0);
T = cluster(myLinkage,'cutoff',th0,'criterion','distance');
hardClusterStats.minTypeConfs = mini0;
hardClusterStats.structuralSplits = myRowSplits;
hardClusterStats.geneticSplits = myColumnSplits;
hardClusterStats.minCutForMinTypeConfs = th0;
hardClusterStats.cutWidthForMinTypeConfs = width*cutStep;
hardClusterStats.ClusterCountAtMinCut = numel(unique(T));
hardClusterStats.hardClusters = T;
hardClusterStats.riAtMinCut = randIndices(pos0);
hardClusterStats.ariAtMinCut = adjustedRandIndices(pos0);
hardClusterStats.ri = randIndices;
hardClusterStats.ari = adjustedRandIndices;
hardClusterStats.allTypeConfs = typeConfs;
hardClusterStats.allCutWidths = widths*cutStep;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [AR,RI,MI,HI,rowSplits,columnSplits,typeConfusions] = reportConfusionsAndRI(c1,c2)
% This code is slightly modified from the original by Uygar Sümbül.
% The copyright information for the original file is generated below.
% The original file can be downloaded from the following website as of Dec. 23, 2013:
% http://www.mathworks.com/matlabcentral/fileexchange/13916-simple-tool-for-estimating-the-number-of-clusters/content/valid_RandIndex.m
%RANDINDEX - calculates Rand Indices to compare two partitions
% ARI=RANDINDEX(c1,c2), where c1,c2 are vectors listing the
% class membership, returns the "Hubert & Arabie adjusted Rand index".
% [AR,RI,MI,HI]=RANDINDEX(c1,c2) returns the adjusted Rand index,
% the unadjusted Rand index, "Mirkin's" index and "Hubert's" index.
%
% See L. Hubert and P. Arabie (1985) "Comparing Partitions" Journal of
% Classification 2:193-218
%(C) David Corney (2000) D.Corney@cs.ucl.ac.uk
if nargin < 2 | min(size(c1)) > 1 | min(size(c2)) > 1
error('RandIndex: Requires two vector arguments')
return
end
C=Contingency(c1,c2); %form contingency matrix
n=sum(sum(C));
nis=sum(sum(C,2).^2); %sum of squares of sums of rows
njs=sum(sum(C,1).^2); %sum of squares of sums of columns
t1=nchoosek(n,2); %total number of pairs of entities
t2=sum(sum(C.^2)); %sum over rows & columnns of nij^2
t3=.5*(nis+njs);
%Expected index (for adjustment)
nc=(n*(n^2+1)-(n+1)*nis-(n+1)*njs+2*(nis*njs)/n)/(2*(n-1));
A=t1+t2-t3; %no. agreements
D= -t2+t3; %no. disagreements
if t1==nc
AR=0; %avoid division by zero; if k=1, define Rand = 0
else
AR=(A-nc)/(t1-nc); %adjusted Rand - Hubert & Arabie 1985
end
RI=A/t1; %Rand 1971 %Probability of agreement
MI=D/t1; %Mirkin 1970 %p(disagreement)
HI=(A-D)/t1; %Hubert 1977 %p(agree)-p(disagree)
tmp1=sum(C>0,1); tmp2=sum(C>0,2);
rowSplits = sum(tmp1)-nnz(tmp1); columnSplits = sum(tmp2)-nnz(tmp2);
typeConfusions = rowSplits+columnSplits; % number of splits in both directions!
%for kk=1:size(C,2)
% [maxi,pos]=max(C(:,kk)); if pos>kk; tmp = C(pos,:); C(pos,:) = C(kk,:); C(kk,:) = tmp; end;
%end
%tmp = diag(diag(C)); tmp = [tmp; zeros(size(C,1)-size(tmp,1),size(tmp,2))]; tmp = [tmp zeros(size(tmp,1), size(C,2)-size(tmp,2))]; tmp = C-tmp;
%typeConfusions = nnz(tmp);
function Cont=Contingency(Mem1,Mem2)
if nargin < 2 | min(size(Mem1)) > 1 | min(size(Mem2)) > 1
error('Contingency: Requires two vector arguments')
return
end
Cont=zeros(max(Mem1),max(Mem2));
for i = 1:length(Mem1);
Cont(Mem1(i),Mem2(i))=Cont(Mem1(i),Mem2(i))+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function L = elinkage(y, alpha)
%ELINKAGE Minimum energy clustering for Matlab
%
% y data matrix with observations in rows,
% or distances produced by pdist
% alpha power of Euclidean distance, 0<alpha<=2
% default alpha=1
%
% Agglomerative hierarchical clustering by minimum
% energy method, using L-W recursion, for Matlab 7.0.
% See the contributed package "energy" for R, and its
% reference manual, for details:
% http://cran.us.r-project.org/src/contrib/PACKAGES.html
% http://cran.us.r-project.org/doc/packages/energy.pdf
% Reference:
% Szekely, G.J. and Rizzo, M.L. (2005), Hierarchical
% clustering via joint between-within distances:
% extending Ward's minimum variance method, Journal
% of Classification, Vol. 22 (2).
% Email:
% gabors @ bgnet.bgsu.edu, mrizzo @ bgnet.bgsu.edu
% Software developed and maintained by:
% Maria Rizzo, mrizzo @ bgnet.bgsu.edu
% Matlab version 1.0.0 created: 12-Dec-2005
% License: GPL 2.0 or later
%%% initialization
if nargin < 2
alpha = 1.0;
end;
%%% if n==1, y is distance, output from pdist()
[n, d] = size(y); %n obs. in rows
if n == 1 %distances in y
ed = 2 .* squareform(y) .^alpha;
else %data matrix
ed = 2 .* (squareform(pdist(y)).^alpha);
end;
n = length(ed);
clsizes = ones(1, n); %cluster sizes
clindex = 1:n; %cluster indices
L = zeros(n-1, 3); %linkage return value
nclus = n;
for merges = 1:(n-2);
%find clusters at minimum e-distance
b = ones(nclus, 1) .* (ed(1,2) + 1);
B = ed + diag(b);
[colmins, mini] = min(B);
[minmin, j] = min(colmins);
i = mini(j);
%cluster j will be merged into i
%update the linkage matrix
L(merges, 1) = clindex(i);
L(merges, 2) = clindex(j);
L(merges, 3) = ed(i, j);
%update the cluster distances
m1 = clsizes(i);
m2 = clsizes(j);
m12 = m1 + m2;
for k = 1:nclus;
if (k ~= i && k ~= j);
m3 = clsizes(k);
m = m12 + m3;
ed(i, k) = ((m1+m3)*ed(i,k) + (m2+m3)*ed(j,k) - m3*ed(i,j))/m;
ed(k, i) = ed(i, k);
end;
end;
%update cluster data, merge j into i and delete j
ed(j, :) = [];
ed(:, j) = [];
clsizes(i) = m12;
clsizes(j) = [];
clindex(i) = n + merges;
clindex(j) = [];
nclus = n - merges;
%order the leaves so that ht incr left to right
if L(merges, 1) > L(merges, 2);
ij = L(merges, 1);
L(merges, 1) = L(merges, 2);
L(merges, 2) = ij;
end;
end;
%handle the final merge
L(n-1, :) = [clindex(1), clindex(2), ed(1, 2)];
if L(n-1, 1) > L(n-1, 2);
ij = L(n-1, 1);
L(n-1, 1) = L(n-1, 2);
L(n-1, 2) = ij;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function L = elinkage(y, alpha)
%ELINKAGE Minimum energy clustering for Matlab
%
% y data matrix with observations in rows,
% or distances produced by pdist
% alpha power of Euclidean distance, 0<alpha<=2
% default alpha=1
%
% Agglomerative hierarchical clustering by minimum
% energy method, using L-W recursion, for Matlab 7.0.
% See the contributed package "energy" for R, and its
% reference manual, for details:
% http://cran.us.r-project.org/src/contrib/PACKAGES.html
% http://cran.us.r-project.org/doc/packages/energy.pdf
% Reference:
% Szekely, G.J. and Rizzo, M.L. (2005), Hierarchical
% clustering via joint between-within distances:
% extending Ward's minimum variance method, Journal
% of Classification, Vol. 22 (2).
% Email:
% gabors @ bgnet.bgsu.edu, mrizzo @ bgnet.bgsu.edu
% Software developed and maintained by:
% Maria Rizzo, mrizzo @ bgnet.bgsu.edu
% Matlab version 1.0.0 created: 12-Dec-2005
% License: GPL 2.0 or later
%%% initialization
if nargin < 2
alpha = 1.0;
end;
%%% if n==1, y is distance, output from pdist()
[n, d] = size(y); %n obs. in rows
if n == 1 %distances in y
ed = 2 .* squareform(y) .^alpha;
else %data matrix
ed = 2 .* (squareform(pdist(y)).^alpha);
end;
n = length(ed);
clsizes = ones(1, n); %cluster sizes
clindex = 1:n; %cluster indices
L = zeros(n-1, 3); %linkage return value
nclus = n;
for merges = 1:(n-2);
%find clusters at minimum e-distance
b = ones(nclus, 1) .* (ed(1,2) + 1);
B = ed + diag(b);
[colmins, mini] = min(B);
[minmin, j] = min(colmins);
i = mini(j);
%cluster j will be merged into i
%update the linkage matrix
L(merges, 1) = clindex(i);
L(merges, 2) = clindex(j);
L(merges, 3) = ed(i, j);
%update the cluster distances
m1 = clsizes(i);
m2 = clsizes(j);
m12 = m1 + m2;
for k = 1:nclus;
if (k ~= i && k ~= j);
m3 = clsizes(k);
m = m12 + m3;
ed(i, k) = ((m1+m3)*ed(i,k) + (m2+m3)*ed(j,k) - m3*ed(i,j))/m;
ed(k, i) = ed(i, k);
end;
end;
%update cluster data, merge j into i and delete j
ed(j, :) = [];
ed(:, j) = [];
clsizes(i) = m12;
clsizes(j) = [];
clindex(i) = n + merges;
clindex(j) = [];
nclus = n - merges;
%order the leaves so that ht incr left to right
if L(merges, 1) > L(merges, 2);
ij = L(merges, 1);
L(merges, 1) = L(merges, 2);
L(merges, 2) = ij;
end;
end;
%handle the final merge
L(n-1, :) = [clindex(1), clindex(2), ed(1, 2)];
if L(n-1, 1) > L(n-1, 2);
ij = L(n-1, 1);
L(n-1, 1) = L(n-1, 2);
L(n-1, 2) = ij;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [zgrid,xgrid,ygrid] = gridfit(x,y,z,xnodes,ynodes,varargin)
%--------------------------%
%Copyright (c) 2006, John D'Errico
%All rights reserved.
%Redistribution and use in source and binary forms, with or without
%modification, are permitted provided that the following conditions are
%met:
% * Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the distribution
%THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%POSSIBILITY OF SUCH DAMAGE.
%--------------------------%
% gridfit: estimates a surface on a 2d grid, based on scattered data
% Replicates are allowed. All methods extrapolate to the grid
% boundaries. Gridfit uses a modified ridge estimator to
% generate the surface, where the bias is toward smoothness.
%
% Gridfit is not an interpolant. Its goal is a smooth surface
% that approximates your data, but allows you to control the
% amount of smoothing.
%
% usage #1: zgrid = gridfit(x,y,z,xnodes,ynodes);
% usage #2: [zgrid,xgrid,ygrid] = gridfit(x,y,z,xnodes,ynodes);
% usage #3: zgrid = gridfit(x,y,z,xnodes,ynodes,prop,val,prop,val,...);
%
% Arguments: (input)
% x,y,z - vectors of equal lengths, containing arbitrary scattered data
% The only constraint on x and y is they cannot ALL fall on a
% single line in the x-y plane. Replicate points will be treated
% in a least squares sense.
%
% ANY points containing a NaN are ignored in the estimation
%
% xnodes - vector defining the nodes in the grid in the independent
% variable (x). xnodes need not be equally spaced. xnodes
% must completely span the data. If they do not, then the
% 'extend' property is applied, adjusting the first and last
% nodes to be extended as necessary. See below for a complete
% description of the 'extend' property.
%
% If xnodes is a scalar integer, then it specifies the number
% of equally spaced nodes between the min and max of the data.
%
% ynodes - vector defining the nodes in the grid in the independent
% variable (y). ynodes need not be equally spaced.
%
% If ynodes is a scalar integer, then it specifies the number
% of equally spaced nodes between the min and max of the data.
%
% Also see the extend property.
%
% Additional arguments follow in the form of property/value pairs.
% Valid properties are:
% 'smoothness', 'interp', 'regularizer', 'solver', 'maxiter'
% 'extend', 'tilesize', 'overlap'
%
% Any UNAMBIGUOUS shortening (even down to a single letter) is
% valid for property names. All properties have default values,
% chosen (I hope) to give a reasonable result out of the box.
%
% 'smoothness' - scalar or vector of length 2 - determines the
% eventual smoothness of the estimated surface. A larger
% value here means the surface will be smoother. Smoothness
% must be a non-negative real number.
%
% If this parameter is a vector of length 2, then it defines
% the relative smoothing to be associated with the x and y
% variables. This allows the user to apply a different amount
% of smoothing in the x dimension compared to the y dimension.
%
% Note: the problem is normalized in advance so that a
% smoothness of 1 MAY generate reasonable results. If you
% find the result is too smooth, then use a smaller value
% for this parameter. Likewise, bumpy surfaces suggest use
% of a larger value. (Sometimes, use of an iterative solver
% with too small a limit on the maximum number of iterations
% will result in non-convergence.)
%
% DEFAULT: 1
%
%
% 'interp' - character, denotes the interpolation scheme used
% to interpolate the data.
%
% DEFAULT: 'triangle'
%
% 'bilinear' - use bilinear interpolation within the grid
% (also known as tensor product linear interpolation)
%
% 'triangle' - split each cell in the grid into a triangle,
% then linear interpolation inside each triangle
%
% 'nearest' - nearest neighbor interpolation. This will
% rarely be a good choice, but I included it
% as an option for completeness.
%
%
% 'regularizer' - character flag, denotes the regularization
% paradignm to be used. There are currently three options.
%
% DEFAULT: 'gradient'
%
% 'diffusion' or 'laplacian' - uses a finite difference
% approximation to the Laplacian operator (i.e, del^2).
%
% We can think of the surface as a plate, wherein the
% bending rigidity of the plate is specified by the user
% as a number relative to the importance of fidelity to
% the data. A stiffer plate will result in a smoother
% surface overall, but fit the data less well. I've
% modeled a simple plate using the Laplacian, del^2. (A
% projected enhancement is to do a better job with the
% plate equations.)
%
% We can also view the regularizer as a diffusion problem,
% where the relative thermal conductivity is supplied.