20220922_spatstat_cycIF.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ec74adfa",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: spatstat.data\n",
      "\n",
      "Loading required package: spatstat.geom\n",
      "\n",
      "spatstat.geom 3.0-3\n",
      "\n",
      "Loading required package: spatstat.random\n",
      "\n",
      "spatstat.random 3.0-1\n",
      "\n",
      "Loading required package: spatstat.explore\n",
      "\n",
      "Loading required package: nlme\n",
      "\n",
      "spatstat.explore 3.0-5\n",
      "\n",
      "Loading required package: spatstat.model\n",
      "\n",
      "Loading required package: rpart\n",
      "\n",
      "spatstat.model 3.0-2\n",
      "\n",
      "Loading required package: spatstat.linnet\n",
      "\n",
      "spatstat.linnet 3.0-3\n",
      "\n",
      "\n",
      "spatstat 3.0-2 \n",
      "For an introduction to spatstat, type ‘beginner’ \n",
      "\n",
      "\n",
      "\n",
      "Attaching package: ‘dplyr’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:nlme’:\n",
      "\n",
      "    collapse\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:stats’:\n",
      "\n",
      "    filter, lag\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:base’:\n",
      "\n",
      "    intersect, setdiff, setequal, union\n",
      "\n",
      "\n",
      "Warning message:\n",
      "“package ‘imager’ was built under R version 4.2.1”\n",
      "Loading required package: magrittr\n",
      "\n",
      "Registered S3 method overwritten by 'imager':\n",
      "  method      from\n",
      "  plot.imlist     \n",
      "\n",
      "\n",
      "Attaching package: ‘imager’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:magrittr’:\n",
      "\n",
      "    add\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:spatstat.geom’:\n",
      "\n",
      "    as.imlist, is.imlist\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:stats’:\n",
      "\n",
      "    convolve, spectrum\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:graphics’:\n",
      "\n",
      "    frame\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:base’:\n",
      "\n",
      "    save.image\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(spatstat)\n",
    "library(dplyr)\n",
    "#library(ggplot2)\n",
    "library(imager)\n",
    "#ulimit::memory_limit(6000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6168d68d",
   "metadata": {},
   "outputs": [],
   "source": [
    "mydata <- read.csv(\"./data/20220420_JP-TMAs_IMC-TMAs_MIBI_CombinedCelltypes_all.csv\")\n",
    "mydata <- mydata[grep(\"JP-TMA1-1\", mydata$X), ]\n",
    "mydata[,\"DAPI_X\"] = mydata[,\"DAPI_X\"]/0.325\n",
    "mydata[,\"DAPI_Y\"] = mydata[,\"DAPI_Y\"]/0.325\n",
    "tail(mydata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0a4cf28",
   "metadata": {},
   "outputs": [],
   "source": [
    "#simplify to tumor\n",
    "mydata[,'leiden'][mydata[,'leidencelltype5'] == 'epithelial'] <- 'tumor'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36c32551",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mydata[,'leiden'][mydata[,'leiden'] == 'CD8 T cell'] <- 'CD3 T cell'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e777a2d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mydata[,'leiden'][mydata[,'leiden'] == 'CD4 T cell'] <- 'CD3 T cell'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f658fe50",
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_scene = unique(mydata$slide_scene)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c035c9d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_scene[1:48]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e578ef5",
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_scene[49:95]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3298b354",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# K cross\n",
    "ls_cells = c('Quies. str.','FN+ FB','Vim+ FB','Macrophage','endothelial','CD3 T cell','CD4 T cell','CD8 T cell')\n",
    "ls_from = c('CD3 T cell')#c('tumor','CD4 T cell','CD8 T cell')\n",
    "combined_df <- data.frame(r=double(),\n",
    "                    theo=double(),\n",
    "                    iso=double(),\n",
    "                    slide_scene=character(),\n",
    "                    cell=character())\n",
    "ls_scene = unique(mydata$slide_scene)\n",
    "# loop version 1\n",
    "datalist = list()\n",
    "for (s_scene in ls_scene[1:48]) {\n",
    "  print(s_scene)\n",
    "  file = sprintf(\"./TissueMask/TissueMask_%s.png\", s_scene)#'TissueMask_JP-TMA1-1_scene034.png'\n",
    "  img <- load.image(file)\n",
    "  mat = as.matrix(img>0)\n",
    "  win_t = owin(mask=t(mat)) #transpose rotate along diagonal\n",
    "  win_p <- as.polygonal(win_t) \n",
    "  mydataSubset <- mydata[grep(s_scene, mydata$X), ]\n",
    "  mypattern <- ppp(x=(mydataSubset[,\"DAPI_X\"]),y=((mydataSubset[,\"DAPI_Y\"])),window=win_p)\n",
    "  isin<-inside.owin(x=(mydataSubset[,\"DAPI_X\"]),y=((mydataSubset[,\"DAPI_Y\"])),w=win_p)\n",
    "  marks(mypattern) <- factor(mydataSubset[isin,'leiden']) #factor(mydataSubset$leiden)\n",
    "  #break\n",
    "  #plot(mypattern)\n",
    "  startTime <- Sys.time()\n",
    "  for (s_cell in intersect((unique(mydataSubset$leiden)),ls_cells)) {\n",
    "      #print(s_cell)\n",
    "      for (s_center in intersect((unique(mydataSubset$leiden)),ls_from)) {\n",
    "          #print(s_center)\n",
    "          Kcross_cell <- Kcross(X=mypattern, from=s_center, to=s_cell,correction='isotropic',rmax=100/.325,nlarge=3000)\n",
    "          #plot(Kcross_cell)\n",
    "          Kcross_cell['slide_scene'] = s_scene\n",
    "          Kcross_cell['from_cell'] = s_center\n",
    "          Kcross_cell['to_cell'] = s_cell\n",
    "          combined_df <- rbind(combined_df, Kcross_cell)\n",
    "          }\n",
    "      break\n",
    "  }\n",
    "  endTime <- Sys.time()\n",
    "  print(endTime - startTime)\n",
    "  break\n",
    "    }\n",
    "\n",
    "#write.csv(combined_df,\"./data/cycIF_JP-TMA1_first_CD3_Kcross.csv\", row.names = FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efc0ff1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(mypattern)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c1b92ca",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#summary(mypattern)\n",
    "#attr(mypattern,\"rejects\")\n",
    "#isin<-inside.owin(x=(mydataSubset[,\"DAPI_X\"]),y=((mydataSubset[,\"DAPI_Y\"])),w=win_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8887332e",
   "metadata": {},
   "source": [
    "## K function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c84e4c36",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "combined_df <- data.frame(r=double(),\n",
    "                    theo=double(),\n",
    "                    iso=double(),\n",
    "                    slide_scene=character(),\n",
    "                    cell=character())\n",
    "ls_scene = unique(mydata$slide_scene)\n",
    "# loop version 1\n",
    "datalist = list()\n",
    "for (s_scene in ls_scene) {\n",
    "  print(s_scene)\n",
    "  file = sprintf(\"./TissueMask/TissueMask_%s.png\", s_scene)#'TissueMask_JP-TMA1-1_scene034.png'\n",
    "  img <- load.image(file)\n",
    "  mat = as.matrix(img>0)\n",
    "  win_t = owin(mask=t(mat)) #transpose rotate along diagonal\n",
    "  win_p <- as.polygonal(win_t) \n",
    "  mydataSubset <- mydata[grep(s_scene, mydata$X), ]\n",
    "  mypattern4 <- ppp(x=(mydataSubset[,\"DAPI_X\"]),y=((mydataSubset[,\"DAPI_Y\"])),window=win_p)\n",
    "  for (s_cell in c('CD4 T cell','CD8 T cell','CD20 B cell')) {\n",
    "    mydata_cell = mydataSubset[grep(s_cell, mydataSubset$leiden), ]\n",
    "    if (nrow(mydata_cell) > 0) {\n",
    "        mypattern <- ppp(x=(mydata_cell[,\"DAPI_X\"]),y=((mydata_cell[,\"DAPI_Y\"])),window=win_p)\n",
    "        startTime <- Sys.time()\n",
    "        Kest_cell <- Kest(mypattern,correction = 'isotropic',rmax=100/.325,nlarge=3000)\n",
    "        endTime <- Sys.time()\n",
    "        print(endTime - startTime)\n",
    "        Kest_cell['slide_scene'] = s_scene\n",
    "        Kest_cell['cell'] = s_cell\n",
    "        combined_df <- rbind(combined_df, Kest_cell)\n",
    "        }\n",
    "    }\n",
    "    #add a CD3 T cell\n",
    "    s_cell = 'T cell'\n",
    "    mydata_cell = mydataSubset[grep(s_cell, mydataSubset$leiden), ]\n",
    "    if (nrow(mydata_cell) > 0) {\n",
    "        mypattern <- ppp(x=(mydata_cell[,\"DAPI_X\"]),y=((mydata_cell[,\"DAPI_Y\"])),window=win_p)\n",
    "        startTime <- Sys.time()\n",
    "        Kest_cell <- Kest(mypattern,correction = 'isotropic',rmax=100/.325,nlarge=3000)\n",
    "        plot(Kest_cell)\n",
    "        endTime <- Sys.time()\n",
    "        print(endTime - startTime)\n",
    "        Kest_cell['slide_scene'] = s_scene\n",
    "        Kest_cell['cell'] = s_cell\n",
    "        combined_df <- rbind(combined_df, Kest_cell)\n",
    "        }\n",
    "  #break\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2e5b396",
   "metadata": {},
   "outputs": [],
   "source": [
    "#write.csv(combined_df,\"./data/cycIF_ripleys_k.csv\", row.names = FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "201532a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "combined_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4b33e07",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Kest_cell %>% select(r,theo,iso)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb6a45cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot(envelope(mypattern,Kest))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60fbdf26",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot(density(mypattern))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aea9675a",
   "metadata": {},
   "outputs": [],
   "source": [
    "marks(mypattern) <- mydata[,c('CD45','Ki67')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ac433a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(Smooth(mypattern,\"msr\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52b4350e",
   "metadata": {},
   "source": [
    "## demo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1425d242",
   "metadata": {},
   "outputs": [],
   "source": [
    "if(dev.cur() <= 1) {\n",
    "  dd <- getOption(\"device\")\n",
    "  if(is.character(dd)) dd <- get(dd)\n",
    "  dd()\n",
    "}\n",
    "\n",
    "oldpar <- par(ask = interactive() && dev.interactive(orNone=TRUE))\n",
    "oldoptions <- options(warn=-1)\n",
    "\n",
    "fanfare <- function(stuff) {\n",
    "  plot(c(0,1),c(0,1),type=\"n\",axes=FALSE, xlab=\"\", ylab=\"\")\n",
    "  text(0.5,0.5, stuff, cex=2.5)\n",
    "}\n",
    "\n",
    "par(mar=c(1,1,2,1)+0.1)\n",
    "\n",
    "fanfare(\"Spatstat demonstration\")\n",
    "\n",
    "fanfare(\"I. Types of data\")\n",
    "plot(swedishpines, main=\"Point pattern\")\n",
    "\n",
    "plot(demopat, cols=c(\"green\", \"blue\"), main=\"Multitype point pattern\")\n",
    "\n",
    "plot(longleaf, fg=\"blue\", main=\"Marked point pattern\")\n",
    "\n",
    "plot(finpines, main=\"Point pattern with multivariate marks\")\n",
    "\n",
    "a <- psp(runif(20),runif(20),runif(20),runif(20), window=owin())\n",
    "plot(a, main=\"Line segment pattern\")\n",
    "marks(a) <- sample(letters[1:4], 20, replace=TRUE)\n",
    "plot(a, main=\"Multitype line segment pattern\")\n",
    "marks(a) <- runif(20)\n",
    "plot(a, main=\"Marked line segment pattern\", style=\"width\")\n",
    "\n",
    "plot(owin(), main=\"Rectangular window\")\n",
    "plot(letterR, main=\"Polygonal window\")\n",
    "plot(as.mask(letterR), main=\"Binary mask window\")\n",
    "\n",
    "Z <- as.im(function(x,y){ sqrt((x - 1)^2 + (y-1)^2)}, square(2))\n",
    "plot(Z, main=\"Pixel image\")\n",
    "\n",
    "X <- runifpoint(42)\n",
    "plot(dirichlet(X), main=\"Tessellation\")\n",
    "\n",
    "plot(rpoispp3(100), main=\"Three-dimensional point pattern\")\n",
    "\n",
    "plot(simplenet, main=\"Linear network (linnet)\")\n",
    "\n",
    "X <- rpoislpp(20, simplenet)\n",
    "plot(X,\n",
    "     main=\"Point pattern on linear network (lpp)\",\n",
    "     show.window=FALSE)\n",
    "\n",
    "V <- as.linim(function(x,y,seg,tp){x^2-y^2}, L=simplenet)\n",
    "plot(V, main=\"Pixel image on a linear network\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e952fdf0",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "fanfare(\"II. Graphics\")\n",
    "\n",
    "plot(letterR, col=\"green\", border=\"red\", lwd=2, main=\"Polygonal window with colour fill\")\n",
    "plot(letterR, hatch=TRUE, spacing=0.15, angle=30, main=\"Polygonal window with line shading\")\n",
    "plot(letterR, hatch=TRUE, hatchargs=list(texture=8, spacing=0.12),\n",
    "     main=\"Polygonal window with texture fill\")\n",
    "\n",
    "plot(amacrine, chars=c(1,16),\n",
    "     main=\"plot(X, chars = c(1,16))\")\n",
    "plot(amacrine, cols=c(\"red\",\"blue\"), chars=16,\n",
    "     main=\"plot(X, cols=c(\\\"red\\\", \\\"blue\\\"))\")\n",
    "\n",
    "opa <- par(mfrow=c(1,2))\n",
    "plot(longleaf, markscale=0.03, main=\"markscale=0.03\")\n",
    "plot(longleaf, markscale=0.09, main=\"markscale=0.09\")           \n",
    "par(opa)\n",
    "\n",
    "plot(longleaf, pch=21, cex=1,\n",
    "     bg=colourmap(terrain.colors(128), range=c(0,80)),\n",
    "     main=\"colourmap for numeric mark values\")\n",
    "\n",
    "Z <- as.im(function(x,y) { r <- sqrt(x^2+y^2); r * exp(-r) },\n",
    "           owin(c(-5,5),c(-5,5)))\n",
    "plot(Z, main=\"pixel image: image plot\")\n",
    "plot(Z, main=\"pixel image: image plot (heat colours)\", col=heat.colors(256))\n",
    "plot(Z, main=\"pixel image: logarithmic colour map\", \n",
    "     log=TRUE, col=rainbow(128, end=5/6))\n",
    "contour(Z, main=\"pixel image: contour plot\", axes=FALSE)\n",
    "plot(Z, main=\"pixel image: image + contour plot\")\n",
    "contour(Z, add=TRUE)\n",
    "persp(Z, colmap=terrain.colors(128), shade=0.3, phi=30,theta=100,\n",
    "      main=\"pixel image: perspective plot\")\n",
    "\n",
    "ct <- colourmap(rainbow(20), breaks=seq(-1,1,length=21))\n",
    "plot(ct, main=\"Colour map for real numbers\")\n",
    "\n",
    "ca <- colourmap(rainbow(8), inputs=letters[1:8])\n",
    "plot(ca, main=\"Colour map for discrete values\")\n",
    "\n",
    "Z <- as.im(nnfun(runifpoint(8)))\n",
    "plot(Z, main=\"colour image for discrete values\")\n",
    "textureplot(Z, main=\"texture plot for discrete values\")\n",
    "\n",
    "W <- owin(c(1,5),c(0,4.5))\n",
    "Lout <- scaletointerval(distmap(rebound.owin(letterR, W)))\n",
    "Lin <- scaletointerval(distmap(complement.owin(letterR, W)))\n",
    "L <- scaletointerval(eval.im(Lin-Lout))\n",
    "D <- scaletointerval(density(runifpoint(30, W), adjust=0.3))\n",
    "X <- scaletointerval(as.im(function(x,y){ x }, W=W))\n",
    "plot(listof(L=L, D=D, X=X), main=\"Multiple images\")\n",
    "pairs(L, D, X, main=\"Multiple images: pairs plot\")\n",
    "persp(L, colin=D,\n",
    "      theta=-24, phi=35, box=FALSE, apron=TRUE, \n",
    "      main=\"Two images:\\nperspective + colours\",\n",
    "      shade=0.4, ltheta=225, lphi=10)\n",
    "plot(rgbim(D,X,L,maxColorValue=1), valuesAreColours=TRUE,\n",
    "     main=\"Three images: RGB display\")\n",
    "plot(hsvim(D,L,X), valuesAreColours=TRUE,\n",
    "     main=\"Three images: HSV display\")\n",
    "\n",
    "V <- as.linim(function(x,y,seg,tp){(y/1000)^2-(x/1000)^3}, L=domain(chicago))\n",
    "plot(V, main=\"Pixel image on a linear network (colour plot)\")\n",
    "plot(V, style=\"w\", main=\"Pixel image on a linear network (width plot)\")\n",
    "persp(V, phi=25, theta=-35, expand=0.2, main=\"Pixel image on a linear network (perspective plot)\")\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "877c6181",
   "metadata": {},
   "outputs": [],
   "source": [
    "fanfare(\"III. Conversion between types\")\n",
    "\n",
    "W <- as.owin(chorley)\n",
    "plot(W, \"window W\")\n",
    "\n",
    "plot(as.mask(W))\n",
    "plot(as.mask(W, dimyx=1000))\n",
    "\n",
    "plot(as.im(W, value=3))\n",
    "plot(as.im(W, value=3, na.replace=0), ribbon=TRUE)\n",
    "\n",
    "plot(as.im(function(x,y) {x^2 + y}, W=square(1)),\n",
    "     main=\"as.im(function(x,y){x^2+y})\")\n",
    "\n",
    "V <- delaunay(runifpoint(12))\n",
    "plot(V, main=\"Tessellation V\")\n",
    "plot(as.im(V, dimyx=256), main=\"as.im(V)\")\n",
    "plot(as.owin(V))\n",
    "\n",
    "X <- swedishpines\n",
    "plot(X, \"point pattern X\")\n",
    "\n",
    "plot(as.im(X), col=c(\"white\",\"red\"), ribbon=FALSE, xlab=\"\", ylab=\"\")\n",
    "plot(as.owin(X), add=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "036f6c1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "fanfare(\"IV. Subsetting and splitting data\")\n",
    "\n",
    "plot(X, \"point pattern X\")\n",
    "subset <- 1:20\n",
    "plot(X[subset], main=\"subset operation: X[subset]\")\n",
    "subwindow <- owin(poly=list(x=c(0,96,96,40,40),y=c(0,0,100,100,50)))\n",
    "plot(X[subwindow], main=\"subset operation: X[subwindow]\")\n",
    "\n",
    "plot(lansing, \"Lansing Woods data\")\n",
    "plot(split(lansing),\n",
    "     main=\"split operation: split(X)\",\n",
    "     mar.panel=c(0,0,2,0), hsep=1, pch=3)\n",
    "\n",
    "plot(longleaf, main=\"Longleaf Pines data\")\n",
    "plot(cut(longleaf, breaks=3),\n",
    "     main=c(\"cut operation\", \"cut(longleaf, breaks=3)\"))\n",
    "\n",
    "Z <- dirichlet(runifpoint(16))\n",
    "X <- runifpoint(100)\n",
    "\n",
    "plot(cut(X,Z), main=\"points cut by tessellation\", leg.side=\"left\")\n",
    "plot(Z, add=TRUE)\n",
    "\n",
    "plot(split(X, Z),\n",
    "     main=\"points split by tessellation\",\n",
    "     mar.panel=c(0,0,2,2), hsep=1)\n",
    "\n",
    "W <- square(1)\n",
    "X <- as.im(function(x,y){sqrt(x^2+y^2)}, W)\n",
    "Y <- dirichlet(runifpoint(12, W))\n",
    "plot(split(X,Y), main=\"image split by tessellation\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "caa99b95",
   "metadata": {},
   "outputs": [],
   "source": [
    "fanfare(\"V. Tessellations\")\n",
    "\n",
    "plot(quadrats(square(1), nx=5, ny=3))\n",
    "plot(hextess(square(1), 0.07))\n",
    "plot(quantess(letterR, \"x\", 7))\n",
    "plot(polartess(letterR, nangular=6, radii=(0:4)/2, origin=c(2.8, 1.5)),\n",
    "     do.col=TRUE)\n",
    "plot(delaunay(cells))\n",
    "plot(cells, add=TRUE)\n",
    "plot(dirichlet(cells))\n",
    "plot(cells, add=TRUE)\n",
    "plot(rpoislinetess(2.5), do.col=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75f72403",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "fanfare(\"VI. Exploratory data analysis\")\n",
    "\n",
    "par(mar=c(3,3,3,2)+0.1)\n",
    "\n",
    "plot(swedishpines, main=\"Quadrat counts\", pch=\"+\")\n",
    "tab <- quadratcount(swedishpines, 4)\n",
    "plot(tab, add=TRUE, lty=2, cex=2, col=\"blue\")\n",
    "\n",
    "par(mar=c(5,3,3,2)+0.1)\n",
    "\n",
    "plot(swedishpines, main=\"\", pch=\"+\")\n",
    "title(main=expression(chi^2 * \" test\"), cex.main=2)\n",
    "tes <- quadrat.test(swedishpines, 3)\n",
    "tes\n",
    "plot(tes, add=TRUE, col=\"red\", cex=1.5, lty=2, lwd=3)\n",
    "title(sub=paste(\"p-value =\", signif(tes$p.value,3)), cex.sub=1.4)\n",
    "\n",
    "par(mar=c(4,4,3,2)+0.1)\n",
    "\n",
    "tesk <- cdf.test(nztrees, \"x\")\n",
    "tesk\n",
    "plot(tesk)\n",
    "\n",
    "\n",
    "mur <- lapply(murchison, rescale, s=1000)\n",
    "mur <- lapply(mur, \"unitname<-\", value=\"km\")\n",
    "X <- mur$gold\n",
    "D <- distfun(mur$faults)\n",
    "plot(X, main=\"Murchison gold deposits\", cols=\"blue\")\n",
    "plot(mur$faults, add=TRUE, col=\"red\")\n",
    "rh <- rhohat(X,D)\n",
    "plot(rh,\n",
    "     main=\"Smoothed rate estimate\",\n",
    "     xlab=\"Distance to nearest fault (km)\",\n",
    "     legend=FALSE)\n",
    "plot(predict(rh), main=\"predict(rhohat(X,D))\")\n",
    "\n",
    "Z <- density(cells, 0.07)\n",
    "plot(Z, main=\"Kernel smoothed intensity of point pattern\")\n",
    "plot(cells, add=TRUE)\n",
    "\n",
    "plot(redwood, main=\"Redwood data\")\n",
    "te <- scan.test(redwood, 0.1, method=\"poisson\")\n",
    "plot(te, main=c(\"Scan Statistic for redwood data\",\n",
    "              paste(\"p-value =\", signif(te$p.value,3))))\n",
    "plot(redwood, add=TRUE)\n",
    "te\n",
    "\n",
    "X <- unique(unmark(shapley))\n",
    "plot(X, \"Shapley galaxy concentration\", pch=\".\")\n",
    "coco <-colourmap(rev(rainbow(128, end=2/3)), range=c(0,1))\n",
    "pa <- function(i, ...) {\n",
    "  if(i == 1) list(chars=c(\".\", \"+\"), cols=1:2) else\n",
    "             list(size=0.5, pch=16, col=coco)\n",
    "}\n",
    "plot(nnclean(X, k=17), panel.args=pa,\n",
    "     mar.panel=c(0,1,1,0), nrows=2,\n",
    "     main=\"Byers-Raftery nearest neighbour cleaning\",\n",
    "     cex.title=1.2)\n",
    "Y <- sharpen(X, sigma=0.5, edgecorrect=TRUE)\n",
    "plot(Y, main=\"Choi-Hall data sharpening\", pch=\".\")\n",
    "\n",
    "  owpa <- par(mfrow=c(1,2))\n",
    "  W <- grow.rectangle(as.rectangle(letterR), 1)\n",
    "  X <- superimpose(runifpoint(300, letterR),\n",
    "                   runifpoint(50, W), W=W)\n",
    "  plot(W, main=\"clusterset(X, 'm')\")\n",
    "  plot(clusterset(X, 'marks', fast=TRUE), add=TRUE, chars=c(\"o\", \"+\"), cols=1:2)\n",
    "  plot(letterR, add=TRUE)\n",
    "  plot(W, main=\"clusterset(X, 'd')\")\n",
    "  plot(clusterset(X, 'domain', exact=FALSE), add=TRUE)\n",
    "  plot(letterR, add=TRUE)\n",
    "  par(owpa)\n",
    "\n",
    "D <- density(a, sigma=0.05)\n",
    "plot(D, main=\"Kernel smoothed intensity of line segment pattern\")\n",
    "plot(a, add=TRUE)\n",
    "\n",
    "X <- runifpoint(42)\n",
    "plot(dirichlet(X))\n",
    "plot(X, add=TRUE)\n",
    "\n",
    "plot(delaunay(X))\n",
    "plot(X, add=TRUE)\n",
    "\n",
    "parsave <- par(mfrow=c(1,1), mar=0.2+c(0,1,3,1))\n",
    "plot(listof(\"Longleaf Pines data\"=longleaf,\n",
    "            \"Nearest mark\"=nnmark(longleaf),\n",
    "            \"Kernel smoothing of marks\"=Smooth(longleaf,10),\n",
    "            \"Inverse distance weighted\\nsmoothing of marks\"=idw(longleaf)),\n",
    "     equal.scales=TRUE, halign=TRUE, valign=TRUE,\n",
    "     main=\"\", mar.panel=0.2+c(0,0,2,2))\n",
    "par(parsave)\n",
    "\n",
    "fryplot(cells, main=c(\"Fry plot\",\"cells data\"), pch=\"+\")\n",
    "miplot(longleaf, main=\"Morishita Index plot\", pch=16, col=\"blue\")\n",
    "\n",
    "plot(swedishpines, main=\"Swedish Pines data\")\n",
    "K <- Kest(swedishpines)\n",
    "plot(K, main=\"K function for Swedish Pines\", legendmath=TRUE)\n",
    "\n",
    "en <- envelope(swedishpines, fun=Kest, nsim=10, correction=\"translate\")\n",
    "plot(en, main=\"Envelopes of K function based on CSR\", shade=c(\"hi\", \"lo\"))\n",
    "\n",
    "pc <- pcf(swedishpines)\n",
    "plot(pc, main=\"Pair correlation function\")\n",
    "\n",
    "plot(swedishpines, main=\"nearest neighbours\")\n",
    "m <- nnwhich(swedishpines)\n",
    "b <- swedishpines[m]\n",
    "arrows(swedishpines$x, swedishpines$y, b$x, b$y,\n",
    "       angle=12, length=0.1, col=\"red\")\n",
    "\n",
    "plot(swedishpines %mark% nndist(swedishpines),\n",
    "     markscale=1, main=\"Stienen diagram\", legend=FALSE, fg=\"blue\")\n",
    "\n",
    "plot(Gest(swedishpines),\n",
    "     main=c(\"Nearest neighbour distance function G\", \"Gest(swedishpines)\"),\n",
    "     legendmath=TRUE)\n",
    "\n",
    "Z <- distmap(swedishpines, dimyx=512)\n",
    "plot(swedishpines$window, main=\"Distance map\")\n",
    "plot(Z, add=TRUE)\n",
    "points(swedishpines)\n",
    "\n",
    "plot(Fest(swedishpines),\n",
    "     main=c(\"Empty space function F\", \"Fest(swedishpines)\"),\n",
    "     legendmath=TRUE)\n",
    "\n",
    "W <- rebound.owin(letterR, square(5))\n",
    "plot(distmap(W), main=\"Distance map\")\n",
    "plot(W, add=TRUE)\n",
    "\n",
    "a <- psp(runif(20),runif(20),runif(20),runif(20), window=owin())\n",
    "contour(distmap(a), main=\"Distance map\")\n",
    "plot(a, add=TRUE,col=\"red\")\n",
    "\n",
    "plot(Jest(swedishpines), main=c(\"J-function\", \"J(r)=(1-G(r))/(1-F(r))\"))\n",
    "\n",
    "X <- swedishpines\n",
    "X <- X[sample(1:npoints(X))]\n",
    "Z <- nnfun(X)\n",
    "plot(as.owin(X), main=\"Nearest neighbour map\")\n",
    "plot(Z, add=TRUE)\n",
    "points(X)\n",
    "\n",
    "plot(allstats(swedishpines))\n",
    "\n",
    "Fig4b <- residualspaper$Fig4b\n",
    "\n",
    "plot(Fig4b, main=\"Inhomogeneous point pattern\")\n",
    "plot(Kinhom(Fig4b), main=\"Inhomogeneous K-function\")\n",
    "plot(pcfinhom(Fig4b, stoyan=0.1), main=\"Inhomogeneous pair correlation\")\n",
    "plot(Ginhom(Fig4b, sigma=0.06), main=\"Inhomogeneous G-function\")\n",
    "plot(Jinhom(Fig4b, sigma=0.06), main=\"Inhomogeneous J-function\")\n",
    "\n",
    "X <- unmark(bronzefilter)\n",
    "plot(X, \"Bronze filter data\")\n",
    "lam <- predict(ppm(X ~x))\n",
    "plot(Kscaled(X, lam), xlim=c(0, 1.5), main=\"Locally-scaled K function\")\n",
    "\n",
    "plot(urkiola)\n",
    "plot(split(urkiola), cex=0.5)\n",
    "plot(density(split(urkiola)))\n",
    "contour(density(split(urkiola)), panel.begin=as.owin(urkiola))\n",
    "plot(relrisk(urkiola), main=\"Relative risk (cross-validated)\")\n",
    "\n",
    "plot(bramblecanes)\n",
    "br <- rescale(bramblecanes)\n",
    "plot(alltypes(br, \"K\"), mar.panel=c(4,5,2,2)+0.1)\n",
    "\n",
    "ama <- rescale(amacrine)\n",
    "plot(alltypes(ama, Lcross, envelope=TRUE, nsim=9), . - r ~ r, ylim=c(-25, 5))\n",
    "\n",
    "ponderosa.extra$plotit(main=\"Ponderosa Pines\")\n",
    "\n",
    "L <- localL(ponderosa)\n",
    "pL <- plot(L, lty=1, col=1, legend=FALSE,\n",
    "           main=c(\"neighbourhood density functions\",\n",
    "             \"for Ponderosa Pines\"), cex.main=0.8)\n",
    "\n",
    "parsave <- par(mfrow=c(1,2))\n",
    "ponderosa.extra$plotit()\n",
    "par(pty=\"s\")\n",
    "plot(L, iso007 ~ r, main=\"point B\")\n",
    "\n",
    "par(mar=0.2+c(1,1,3,1))\n",
    "ponderosa.extra$plotit()\n",
    "L12 <- localL(ponderosa, rvalue=12)\n",
    "P12 <- ponderosa %mark% L12\n",
    "Z12 <- Smooth(P12, sigma=5, dimyx=128)\n",
    "plot(Z12, col=topo.colors(128),\n",
    "     main=c(\"smoothed\", \"neighbourhood density\"),\n",
    "     cex.main=0.8)\n",
    "contour(Z12, add=TRUE)\n",
    "points(ponderosa, pch=16, cex=0.5)\n",
    "\n",
    "plot(amacrine, main=\"Amacrine cells data\", cex.main=0.8)\n",
    "par(pty=\"s\")\n",
    "mkc <- markcorr(amacrine, \n",
    "                correction=\"translate\", method=\"density\",\n",
    "                kernel=\"epanechnikov\")\n",
    "plot(mkc, main=\"Mark correlation function\", legend=FALSE, cex.main=0.8)\n",
    "par(parsave)\n",
    "\n",
    "par(mar=0.2+c(4,4,3,1))\n",
    "plot(alltypes(amacrine, markconnect), \n",
    "     title=\"Mark connection functions for amacrine cells\")\n",
    "\n",
    "parsave <- par(mfrow=c(1,2))\n",
    "\n",
    "parspruce2 <- par(mar=0.2+c(0,2,2,0))\n",
    "plot(spruces, cex.main=0.8, markscale=10)\n",
    "par(pty=\"s\", mar=0.2+c(2,3,2,0))\n",
    "plot(markcorr(spruces), main=\"Mark correlation\", legendpos=\"bottomright\")\n",
    "\n",
    "par(parspruce2)\n",
    "plot(spruces, cex.main=0.8, markscale=10)\n",
    "par(pty=\"s\", mar=0.2+c(2,3,2,0))\n",
    "plot(markvario(spruces), main=\"Mark variogram\", legendpos=\"topright\")\n",
    "par(parsave)\n",
    "\n",
    "plot(listof(\"Emark(spruces)\"=Emark(spruces),\n",
    "            \"Vmark(spruces)\"=Vmark(spruces)),\n",
    "     main=\"Independence diagnostics\", ylim.covers=0,\n",
    "     legendpos=\"bottom\")\n",
    "\n",
    "par3 <- par(mfrow=c(1,2))\n",
    "X <- rpoispp3(100)\n",
    "plot(X, main=\"3D point pattern X\")\n",
    "plot(K3est(X), main=\"K-function in 3D\")\n",
    "plot(X, main=\"3D point pattern X\")\n",
    "plot(G3est(X), main=\"G-function in 3D\", legendpos=\"bottomright\")\n",
    "par(par3)\n",
    "\n",
    "par(mfrow=c(1,3))\n",
    "X <- unmark(chicago)\n",
    "plot(X, col=\"green\", cols=\"red\", pch=16,\n",
    "     main=\"Chicago Street Crimes\", cex.main=0.75,\n",
    "     show.window=FALSE)\n",
    "plot(density(X, 100, distance=\"e\"), main=\"Kernel density estimate (Euclidean)\")\n",
    "plot(density(X, 100, distance=\"p\"), main=\"Kernel density estimate (shortest path)\")\n",
    "\n",
    "plot(X, col=\"green\", cols=\"red\", pch=16,\n",
    "     main=\"Chicago Street Crimes\", cex.main=0.75,\n",
    "     show.window=FALSE)\n",
    "plot(linearK(X, correction=\"none\"), main=\"Network K-function\", cex.main=0.75)\n",
    "plot(linearK(X, correction=\"Ang\"), main=\"Corrected K-function\", cex.main=0.75)\n",
    "\n",
    "par(mfrow=c(1,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a4b1e48",
   "metadata": {},
   "outputs": [],
   "source": [
    "fanfare(\"VII. Model-fitting\")\n",
    "\n",
    "parsave <- par(mar=0.2+c(1,1,3,2))\n",
    "plot(japanesepines)\n",
    "fit <- ppm(japanesepines ~1)\n",
    "print(fit)\n",
    "fit <- ppm(japanesepines ~polynom(x,y,2))\n",
    "print(fit)\n",
    "plot(fit, how=\"image\", se=FALSE, main=c(\"Inhomogeneous Poisson model\",\n",
    "                               \"fit by maximum likelihood\",\n",
    "                               \"Fitted intensity\"))\n",
    "plot(fit, how=\"image\", trend=FALSE,\n",
    "     main=c(\"Standard error\", \"of fitted intensity\"))\n",
    "\n",
    "plot(leverage(fit))\n",
    "plot(influence(fit))\n",
    "\n",
    "plot(mur$gold, main=\"Murchison gold deposits\", cols=\"blue\")\n",
    "plot(mur$faults, add=TRUE, col=\"red\")\n",
    "fit <- ppm(mur$gold ~D, covariates=list(D=distfun(mur$faults)))\n",
    "par(mar=0.2+c(4,4,4,2))\n",
    "plot(parres(fit, \"D\"),\n",
    "     main=\"Partial residuals from loglinear Poisson model\",\n",
    "     xlab=\"Distance to nearest fault (km)\",\n",
    "     ylab=\"log intensity of gold\", legend=FALSE)\n",
    "legend(\"bottomleft\", legend=c(\"partial residual\", \"loglinear fit\"),\n",
    "       col=c(1,4), lty=c(1,4))\n",
    "\n",
    "par(mar=rep(0.2, 4), mfrow=c(1,1))\n",
    "fitT <- kppm(redwood ~1, clusters=\"Thomas\")\n",
    "simT <- simulate(fitT)[[1]]\n",
    "plot(listof(redwood, simT),\n",
    "     main.panel=c(\"Redwood\", \"simulation from\\nfitted Thomas model\"),\n",
    "     main=\"\", mar.panel=0.2, equal.scales=TRUE)\n",
    "\n",
    "mop <- par(mfrow=c(1,2), pty=\"s\", mar=rep(4.4, 4))\n",
    "plot(fitT, xname=c(\"Thomas model\", \"minimum contrast fit\"), pause=FALSE)\n",
    "par(mop)\n",
    "\n",
    "oop <- par(pty=\"s\", mar=0.2+c(4,4,4,2))\n",
    "os <- objsurf(fitT)\n",
    "plot(os, main=\"Minimum contrast objective function\", col=terrain.colors(128))\n",
    "contour(os, add=TRUE)\n",
    "par(oop)\n",
    "\n",
    "parra <- par(mfrow=c(1,2), mar=0.2+c(3,3,4,2))\n",
    "plot(swedishpines)\n",
    "fit <- ppm(swedishpines ~1, Strauss(r=7))\n",
    "print(fit)\n",
    "plot(fit, how=\"image\", main=c(\"Strauss model\",\n",
    "                               \"fit by maximum pseudolikelihood\",\n",
    "                               \"Conditional intensity plot\"))\n",
    "# fitted interaction\n",
    "plot(swedishpines)\n",
    "fit <- ppm(swedishpines ~1, PairPiece(c(3,5,7,9,11,13)))\n",
    "plot(fitin(fit), legend=FALSE,\n",
    "     main=c(\"Pairwise interaction model\",\n",
    "            \"fit by maximum pseudolikelihood\"))\n",
    "\n",
    "# simulation\n",
    "par(mfrow=c(1,1), mar=0.5+c(0,0,2,0))\n",
    "Xsim <- rmh(model=fit,\n",
    "            start=list(n.start=80),\n",
    "            control=list(nrep=100))\n",
    "plot(listof(swedishpines, Xsim),\n",
    "     main=\"\",\n",
    "     main.panel=c(\"Swedish Pines\",\n",
    "       \"Simulation from\\nfitted Strauss model\"),\n",
    "     mar.panel=c(0,0,3,0),hsep=1,equal.scales=TRUE)\n",
    "\n",
    "# model compensator\n",
    "par(parra)\n",
    "par(mar=0.2+c(4,4,3,1))\n",
    "plot(swedishpines)\n",
    "fit <- ppm(swedishpines ~1, Strauss(r=7))\n",
    "plot(Kcom(fit), cbind(iso, icom, pois) ~ r,\n",
    "     legend=FALSE, main=\"model compensators\")\n",
    "legend(\"topleft\", legend=c(\"empirical K function\",\n",
    "                    \"Strauss model compensator of K\",\n",
    "                    \"Poisson theoretical K\"), lty=1:3, col=1:3, inset=0.05)\n",
    "\n",
    "par(parsave)\n",
    "\n",
    "# Multitype data\n",
    "dpat <- rescale(demopat, 8)\n",
    "unitname(dpat) <- c(\"mile\", \"miles\")\n",
    "dpat\n",
    "\n",
    "plot(dpat, cols=c(\"red\", \"blue\"))\n",
    "fit <- ppm(dpat ~marks + polynom(x,y,2), Poisson())\n",
    "plot(fit, trend=TRUE, se=TRUE)\n",
    "\n",
    "# Linear network data\n",
    "plot(spiders)\n",
    "fit <- lppm(spiders ~ polynom(x,y,2))\n",
    "anova(fit, test=\"Chi\")\n",
    "(fit <- step(fit))\n",
    "lam <- predict(fit)\n",
    "plot(lam, main=\"Point process model on linear network\", ribscale=1000)\n",
    "plot(spiders, add=TRUE, pch=16, show.network=FALSE)\n",
    "\n",
    "fanfare(\"VIII. Simulation\")\n",
    "\n",
    "plot(letterR, main=\"Poisson random points\")\n",
    "lambda <- 10/area.owin(letterR)\n",
    "points(rpoispp(lambda, win=letterR))\n",
    "points(rpoispp(9 * lambda, win=letterR))\n",
    "points(rpoispp(90 * lambda, win=letterR))\n",
    "plot(rpoispp(100))\n",
    "plot(rpoispp(function(x,y){1000 * exp(-3*x)}, 1000),\n",
    "     main=\"rpoispp(function)\")\n",
    "\n",
    "plot(rMaternII(200, 0.05))\n",
    "plot(rSSI(0.05, 200))\n",
    "plot(rThomas(10, 0.2, 5))\n",
    "plot(rMatClust(10, 0.05, 4))\n",
    "plot(rCauchy(30, 0.01, 5))\n",
    "plot(rVarGamma(30, 2, 0.02, 5))\n",
    "plot(rGaussPoisson(30, 0.05, 0.5))\n",
    "\n",
    "if(require(RandomFields) && RandomFieldsSafe()) {\n",
    "  X <- rLGCP(\"exp\", 4, var=0.2, scale=0.1)\n",
    "  plot(attr(X, \"Lambda\"), main=\"log-Gaussian Cox process\")\n",
    "  plot(X, add=TRUE, pch=16)\n",
    "}\n",
    "\n",
    "plot(rStrauss(200, 0.3, 0.07))\n",
    "plot(rDiggleGratton(200,0.03,0.08))\n",
    "plot(rDGS(300, 0.05))\n",
    "\n",
    "plot(redwood, main=\"random thinning - rthin()\")\n",
    "points(rthin(redwood, 0.5), col=\"green\", cex=1.4)\n",
    "\n",
    "plot(rcell(nx=15))\n",
    "\n",
    "plot(rsyst(nx=5))\n",
    "abline(h=(1:4)/5, lty=2)\n",
    "abline(v=(1:4)/5, lty=2)\n",
    "\n",
    "plot(rstrat(nx=5))\n",
    "abline(h=(1:4)/5, lty=2)\n",
    "abline(v=(1:4)/5, lty=2)\n",
    "\n",
    "X <- rsyst(nx=10)\n",
    "plot(rjitter(X, 0.02))\n",
    "\n",
    "Xg <- rmh(list(cif=\"geyer\", par=list(beta=1.25, gamma=1.6, r=0.2, sat=4.5),\n",
    "               w=c(0,10,0,10)),\n",
    "          control=list(nrep=1e4), start=list(n.start=200))\n",
    "plot(Xg, main=paste(\"Geyer saturation process\\n\",\n",
    "                    \"rmh() with cif=\\\"geyer\\\"\"))\n",
    "\n",
    "L <- as.psp(matrix(runif(20), 5, 4), window=square(1))\n",
    "plot(L, main=\"runifpointOnLines(30, L)\")\n",
    "plot(runifpointOnLines(30, L), add=TRUE, pch=\"+\")\n",
    "\n",
    "plot(L, main=\"rpoisppOnLines(3, L)\")\n",
    "plot(rpoisppOnLines(3, L), add=TRUE, pch=\"+\")\n",
    "\n",
    "plot(runiflpp(20, simplenet))\n",
    "plot(rpoislpp(5, simplenet))\n",
    "\n",
    "plot(rpoisline(10))\n",
    "\n",
    "plot(rlinegrid(30, 0.1))\n",
    "\n",
    "spatstat.options(npixel=256)\n",
    "X <- dirichlet(runifpoint(30))\n",
    "plot(rMosaicSet(X, 0.4), col=\"green\", border=NA)\n",
    "plot(X, add=TRUE)\n",
    "plot(rMosaicField(X, runif))\n",
    "plot(rMosaicSet(rpoislinetess(3), 0.5), col=\"green\", border=NA, main=\"Switzer's random set\")\n",
    "spatstat.options(npixel=100)\n",
    "\n",
    "plot(Halton(512, c(2,3)), main=\"quasirandom pattern\")\n",
    "plot(Halton(16384, c(2,3)), main=\"quasirandom pattern\", pch=\".\")\n",
    "\n",
    "fanfare(\"IX. Geometry\")\n",
    "\n",
    "A <- letterR\n",
    "\n",
    "B <- shift(letterR, c(0.2,0.1))\n",
    "plot(boundingbox(A,B), main=\"shift\", type=\"n\")\n",
    "plot(A, add=TRUE)\n",
    "plot(B, add=TRUE, border=\"red\")\n",
    "\n",
    "B <- rotate(letterR, 0.2)\n",
    "plot(boundingbox(A,B), main=\"rotate\", type=\"n\")\n",
    "plot(A, add=TRUE)\n",
    "plot(B, add=TRUE, border=\"red\")\n",
    "\n",
    "mat <- matrix(c(1.1, 0, 0.3, 1), 2, 2)\n",
    "B <- affine(letterR, mat=mat, vec=c(0.2,-0.1))\n",
    "plot(boundingbox(A,B), main=\"affine\", type=\"n\")\n",
    "plot(A, add=TRUE)\n",
    "plot(B, add=TRUE, border=\"red\")\n",
    "\n",
    "par1x2 <- par(mfrow=c(1,2))\n",
    "L <- rpoisline(10, owin(c(1.5,4.5),c(0.2,3.6)))\n",
    "plot(L, main=\"Line segment pattern\")\n",
    "plot(L$window, main=\"L[window]\", type=\"n\")\n",
    "plot(L[letterR], add=TRUE)\n",
    "plot(letterR, add=TRUE, border=\"red\")\n",
    "par(par1x2)\n",
    "\n",
    "a <- psp(runif(20),runif(20),runif(20),runif(20), window=owin())\n",
    "plot(a, main=\"Self-crossing points\")\n",
    "plot(selfcrossing.psp(a), add=TRUE, col=\"red\")\n",
    "\n",
    "a <- as.psp(matrix(runif(20), 5, 4), window=square(1))\n",
    "b <- rstrat(square(1), 5)\n",
    "plot(a, lwd=3, col=\"green\", main=\"project points to segments\")\n",
    "plot(b, add=TRUE, col=\"red\", pch=16)\n",
    "v <- project2segment(b, a)\n",
    "Xproj <- v$Xproj\n",
    "plot(Xproj, add=TRUE, pch=16)\n",
    "arrows(b$x, b$y, Xproj$x, Xproj$y, angle=10, length=0.15, col=\"red\")\n",
    "\n",
    "plot(a, main=\"pointsOnLines(L)\")\n",
    "plot(pointsOnLines(a, np=100), add=TRUE, pch=\"+\")\n",
    "\n",
    "parry <- par(mfrow=c(1,3), mar=0.3+c(1,1,3,1))\n",
    "X <- tess(xgrid=seq(2, 4, length=10), ygrid=seq(0, 3.5, length=8))\n",
    "plot(X, cex.main=0.75)\n",
    "plot(letterR, cex.main=0.75)\n",
    "plot(intersect.tess(X, letterR), cex.main=0.75)\n",
    "\n",
    "X <- dirichlet(runifpoint(10))\n",
    "plot(X)\n",
    "L <- infline(0.3,0.5)\n",
    "plot(owin(), main=\"L\", cex.main=0.75)\n",
    "plot(L, col=\"red\", lwd=2, cex.main=0.75)\n",
    "plot(chop.tess(X,L), cex.main=0.75)\n",
    "par(parry)\n",
    "\n",
    "W <- chorley$window\n",
    "plot(W, main=\"simplify.owin\")\n",
    "WS <- simplify.owin(W, 2)\n",
    "plot(WS, add=TRUE, border=\"green\")\n",
    "\n",
    "nopa <- par(mfrow=c(2,2))\n",
    "Rbox <- grow.rectangle(as.rectangle(letterR), 0.3)\n",
    "\n",
    "v <- erosion.owin(letterR, 0.25)\n",
    "plot(Rbox, type=\"n\", main=\"erode.owin\", cex.main=0.75)\n",
    "plot(letterR, add=TRUE, col=\"red\", cex.main=0.75)\n",
    "plot(v, add=TRUE, col=\"blue\")\n",
    "\n",
    "v <- dilation.owin(letterR, 0.25)\n",
    "plot(Rbox, type=\"n\", main=\"dilate.owin\", cex.main=0.75)\n",
    "plot(v, add=TRUE, col=\"blue\")\n",
    "plot(letterR, add=TRUE, col=\"red\")\n",
    "\n",
    "v <- closing.owin(letterR, 0.3)\n",
    "plot(Rbox, type=\"n\", main=\"closing.owin\", cex.main=0.75)\n",
    "plot(v, add=TRUE, col=\"blue\")\n",
    "plot(letterR, add=TRUE, col=\"red\")\n",
    "\n",
    "v <- opening.owin(letterR, 0.3)\n",
    "plot(Rbox, type=\"n\", main=\"opening.owin\", cex.main=0.75)\n",
    "plot(letterR, add=TRUE, col=\"red\")\n",
    "plot(v, add=TRUE, col=\"blue\")\n",
    "par(nopa)\n",
    "\n",
    "laslett(heather$fine, main=\"Laslett's Transform\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a653364a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fanfare(\"X. Operations on pixel images\")\n",
    "\n",
    "Z <- distmap(swedishpines, dimyx=512)\n",
    "plot(Z, main=\"An image Z\")\n",
    "plot(levelset(Z, 4))\n",
    "plot(cut(Z, 5))\n",
    "plot(eval.im(sqrt(Z) - 3))\n",
    "plot(solutionset(abs(Z - 6) <= 1))\n",
    "nopa <- par(mfrow=c(1,2))\n",
    "plot(Z)\n",
    "segments(0,0,96,100,lwd=2)\n",
    "plot(transect.im(Z))\n",
    "par(nopa)\n",
    "\n",
    "d <- distmap(cells, dimyx=256)\n",
    "W <- levelset(d, 0.06)\n",
    "nopa <- par(mfrow=c(1,2))\n",
    "plot(W)\n",
    "plot(connected(W))\n",
    "par(nopa)\n",
    "\n",
    "Z <- as.im(function(x,y) { 4 * x^2 + 3 * y }, letterR)\n",
    "plot(Z)\n",
    "plot(letterR, add=TRUE)\n",
    "\n",
    "plot(blur(Z, 0.3, bleed=TRUE))\n",
    "plot(letterR, add=TRUE)\n",
    "\n",
    "plot(blur(Z, 0.3, bleed=FALSE))\n",
    "plot(letterR, add=TRUE)\n",
    "          \n",
    "plot(blur(Z, 0.3, bleed=FALSE))\n",
    "plot(letterR, add=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4f99e3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "          \n",
    "fanfare(\"XI. Programming tools\")\n",
    "\n",
    "showoffK <- function(Y, current, ..., fullpicture,rad) { \n",
    "\tplot(fullpicture,\n",
    "             main=c(\"Animation using `applynbd'\", \"explaining the K function\"))\n",
    "\tpoints(Y, cex=2)\n",
    "        u <- current\n",
    "\tpoints(u[1],u[2],pch=\"+\",cex=3)\n",
    "\ttheta <- seq(0,2*pi,length=100)\n",
    "\tpolygon(u[1]+ rad * cos(theta),u[2]+rad*sin(theta))\n",
    "\ttext(u[1]+rad/3,u[2]+rad/2,Y$n,cex=3)\n",
    "        if(runif(1) < 0.2) Sys.sleep(runif(1, max=0.4))\n",
    "\treturn(npoints(Y))\n",
    "}\n",
    "par(ask=FALSE)\n",
    "applynbd(redwood, R=0.2, showoffK, fullpicture=redwood, rad=0.2, exclude=TRUE)\n",
    "\n",
    "par(oldpar)\n",
    "options(oldoptions)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5474f477",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}