Error unclear with infinite determinant in retistruct.reconstruct

[bc]

I'm getting an odd error with reconstruction that I haven't come across. What do I need to change to my data to make it work in this context? I attached the rad object rad.rds.zip and the original files Ntae_381.zip ; here's how to load it into your environment:

rad <- readRDS('rad.rds')

> rad
$P
      [,1] [,2]
 [1,]  307  587
 [2,]  353  414
 [3,]  439  560
 [4,]  544  507
 [5,]  592  412
 [6,]  614  352
 [7,]  616  299
 [8,]  616  268
 [9,]  549  239
[10,]  517  216
[11,]  549  159
[12,]  573  131
[13,]  568   99
[14,]  516   77
[15,]  487   78
[16,]  373  139
[17,]  448    7
[18,]  285    0
[19,]  218  196
[20,]  130   36
[21,]   28  137
[22,]   27  294
[23,]  147  330
[24,]   18  387
[25,]   61  499
[26,]  136  455
[27,]   64  530
[28,]  196  568

$h
 [1]  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] 26 27 28

$gf
 [1] 28  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
[26] 25 26 27

$gb
 [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
[26] 27 28  1

$scale
[1] NA

$dataset
[1] "/Users/briancohn/Documents/GitHub/bc/retina/inst/extdata/Ntae_381"

$Ds
$Ds$x
            [,1]      [,2]
  [1,] 303.17647  20.00000
  [2,] 335.82353  20.00000
  [3,] 368.47059  20.00000
  [4,] 401.11765  20.00000
  [5,] 270.52941  52.47059
  [6,] 303.17647  52.47059
  [7,] 335.82353  52.47059
  [8,] 368.47059  52.47059
  [9,] 401.11765  52.47059
 [10,] 237.88235  84.94118
 [11,] 270.52941  84.94118
 [12,] 303.17647  84.94118
 [13,] 335.82353  84.94118
 [14,] 368.47059  84.94118
 [15,]  74.64706 117.41176
 [16,] 107.29412 117.41176
 [17,] 139.94118 117.41176
 [18,] 237.88235 117.41176
 [19,] 270.52941 117.41176
 [20,] 303.17647 117.41176
 [21,] 335.82353 117.41176
 [22,] 368.47059 117.41176
 [23,] 466.41176 117.41176
 [24,] 499.05882 117.41176
 [25,] 531.70588 117.41176
 [26,]  74.64706 149.88235
 [27,] 107.29412 149.88235
 [28,] 139.94118 149.88235
 [29,] 172.58824 149.88235
 [30,] 237.88235 149.88235
 [31,] 270.52941 149.88235
 [32,] 303.17647 149.88235
 [33,] 335.82353 149.88235
 [34,] 368.47059 149.88235
 [35,] 401.11765 149.88235
 [36,] 433.76471 149.88235
 [37,] 466.41176 149.88235
 [38,] 499.05882 149.88235
 [39,] 531.70588 149.88235
 [40,] 107.29412 182.35294
 [41,] 139.94118 182.35294
 [42,] 172.58824 182.35294
 [43,] 237.88235 182.35294
 [44,] 270.52941 182.35294
 [45,] 303.17647 182.35294
 [46,] 335.82353 182.35294
 [47,] 368.47059 182.35294
 [48,] 401.11765 182.35294
 [49,] 433.76471 182.35294
 [50,] 466.41176 182.35294
 [51,] 499.05882 182.35294
 [52,]  74.64706 214.82353
 [53,] 107.29412 214.82353
 [54,] 139.94118 214.82353
 [55,] 172.58824 214.82353
 [56,] 205.23529 214.82353
 [57,] 237.88235 214.82353
 [58,] 270.52941 214.82353
 [59,] 303.17647 214.82353
 [60,] 335.82353 214.82353
 [61,] 368.47059 214.82353
 [62,] 401.11765 214.82353
 [63,] 433.76471 214.82353
 [64,] 466.41176 214.82353
 [65,] 499.05882 214.82353
 [66,]  42.00000 247.29412
 [67,]  74.64706 247.29412
 [68,] 107.29412 247.29412
 [69,] 139.94118 247.29412
 [70,] 172.58824 247.29412
 [71,] 205.23529 247.29412
 [72,] 270.52941 247.29412
 [73,] 303.17647 247.29412
 [74,] 335.82353 247.29412
 [75,] 368.47059 247.29412
 [76,] 401.11765 247.29412
 [77,] 433.76471 247.29412
 [78,] 466.41176 247.29412
 [79,] 499.05882 247.29412
 [80,]  42.00000 279.76471
 [81,]  74.64706 279.76471
 [82,] 107.29412 279.76471
 [83,] 139.94118 279.76471
 [84,] 172.58824 279.76471
 [85,] 205.23529 279.76471
 [86,] 303.17647 279.76471
 [87,] 335.82353 279.76471
 [88,] 368.47059 279.76471
 [89,] 433.76471 279.76471
 [90,] 466.41176 279.76471
 [91,] 531.70588 279.76471
 [92,] 564.35294 279.76471
 [93,] 597.00000 279.76471
 [94,] 139.94118 312.23529
 [95,] 172.58824 312.23529
 [96,] 205.23529 312.23529
 [97,] 270.52941 312.23529
 [98,] 303.17647 312.23529
 [99,] 368.47059 312.23529
[100,] 401.11765 312.23529
[101,] 466.41176 312.23529
[102,] 499.05882 312.23529
[103,] 531.70588 312.23529
[104,] 564.35294 312.23529
[105,] 139.94118 344.70588
[106,] 172.58824 344.70588
[107,] 205.23529 344.70588
[108,] 237.88235 344.70588
[109,] 270.52941 344.70588
[110,] 303.17647 344.70588
[111,] 335.82353 344.70588
[112,] 368.47059 344.70588
[113,] 401.11765 344.70588
[114,] 433.76471 344.70588
[115,] 466.41176 344.70588
[116,] 499.05882 344.70588
[117,] 531.70588 344.70588
[118,] 564.35294 344.70588
[119,] 107.29412 377.17647
[120,] 139.94118 377.17647
[121,] 172.58824 377.17647
[122,] 205.23529 377.17647
[123,] 237.88235 377.17647
[124,] 270.52941 377.17647
[125,] 303.17647 377.17647
[126,] 335.82353 377.17647
[127,] 368.47059 377.17647
[128,] 401.11765 377.17647
[129,] 433.76471 377.17647
[130,] 466.41176 377.17647
[131,] 499.05882 377.17647
[132,] 531.70588 377.17647
[133,]  42.00000 409.64706
[134,]  74.64706 409.64706
[135,] 107.29412 409.64706
[136,] 139.94118 409.64706
[137,] 172.58824 409.64706
[138,] 205.23529 409.64706
[139,] 237.88235 409.64706
[140,] 270.52941 409.64706
[141,] 303.17647 409.64706
[142,] 335.82353 409.64706
[143,] 368.47059 409.64706
[144,] 401.11765 409.64706
[145,] 433.76471 409.64706
[146,] 466.41176 409.64706
[147,] 499.05882 409.64706
[148,] 531.70588 409.64706
[149,]  74.64706 442.11765
[150,] 107.29412 442.11765
[151,] 139.94118 442.11765
[152,] 172.58824 442.11765
[153,] 205.23529 442.11765
[154,] 237.88235 442.11765
[155,] 270.52941 442.11765
[156,] 335.82353 442.11765
[157,] 401.11765 442.11765
[158,] 433.76471 442.11765
[159,] 466.41176 442.11765
[160,] 499.05882 442.11765
[161,] 531.70588 442.11765
[162,]  42.00000 474.58824
[163,] 139.94118 474.58824
[164,] 172.58824 474.58824
[165,] 205.23529 474.58824
[166,] 237.88235 474.58824
[167,] 270.52941 474.58824
[168,] 303.17647 474.58824
[169,] 401.11765 474.58824
[170,] 433.76471 474.58824
[171,] 466.41176 474.58824
[172,] 499.05882 474.58824
[173,] 531.70588 474.58824
[174,] 107.29412 507.05882
[175,] 139.94118 507.05882
[176,] 172.58824 507.05882
[177,] 205.23529 507.05882
[178,] 237.88235 507.05882
[179,] 270.52941 507.05882
[180,] 303.17647 507.05882
[181,] 433.76471 507.05882
[182,] 466.41176 507.05882
[183,] 499.05882 507.05882
[184,] 205.23529 539.52941
[185,] 237.88235 539.52941
[186,] 270.52941 539.52941
[187,] 303.17647 539.52941
[188,] 270.52941 572.00000
[189,] 303.17647 572.00000
[190,] 334.00000 457.00000
[191,] 294.00000 443.00000
[192,] 236.00000 409.00000
[193,] 220.00000 392.00000
[194,] 211.00000 342.00000
[195,] 211.00000 308.00000
[196,] 211.00000 272.00000
[197,] 214.00000 249.00000
[198,] 242.00000 248.00000
[199,] 269.00000 249.00000
[200,] 269.00000 261.00000
[201,] 269.00000 275.00000
[202,] 269.00000 289.00000
[203,] 262.00000 301.00000
[204,] 251.00000 316.00000
[205,] 242.00000 324.00000
[206,] 241.00000 332.00000
[207,] 241.00000 338.00000
[208,] 241.00000 345.00000
[209,] 241.00000 350.00000
[210,] 241.00000 354.00000
[211,] 241.00000 365.00000
[212,] 246.00000 375.00000
[213,] 247.00000 381.00000
[214,] 249.00000 388.00000
[215,] 252.00000 391.00000
[216,] 259.00000 397.00000
[217,] 269.00000 403.00000
[218,] 275.00000 405.00000
[219,] 288.00000 424.00000
[220,] 314.00000 434.00000
[221,] 323.00000 437.00000
[222,] 328.00000 440.00000
[223,] 331.00000 440.00000

$Ss
named list()

$cols
$cols$OD
[1] "blue"

$cols$default
[1] "orange"

$cols$x
[1] "cyan"

$raw
$raw$outline
        x   y
 [1,] 307  17
 [2,] 353 190
 [3,] 439  44
 [4,] 544  97
 [5,] 592 192
 [6,] 614 252
 [7,] 616 305
 [8,] 616 336
 [9,] 549 365
[10,] 517 388
[11,] 549 445
[12,] 573 473
[13,] 568 505
[14,] 516 527
[15,] 487 526
[16,] 373 465
[17,] 448 597
[18,] 285 604
[19,] 218 408
[20,] 130 568
[21,]  28 467
[22,]  27 310
[23,] 147 274
[24,]  18 217
[25,]  61 105
[26,] 136 149
[27,]  64  74
[28,] 196  36

$V0
V0 V0 V0 V0 V0 V0
26  2 10 16 19 23

$VB
VB VB VB VB VB VB
27  3 12 17 20 24

$VF
VF VF VF VF VF VF
25  1  8 14 18 22

$phi0
phi0
   0

$i0
Nasal
   12

$DVflip
DVflip
     0

$side
[1] "Left"

attr(,"class")
[1] "retinalDataset"   "dataset"          "annotatedOutline" "outline"

Browse[1]> retistruct.reconstruct(rad, report = do_not_print, plot.3d = TRUE)

```r
[1] "Triangulating..."
[1] "Stitching..."
[1] "Triangulating..."
[1] "Merging points..."
[1] "Projecting to sphere..."

Error in stretchMesh(Cut, Lt, Rsett, circle(L = L.Rsett)) :
  det(D) is infinite

thanks!

[davidcsterratt]

Sorry about this - it's due to some flakiness in the stitching algorithm. The problem is that for every point along each side or each tear, the stitching algorithm creates a point in the opposing tear unconditionally. This has the unfortunate side-effect that it can create very short edges.

retistruct.reconstruct(rad)
debugger()
# Choose frame 2
range(r$L)
[1] 1.421085e-14 5.998541e+01
flatplot(r)
points(r$P[r$Cu[which.min(r$L),],])

You should be able to see the two points on either side of the edge are indistinguishable.

This huge range then causes problems for the code that stretches out the mesh points as evenly as it can before the final energy minimisation.

I think the solution is to modify the stitching algorithm so that new points are only created when there are no points within a tolerance (say 1% of tear distance) of the corresponding fractional distance on the opposing tear. This change needs a bit of thought, as the stitching algorithm is probably the most fiddly bit of Retistruct code.

[bc]

ah I see. There's a step in retistruct that simplifies the outline. Perhaps that step could help us avoid this issue?

[davidcsterratt]

It might help, but I think there is still a fundamental problem. This is an issue for others too, so I'm just putting my mind to it now.

[davidcsterratt]

I'm working on making the stitching algorithm only add new points when they lie outwith a given tolerance of an existing point. As with everything to do with the stitching algorithm, this is fiddly, so it's taking longer than anticipated.

[davidcsterratt]

This fix is in the master (i.e. development) branch that I have just refactored (see #25). I've tested your example and it seems to work. There are a number of things that aren't working in this refactored code - hopefully I will be able to clear up the mess next week.