Ability to use depthmap information

It is now possible to specify the depth at various points of the
flatmount image, for more accurate reconstruction. This is achieved
via a single-channel TIFF file called depthmap.tiff, and a
corresponding CSV file depthmap.csv, which specifies the background
value.

It is also possible to switch easily between 3D views of flatmount and
reconstruction in the GUI.

Author: davidcsterratt

doc/retistruct-user-guide.tex

@@ -1,6 +1,6 @@
 \documentclass{book}
 
-\newcommand{\svninfo}{Retistruct version: 0.6.2 , Date: 2019-12-13}
+\newcommand{\svninfo}{Retistruct version: 0.7.0 , Date: 2020-08-26}
 \pagestyle{myheadings}
 \markboth{\svninfo}{\svninfo}
 
@@ -245,7 +245,7 @@ \section{Reading files using the \textsc{ImageJ} ROI format}
 \textbf{Outline colour}.
 \end{enumerate}
 
-\subsection{Reading in images comprising multiple fragments}
+\subsection{Reading in images comprising multiple fragments [Optional]}
 \label{retistruct-user-guide:sec:read-imag-compr}
 
 Images from retinae that been cut into multiple fragments can also be
@@ -259,9 +259,38 @@ \subsection{Specifying the scale} Optionally, the scale of the image
 lines, the first with the heading \texttt{Scale} and the second with
 the length of the side of a pixel in micrometres. For example:
 \begin{verbatim}
-"Scale"
+"XY"
 1.5
 \end{verbatim}
+Note that for compatibility with older versions of
+\texttt{Retistruct}, this column can also be called ``Scale''.
+
+\subsection{Specifying depth information [Optional]}
+\label{retistruct-user-guide:sec:specifying-depth}
+
+A ``flat-mount'' needn't be flat. If depth information is available,
+this should be supplied in the form of a single-channel TIFF file
+called \texttt{depthmap.tiff} with the same dimensions as
+\texttt{image.png}. The value of each pixel in the TIFF file specifies
+the height of the corresponding pixel in the image. The size of each
+step in the Z-direction should be specified using an extra column
+``Z'' in \texttt{scale.csv}.
+
+Parts of the TIFF image will correspond to parts of the flat-mount that
+are empty; we refer to these parts as the background. Typically the
+background is represented by a value of 0. There will be problems if
+the outline specified in ImageJ strays onto parts of the TIFF depthmap
+image, since it might appear that there are very sharp jumps in the
+depth. To get around this problem, it is necessary to specify the
+value of the background, so that \textsc{Retistruct} can infer the
+depth of any background regions that the outline strays onto by
+looking at the closest depths that are in the foreground. To specify
+the background value, include a file \texttt{depthmap.csv}, in the
+following format:
+\begin{verbatim}
+"Background value"
+0
+\end{verbatim}
 
 \subsection{Reading in data points}
 \label{retistruct-user-guide:sec:read-data-points}
@@ -702,7 +731,7 @@ \subsection{Transform}
 \section{Exporting the plots}
 \label{retistruct-user-guide:sec:exporting-plots}
 
-The \textbf{Bitmap} and \textbf{PDF} buttons above the flatmount view
+The \textbf{Bitmap} and \textbf{PDF} buttons above the flat-mount view
 and the reconstructed view export the image show to a PNG file format
 or an PDF respectively. To change the width of the image in pixels,
 select \textbf{Properties} and set the
@@ -1028,3 +1057,4 @@ \chapter{The IDT Data format}
 %  LocalWords:  Csefalvay's workspace od ReconstructedOutline getDss
 %  LocalWords:  ReconstructedDataset getSss Xenial Github datacounts
 %  LocalWords:  Wulff Az El uncheck lll Azimuthal Conformal RoiSet
+% LocalWords:  depthmap

pkg/retistruct/DESCRIPTION

