Stitch outlines comprising multiple fragments with a hole in the centre

The main work has been to generalis the notion of a rim set to a
boundary set. The boundary set is the part of the outline that is not
stitched. It may comprise only one connected part, which is by
definition the rim, or it may comprise multiple closed paths, the
longest of which is assumed to be rim. Once the rim is identified,
reconstruction can start. There may be future work to do with ensuring
the area of the projection is appropriate, but this will probably have
a small influence on the quailty of the reconstruction.

Author: davidcsterratt

doc/retistruct-user-guide.tex

@@ -1,6 +1,6 @@
 \documentclass{book}
 
-\newcommand{\svninfo}{Retistruct version: 0.7.0 , Date: 2020-08-26}
+\newcommand{\svninfo}{Retistruct version: 0.7.1 , Date: 2020-08-28}
 \pagestyle{myheadings}
 \markboth{\svninfo}{\svninfo}
 
@@ -386,7 +386,7 @@ \section{Editing the retinal mark-up}
   \caption{Adding a correspondence. One correspondence between the two
     upper petals has been added. The red dashed line indicate the ends
     of the correspondence, and the red highlighted edge indicates the
-    edges the correspond.}
+    edges of the correspondence.}
   \label{retistruct-user-guide:fig:correspondence}
 \end{figure}
 

pkg/retistruct/DESCRIPTION

@@ -9,10 +9,10 @@ Description: Reconstructs retinae by morphing a flat surface with cuts (a
     shape). It can estimate the position of a point on the intact adult retina
     to within 8 degrees of arc (3.6% of nasotemporal axis). The coordinates in
     reconstructed retinae can be transformed to visuotopic coordinates.
-Version: 0.7.0
+Version: 0.7.1
 URL: http://davidcsterratt.github.io/retistruct/
 BugReports: https://github.com/davidcsterratt/retistruct/issues
-Date: 2020-08-26
+Date: 2020-08-28
 Depends:
     R (>= 3.5.0)
 Imports:

pkg/retistruct/NEWS

@@ -1,3 +1,8 @@
+CHANGES IN VERSION 0.7.1 - Released 2020/08/28
+
+* Stitching together of fragments that contain a hole in the centre
+  is now possible
+
 CHANGES IN VERSION 0.7.0 - Released 2020/08/26
 
 NEW FEATURES

pkg/retistruct/R/AnnotatedOutline.R

@@ -272,12 +272,30 @@ AnnotatedOutline <- R6Class("AnnotatedOutline",
       return(self$i0)
     },
     ##' @description Get point IDs of points on rim
-    ##' @return Point IDs of points on rim
+    ##' @return Point IDs of points on rim. If the outline has been
+    ##' stitched (see \code{\link{StitchedOutline}}), the point IDs
+    ##' will be ordered in the direction of the
+    ##' forward pointer.
     getRimSet = function() {
+      ## TR <- self$computeTearRelationships(self$tears)
+      ## CR <- self$computeCorrespondenceRelationships(self$corrs)
+      ## Rset <- intersect(TR$Rset, CR$Rset)
+      return(self$getBoundarySets()[["Rim"]])
+    },
+    ##' @description Get point IDs of points on boundaries
+    ##' @return List of Point IDs of points on the boundaries.
+    ##' If the outline has been stitched (see \code{\link{StitchedOutline}}),
+    ##' the point IDs in each
+    ##' element of the list will be ordered in the direction of the
+    ##' forward pointer, and the boundary that is longest will be
+    ##' named as \code{Rim}. If the outline has not been stitched,
+    ##' the list will have one element named \code{Rim}.
+    getBoundarySets = function() {
       TR <- self$computeTearRelationships(self$tears)
       CR <- self$computeCorrespondenceRelationships(self$corrs)
-      Rset <- intersect(TR$Rset, CR$Rset)
-      return(Rset)
+      ## The intersection of these two is the union of the boundary
+      ## sets
+      return(list(Rim=intersect(TR$Rset, CR$Rset)))
     },
     ##' @description Ensure that the fixed point \code{i0} is in the rim, not a tear.
     ##' Alters object in which \code{i0} may have been changed. 
@@ -294,13 +312,15 @@ AnnotatedOutline <- R6Class("AnnotatedOutline",
         }
       }
     },
-    ##' @description Get flat rim length
-    ##' @return The rim length 
-    getFlatRimLength = function() {
-      suppressMessages(r <- self$computeTearRelationships(self$V0, self$VB, self$VF))
-      return(path.length(self$i0, path.next(self$i0, self$gf, r$hf), self$gf, r$hf, self$P) +
-             path.length(self$i0, path.next(self$i0, self$gf, r$hf), self$gb, r$hb, self$P))
-    },
+    ## FIXME: Remove getFlatRimLength? It is not used and I'm not convinced it's correct
+    ##
+    ## ##' @description Get flat rim length
+    ## ##' @return The rim length
+    ## getFlatRimLength = function() {
+    ##   suppressMessages(r <- self$computeTearRelationships(self$V0, self$VB, self$VF))
+    ##   return(path.length(self$i0, path.next(self$i0, self$gf, r$hf), self$gf, r$hf, self$P) +
+    ##          path.length(self$i0, path.next(self$i0, self$gf, r$hf), self$gb, r$hb, self$P))
+    ## },
     ##' @description Label a set of four unlabelled points supposed to refer to a
     ##' correspondence.
     ##' @param pids the vector of point indices

pkg/retistruct/R/ReconstructedOutline.R

@@ -287,7 +287,7 @@ ReconstructedOutline <- R6Class("ReconstructedOutline",
       }
 
       ## Transform the rim set
-      Rset <- order.Rset(self$ol$getRimSet(), self$ol$gf, self$ol$h)
+      Rset <- self$ol$getRimSet()
       Rsett <- unique(ht[Rset])
       i0t <- ht[self$ol$i0]
 
@@ -669,8 +669,8 @@ ReconstructedOutline <- R6Class("ReconstructedOutline",
       }
       return(private$ims)
     },