@@ -12,7 +12,7 @@ Description: Reconstructs retinae by morphing a flat surface with cuts (a
 Version: 0.7.0
 URL: http://davidcsterratt.github.io/retistruct/
 BugReports: https://github.com/davidcsterratt/retistruct/issues
-Date: 2020-08-23
+Date: 2020-08-26
 Depends:
     R (>= 3.5.0)
 Imports:
@@ -25,7 +25,8 @@ Imports:
     RTriangle (>= 1.6-0.9),
     rgl,
     R.matlab,
-    R6
+    R6,
+    tiff
 Suggests:
     spelling,
     testthat,

pkg/retistruct/NAMESPACE

@@ -60,6 +60,7 @@ export(kr.yhat.cart)
 export(line.line.intersection)
 export(list.datasets)
 export(lvsLplot)
+export(morph.dataset.to.parabola)
 export(normalise.angle)
 export(orthographic)
 export(panlabel)

pkg/retistruct/NEWS

@@ -1,11 +1,17 @@
-CHANGES IN VERSION 0.7.0 - Released 2020/08/23
+CHANGES IN VERSION 0.7.0 - Released 2020/08/26
 
 NEW FEATURES
 
-* Ability to stitch together separate petals or "fragments"
+* Stitching together separate petals or "fragments"
 
 * Ability to magnify plots in the GUI
 
+* Depthmaps, i.e. the ability to specify the depth at various points
+  of the flatmount image, for more accurate reconstruction.
+
+* Switch easily between 3D views of flatmount and reconstruction in
+  GUI
+
 CHANGES IN VERSION 0.6.3 - Released 2020/04/03
 
 BUG FIX

pkg/retistruct/R/Outline.R

@@ -10,8 +10,10 @@ Outline <- R6Class("Outline",
     ##' @field P A N-by-2 matrix of points of the \code{Outline}
     ##'   arranged in anticlockwise order
     P=NULL,
-    ##' @field scale The length of one unit of \code{P} in arbitrary units
-    scale=NULL,        
+    ##' @field scale The length of one unit of \code{P} in the X-Y plane in arbitrary units
+    scale=NULL,
+    ##' @field scalez The length of one unit of \code{P} in the Z-direction in arbitrary units
+    scalez=NULL,
     ##' @field units String giving units of scaled P, e.g. \dQuote{um}
     units=NA,
     ##' @field gf For each row of \code{P}, the index of \code{P} that
@@ -25,19 +27,32 @@ Outline <- R6Class("Outline",
     h=NULL,
     ##' @field im An image as a \code{raster} object
     im=NULL,
+    ##' @field dm Depthmap, with same dimensions as \code{im}, which indicates
+    ##' height of each pixel in Z-direction
+    dm=NULL,
     ##' @field A.fragments Areas of fragments
     A.fragments = NULL,
     ##' @description Construct an outline object. This sanitises the
     ##'   input points \code{P}.
     ##' @param fragments A list of N-by-2 matrix of points for each fragment of the \code{Outline}
     ##' @param scale The length of one unit of \code{P} in arbitrary units
     ##' @param im The image as a \code{raster} object
+    ##' @param scalez The length of one unit of \code{P} in the Z-direction in arbitrary units. If \code{NA}, the depthmap is ignored.
+    ##' @param dm Depthmap, with same dimensions as \code{im}, which indicates
+    ##' height of each pixel in Z-direction
     ##' @param units String giving units of scaled P, e.g. \dQuote{um}
-    initialize=function(fragments=list(), scale=NA, im=NULL, units=NA) {
-      self$P <- matrix(0, 0, 3)
-      colnames(self$P) <- c("X", "Y", "FID")
+    initialize=function(fragments=list(), scale=NA, im=NULL, scalez=NA, dm=NULL, units=NA) {
+      self$P <- matrix(0, 0, 4)
+      colnames(self$P) <- c("X", "Y", "Z", "FID")
+      if (!is.null(dm) & !is.null(im)) {
+        if (all(dim(im) != dim(dm))) {
+          stop("Image and depthmap must have the same dimensions")
+        }
+      }
       self$im <- im
+      self$dm <- dm
       self$scale <- scale
+      self$scalez <- scalez
       self$units <- sub(units, "um", "\U00B5m")
       if (!is.list(fragments)) {
         fragments <- list(fragments)
@@ -60,6 +75,9 @@ Outline <- R6Class("Outline",
     ##' @description Image setter
     ##' @param im An image as a \code{raster} object
     replaceImage = function(im) {
+      if (all(dim(im) != dim(self$dm))) {
+        stop("Image and depthmap must have the same dimensions")
+      }
       self$im <- im
     },
     ##' @description Map the point IDs of a \link{Fragment} on the
@@ -84,22 +102,42 @@ Outline <- R6Class("Outline",
       y[pids[nna]] <- pids[x[nna]]
       return(y)
     },
+    ##' @description Get depth of points P
+    ##' @param P matrix containing unscaled X-Y coordinates of points
+    ##' @return Vector of unscaled Z coordinates of points P
+    getDepth = function(P) {
+      if (!is.null(self$dm)) {
+        Z <- interpolate.image(self$dm, P, invert.y=TRUE)
+      } else {
+        Z <- rep(0, nrow(P))
+      }
+      return(Z)
+    },
+
     ##' @description Add points to the outline register of points
-    ##' @param P 2 column matrix of points to add
-    ##' @param fid fragment id of the points
+    ##' @param P 2 or 3 column matrix of points to add
+    ##' @param fid ID of fragment to which to add  the points
     ##' @return The ID of each added point in the register. If points already
     ##'   exist a point will not be created in the register,
     ##'   but an ID will be returned 
     addPoints = function(P, fid) {
-      if (!is.matrix(P)) {
-        if (length(P) == 2) {
+      if (!(is.vector(P) | is.matrix(P))) {
+        stop("P must be a matrix with 2 or 3 columns, or a vector of length 2 or 3")
+      }
+      if (is.vector(P)) {
+        if (length(P) == 2 | length(P) == 3) {
           P <- matrix(P, nrow=1)
         } else {
-          stop("P must be a matrix or vector of length 2")
+          stop("P should be vector of length 2 or 3")
+        }
+      }
+      if (is.matrix(P)) {
+        if (!(ncol(P) == 2 | ncol(P) == 3)) {
+          stop("P should be a matrix with 2 or 3 columns")
         }
       }
-      if (!(ncol(P) == 2)) {
-        stop("P must have 2 (X, Y) columns")
+      if (ncol(P) == 2) {
+        P <- cbind(P, self$getDepth(P))
       }
       pids <- rep(NA, nrow(P))
       for (i in (1:nrow(P))) {
@@ -109,7 +147,7 @@ Outline <- R6Class("Outline",
           self$h[1] <- 1
         } else {
           ## Check point doesn't already exist
-          id <- which(apply(t(self$P[,-3,drop=FALSE]) == P[i,], 2, all))
+          id <- which(apply(t(self$P[,1:2,drop=FALSE]) == P[i,1:2], 2, all))
           if (length(id) > 1) {
             stop(paste("Point register has duplicates", self$P[id,], collapse=", "))
           }
@@ -139,7 +177,7 @@ Outline <- R6Class("Outline",
     ##' @param fid fragment id of the points
     ##' @return Vector of points
     getFragmentPoints = function(fid) {
-      return(self$P[self$getFragmentPointIDs(fid),c("X", "Y")])
+      return(self$P[self$getFragmentPointIDs(fid),c("X", "Y", "Z")])
     },
     ##' @description Get fragment
     ##' @param fid Fragment ID
@@ -167,18 +205,36 @@ Outline <- R6Class("Outline",
       return(unique(self$P[,"FID"]))
     },
     ##' @description Get unscaled mesh points
+    ##' @param pids IDs of point to return
+    ##' @return  Matrix with columns \code{X}, \code{Y} and \code{Z}
+    getPoints = function(pids=NULL) {
+      if (!is.null(pids)) {
+        return(self$P[pids,c("X", "Y", "Z")])
+      }
+      return(self$P[,c("X", "Y", "Z")])
+    },
+    ##' @description Get X-Y coordinates of unscaled mesh points
+    ##' @param pids IDs of point to return
     ##' @return  Matrix with columns \code{X} and \code{Y}
-    getPoints = function() {
+    getPointsXY = function(pids=NULL) {
+      if (!is.null(pids)) {
+        return(self$P[pids,c("X", "Y")])
+      }
       return(self$P[,c("X", "Y")])
     },
     ##' @description Get scaled mesh points
     ##' @return  Matrix with columns \code{X} and \code{Y} which is
     ##'   exactly \code{scale} times the matrix returned by \code{getPoints}
     getPointsScaled = function() {
       if (is.na(self$scale)) {
-        return(self$P[,c("X", "Y")])
+        return(self$P[,c("X", "Y", "Z")])
+      }
+      if (is.na(self$scalez)) {
+        return(cbind(self$scale*self$P[,c("X", "Y")],
+                     Z=0))
       }
-      return(cbind(self$scale*self$P[,c("X", "Y")]))
+      return(cbind(self$scale*self$P[,c("X", "Y")],
+                   Z=self$scalez*self$P[,"Z"]))
     },
     ##' @description Get set of points on rim
     ##' @return Vector of point IDs, i.e. indices of the rows in

pkg/retistruct/R/PathOutline.R

@@ -56,17 +56,19 @@ PathOutline <- R6Class("PathOutline",
       if (!((self$gf[i0] == i1) || (self$gf[i1] == i0))) {
         stop("Points", i0, "and", i1, "are not connected by an edge")
       }
+      if ((f >= 1) | (f <= 0)) {
+        stop("f argument should be between 0 and 1")
+      }
       fid  <- self$getFragmentIDsFromPointIDs(i0)
       fid1  <- self$getFragmentIDsFromPointIDs(i1)
       if (fid != fid1) {
         stop("Fragment IDs of points differ")
       }
-      PXY <- self$getPoints()
       p <-
-        (1 - f)*PXY[i0,] +
-        f*      PXY[i1,]
+        (1 - f)*self$getPointsXY(i0) +
+        f*      self$getPointsXY(i1)
       ## Find the index of any row of P that matches p
-      n <- anyDuplicated(rbind(PXY, p),
+      n <- anyDuplicated(rbind(self$getPointsXY(c(i0, i1)), p),
                          fromLast=TRUE)
       if (n == 0) {
         ## If the point p doesn't exist
@@ -101,6 +103,10 @@ PathOutline <- R6Class("PathOutline",
       Sf <- path.length(VF0, VF1, self$gf, self$hf, self$getPointsScaled())
       Sb <- path.length(VB0, VB1, self$gb, self$hb, self$getPointsScaled())
 
+      report(paste("\nstitchSubpaths(", VF0, ",", VF1, ",", VB0, ",", VB1, ",", epsilon, ")\n"))
+
+      report(paste("Sf =", Sf, ", Sb =", Sb, "\n"))
+
       ## Initialise forward and backward indicies to first points in each
       ## path
       i0 <- VF0

pkg/retistruct/R/ReconstructedOutline.R

@@ -113,7 +113,7 @@ ReconstructedOutline <- R6Class("ReconstructedOutline",
       ## ol$orderRset()
       self$ol <- ol
       self$phi0 <- ol$phi0
-      self$lambda0 <- ol$lambda0
+     self$lambda0 <- ol$lambda0
 
       report("Merging points...")
       self$mergePointsEdges()
@@ -645,7 +645,8 @@ ReconstructedOutline <- R6Class("ReconstructedOutline",
                    as.vector(outer(ys, xs*0, FUN="+")))
 
         ## Find Barycentric coordinates of corners of pixels
-        Ib <- tsearch(self$ol$getPoints()[,"X"], self$ol$getPoints()[,"Y"],
+
+        Ib <- tsearch(self$ol$getPointsXY()[,"X"], self$ol$getPointsXY()[,"Y"],
                       self$ol$T, I[,1], I[,2], bary=TRUE)
         rm(I)
         gc()
@@ -848,7 +849,7 @@ flatplot.ReconstructedOutline <- function(x, axt="n",
       } else {
         col <- grid.min.col
       }
-      P1 <- get.gridline.flat(x$ol$getPoints(), x$ol$T, x$phi, x$lambda, x$Tt,
+      P1 <- get.gridline.flat(x$ol$getPointsXY(), x$ol$T, x$phi, x$lambda, x$Tt,
                               c(0,0,1), sin(Phi*pi/180))
       cols <- c(cols, rep(col, nrow(P1)))
       P <- rbind(P, P1)
@@ -860,7 +861,7 @@ flatplot.ReconstructedOutline <- function(x, axt="n",
         col <- grid.min.col
       }
       Lambda <- Lambda * pi/180
-      P1 <- get.gridline.flat(x$ol$getPoints(), x$ol$T, x$phi, x$lambda, x$Tt,
+      P1 <- get.gridline.flat(x$ol$getPointsXY(), x$ol$T, x$phi, x$lambda, x$Tt,
                               c(sin(Lambda),cos(Lambda),0), 0)
       cols <- c(cols, rep(col, nrow(P1)))
       P <- rbind(P, P1)

pkg/retistruct/R/StitchedOutline.R

@@ -57,8 +57,8 @@ StitchedOutline <- R6Class("StitchedOutline",
       self$h <- r$h
       for (i in 1:nrow(self$tears)) {
         self$stitchSubpaths(self$tears[i,"V0"], self$tears[i,"VF"],
-          self$tears[i,"V0"], self$tears[i,"VB"],
-          epsilon=self$epsilon)
+                            self$tears[i,"V0"], self$tears[i,"VB"],
+                            epsilon=self$epsilon)
       }
 
       ## Link up points on rim
@@ -85,11 +85,6 @@ StitchedOutline <- R6Class("StitchedOutline",
                    paste(self$Rset, collapse=", ")))
       }
 
-      ## VF0 <- self$corrs[,1]
-      ## VF1 <- self$corrs[,2]
-      ## VB0 <- self$corrs[,4]
-      ## VB1 <- self$corrs[,3]
-
       for (i in 1:nrow(self$corrs)) {
         self$stitchSubpaths(self$corrs[i,"VF0"], self$corrs[i,"VF1"],
                             self$corrs[i,"VB0"], self$corrs[i,"VB1"],

pkg/retistruct/R/TriangulatedOutline.R

@@ -69,7 +69,7 @@ TriangulatedOutline <- R6Class("TriangulatedOutline",
       }
       ## Find areas and lengths of connections
       P <- self$getPointsScaled()
-      self$A <- tri.area(cbind(P, 0), self$T)
+      self$A <- tri.area(P, self$T)
       self$L <- vecnorm(P[self$Cu[,1],] - P[self$Cu[,2],])
       self$A.tot <- sum(self$A)
     },