-    ##' @description Get location of tear coordinates in spherical coordinates
-    ##' @return Location of tear coordinates in spherical coordinates
+    ##' @description Get locations of tears in spherical coordinates
+    ##' @return List containing locations of tears in spherical coordinates
     getTearCoords = function() {
       Tss <- list()
       for (TF in self$ol$TFset) {
@@ -680,8 +680,8 @@ ReconstructedOutline <- R6Class("ReconstructedOutline",
       }
       return(Tss)
     },
-    ##' @description Get location of correspondence coordinates in spherical coordinates
-    ##' @return Location of correspondence coordinates in spherical coordinates
+    ##' @description Get locations of correspondences in spherical coordinates
+    ##' @return List containing locations of correspondences in spherical coordinates
     getCorrespondenceCoords = function() {
       Css <- list()
       for (CF in self$ol$CFset) {
@@ -691,6 +691,22 @@ ReconstructedOutline <- R6Class("ReconstructedOutline",
       }
       return(Css)
     },
+    ##' @description Get location of non-rim boundaries in spherical coordinates
+    ##' @return List containing locations of non-rim boundaries in spherical coordinates
+    getNonRimBoundaryCoords = function() {
+      Bsets <- self$ol$getBoundarySets()
+      if (length(Bsets) <= 1) {
+        return(NULL)
+      }
+      Bss <- list()
+      for (B in Bsets[names(Bsets) != "Rim"]) {
+        ## Convert indices to the spherical frame of reference
+        j <- self$ht[B]
+        Bss <- c(Bss, list(cbind(phi=self$phi[j], lambda=self$lambda[j])))
+      }
+      return(Bss)
+
+    },
     ##' @description Get \link{ReconstructedFeatureSet}
     ##' @param type Base type of \link{FeatureSet} as string.
     ##'   E.g. \code{PointSet} returns a \link{ReconstructedPointSet}
@@ -1145,16 +1161,24 @@ projection.ReconstructedOutline <- function(r,
     graphics::polygon(boundary[,"x"], boundary[,"y"], border="black")
   }
 
-  ## Plot rim in visutopic space
+  ## Plot rim in visuotopic space
   rs <- cbind(phi=r$phi0, lambda=seq(0, 2*pi, len=360))
   rs.rot <- rotate.axis(transform(rs, phi0=r$phi0), axisdir*pi/180)
   ## "Home" position for a cyclops looking ahead
   ## r$axisdir = cbind(phi=0, lambda=0)
-
   lines(projection(rs.rot, lambdalim=lambdalim*pi/180, lines=TRUE,
                    proj.centre=pi/180*proj.centre),
         col=getOption("TF.col"))
 
+  ## Plot non-rim boundary in visuotopic space
+  bss <- r$getNonRimBoundaryCoords()
+  for (bs in bss) {
+    bs.rot <- rotate.axis(transform(bs, phi0=r$phi0), axisdir*pi/180)
+    lines(projection(bs.rot, lambdalim=lambdalim*pi/180, lines=TRUE,
+                     proj.centre=pi/180*proj.centre),
+          col=getOption("TF.col"))
+  }
+
   ## Projection of pole
   if (pole) {
     oa.rot <- rotate.axis(transform(cbind(phi=-pi/2, lambda=0), phi0=r$phi0),

pkg/retistruct/R/StitchedOutline.R

@@ -103,6 +103,48 @@ StitchedOutline <- R6Class("StitchedOutline",
         self$h <- self$h[self$h]
       }
       self$CFset <- r$CFset
+    },
+    ##' @description Test if the outline has been stitched
+    ##' @return Boolean, indicating if the outline has been stitched or not
+    isStitched = function() {
+      return((!(is.null(self$TFset))) &
+               ((nrow(self$corrs) == 0) | !(is.null(self$CFset))))
+    },
+    ##' @description Get point IDs of points on boundaries
+    ##' @return List of Point IDs of points on the boundaries.
+    ##' If the outline has been stitched,
+    ##' the point IDs in each
+    ##' element of the list will be ordered in the direction of the
+    ##' forward pointer, and the boundary that is longest will be
+    ##' named as \code{Rim}. If the outline has not been stitched,
+    ##' the list will have one element named \code{Rim}.
+    getBoundarySets = function() {
+      Bsets <- super$getBoundarySets()
+      if (!self$isStitched()) {
+        return(Bsets)
+      }
+      UBset <- Bsets[["Rim"]]
+      ## Now separate out the Bsets
+      Bsets  <- list()
+      Blengths <- c()
+      P <- self$getPointsScaled()
+      while(length(UBset) > 0) {
+        i  <- UBset[1]
+        j  <- path.next(i, self$gf, self$hf)
+        B12 <- path(i, j, self$gf, self$hf)
+        B21 <- path(j, i, self$gf, self$hf)
+        BL12 <- path.length(i, j,
+                            self$gf, self$hf, P)
+        BL21 <- path.length(j, i,
+                            self$gf, self$hf, P)
+        Bset <- c(B12[-1], B21[-1])
+        Bsets <- c(Bsets, list(Bset))
+        Blengths <- c(Blengths, BL12 + BL21)
+        UBset <- setdiff(UBset, Bset)
+      }
+      Bsets <- name.list(Bsets)
+      names(Bsets)[which.max(Blengths)] <- "Rim"
+      return(Bsets)
     }
   )
 )

pkg/retistruct/R/spheristruct.R

@@ -1,14 +1,3 @@
-## order.Rset(Rset, gf, hf)
-##
-## It is nice to create Rset as an ordered set
-order.Rset <- function(Rset, gf, hf) {
-  ## To to this, join the path from the first two members of the set.
-  R12 <- path(Rset[1], Rset[2], gf, hf)
-  R21 <- path(Rset[2], Rset[1], gf, hf)
-  Rset <- c(R12[-1], R21[-1])
-  return(Rset)
-}
-
 ##' Stretch the mesh in the flat retina to a circular outline
 ##'
 ##' @title Stretch mesh

pkg/retistruct/demo/00Index

@@ -5,3 +5,4 @@ orange               Reconstruct a peeled orange!
 wedge		             Demonstration of wedge coordinates
 smi32                Reconstruct a retina labelled with SMI-32
 parabola             Reconstruct a parabolically transformed retina
+hole                 Reconstruct a retina with a holes

pkg/retistruct/demo/hole.R

@@ -0,0 +1,60 @@
+par(mfcol=c(1, 2))
+
+P <- list(rbind(c(1,1.5),
+                c(1.5,1),
+                c(2,1),
+                c(2.5,2),
+                c(3,1),
+                c(4,1),
+                c(1,4)),
+          rbind(c(-1.5,1),
+                c(-1,1.5),
+                c(-1,4),
+                c(-2,3),
+                c(-2,2),
+                c(-3,2),
+                c(-4,1)),
+          rbind(c(-4,-1),
+                c(-1.5,-1),
+                c(-1,-1.5),
+                c(-1,-4)),
+          rbind(c(1,-1.5),
+                c(1.5,-1),
+                c(2,-1),
+                c(2.5,-2),
+                c(3,-1),
+                c(4,-1),
+                c(1,-4)))
+## Stitched outlines
+a <- StitchedOutline$new(P)
+
+expect_false(a$isStitched())
+
+## Set a fixed point
+## One that is in the rim should be fine
+a$setFixedPoint(6, "Nasal")
+expect_equal(a$i0, c(Nasal=6))
+
+## One that is not in the rim should be moved
+a$addTear(c(11, 12, 13))
+a$addTear(c(3, 4, 5))
+a$addTear(c(21, 22, 23))
+
+## Add correspondences
+a$addCorrespondence(c(2, 6, 20, 24))
+a$addCorrespondence(c(1, 7, 9, 10))
+a$addCorrespondence(c(8, 14, 15, 16))
+a$addCorrespondence(c(17, 18, 19, 25))
+
+r <- ReconstructedOutline$new()
+r$loadOutline(a)
+
+flatplot(a)
+flatplot(r)
+
+r$reconstruct()
+
+flatplot(r)
+projection(r)
+
+par(mfcol=c(1, 1))