# HG changeset patch
# User Martin Janda <janda.martin1@gmail.com>
# Date 1542240009 -3600
#      Thu Nov 15 01:00:09 2018 +0100
# Node ID 76d1155e62fdb2d436a6bb52327a971a7fcab9e6
# Parent  3b1f25568ccdb7cd9797506c04da3efa792f821a
New objects imref2d and imref3d to reference images to world coordinates

* inst/@imref2d/*.m: imref2d constructor and method files
* inst/@imref3d/*.m: imref3d constructor and method files

diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref2d/contains.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref2d/contains.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,63 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} {@var{tf} =} contains (@var{r}, @var{xWorld}, @var{yWorld})
+## Determine if image contains points in world coordinate system.
+##
+## Outputs a logical array @var{tf}, where the i-th nonzero value means the
+## point (@var{xWorld}(i), @var{yWorld}(i)) lies within the bounds of
+## an image associated with a spatial referencing object @var{r}.
+##
+## @seealso{imref2d, imref3d}
+## @end deftypefn
+
+function tf = contains (r, xWorld, yWorld)
+  if (nargin != 3)
+    print_usage();
+  endif
+
+  validateattributes (xWorld, {"numeric"}, ...
+  {"real"}, "imref2d", "xWorld");
+  validateattributes (yWorld, {"numeric"}, ...
+  {"real"}, "imref2d", "yWorld");
+  
+  if (! all (size (xWorld) == size (yWorld)))
+    error ("Octave:invalid-input-arg", ...
+    "xWorld and yWorld must be of the same size");
+  endif
+
+  xWorldLimits = r.XWorldLimits;
+  yWorldLimits = r.YWorldLimits;
+  containsX = xWorld >= xWorldLimits(1) & xWorld <= xWorldLimits(2);
+  containsY = yWorld >= yWorldLimits(1) & yWorld <= yWorldLimits(2);
+  tf = containsX & containsY;
+endfunction
+
+%!error id=Octave:invalid-fun-call contains (imref2d)
+%!error id=Octave:invalid-fun-call contains (imref2d, 1)
+%!error id=Octave:invalid-fun-call contains (imref2d, 1, 2, 3)
+%!error id=Octave:invalid-input-arg contains (imref2d, 1, [2, 3])
+%!error id=Octave:invalid-input-arg contains (imref2d, [1, 2], 3)
+%!error id=Octave:expected-real contains (imref2d, 0, j)
+%!error id=Octave:expected-real contains (imref2d, j, 0)
+
+%!assert (contains (imref2d, [], []), logical( zeros (0, 0)))
+%!assert (contains (imref2d, [1, 2; 3, 4], [5, -6; 7, 8]), logical (zeros (2, 2)))
+
+%!test
+%! r = imref2d ([256, 256]);
+%! assert (contains(r, [5, 8, 8], [5, 10, 257]), logical([1, 1, 0]))
\ No newline at end of file
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref2d/disp.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref2d/disp.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,35 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} disp (@var{r})
+##
+## @seealso{}
+## @end deftypefn
+
+function disp (r)
+  printf("%s with properties:\n", class (r));
+  printf("\n");
+  printf("         XWorldLimits: [%d %d]\n", r.XWorldLimits);
+  printf("         YWorldLimits: [%d %d]\n", r.YWorldLimits);
+  printf("            ImageSize: [%d %d]\n", r.ImageSize);
+  printf("  PixelExtentInWorldX: %d\n", r.PixelExtentInWorldX);
+  printf("  PixelExtentInWorldY: %d\n", r.PixelExtentInWorldY);
+  printf("  ImageExtentInWorldX: %d\n", r.ImageExtentInWorldX);
+  printf("  ImageExtentInWorldY: %d\n", r.ImageExtentInWorldY);
+  printf("     XIntrinsicLimits: [%d %d]\n", r.XIntrinsicLimits);
+  printf("     YIntrinsicLimits: [%d %d]\n", r.YIntrinsicLimits);
+endfunction
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref2d/imref2d.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref2d/imref2d.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,270 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} {@var{r} =} imref2d
+## @deftypefnx {} {@var{r} =} imref2d (@var{imageSize})
+## @deftypefnx {} {@var{r} =} imref2d (@var{imageSize}, @var{pixelExtentInWorldX}, @var{pixelExtentInWorldY})
+## @deftypefnx {} {@var{r} =} imref2d (@var{imageSize}, @var{xWorldLimits}, @var{yWorldLimits})
+## Reference 2-D image to world coordinates.
+##
+## Creates an imref2d object referencing a 2-D m-by-n image with the size
+## @var{imageSize} to world coordinates.  The world extent is either given
+## by @var{xWorldLimits} and @var{yWorldLimits} or computed from
+## @var{pixelExtentInWorldX} and @var{pixelExtentInWorldY}.
+## @var{imageSize} is [2, 2] by default.
+##
+## Intrinsic coordinates are x = 1.0, y = 1.0 in the center of the top left
+## pixel and x = n, y = m in the center of the bottom right pixel.
+## Spatial resolution in each dimension can be different.
+##
+## imref2d object has the following properties:
+##
+## ImageSize - two element integer vector with image height and width
+## in pixels.
+##
+## XWorldLimits - limits of the image along the x-axis in world units
+## specified as a two element real vector.
+##
+## YWorldLimits - limits of the image along the y-axis in world units
+## specified as a two element real vector.
+##
+## PixelExtentInWorldX - pixel extent along the x-axis in world units
+## specified as a real scalar.
+##
+## PixelExtentInWorldY - pixel extent along the y-axis in world units
+## specified as a real scalar.
+##
+## ImageExtentInWorldX - image extent along the x-axis in world units
+## specified as a real scalar.
+##
+## ImageExtentInWorldY - image extent along the y-axis in world units
+## specified as a real scalar.
+##
+## XIntrinsicLimits - limits of the image along the x-axis in intrinsic
+## units, equals to @code{[n - 0.5, n + 0.5]}.
+##
+## YIntrinsicLimits - limits of the image along the y-axis in intrinsic
+## units, equals to @code{[m - 0.5, m + 0.5]}.
+##
+## @seealso{imref3d}
+## @end deftypefn
+
+function r = imref2d (imageSize, varargin)
+  if (nargin > 1 && nargin != 3)
+    print_usage();
+  endif
+
+  if (nargin == 0)
+    imageSize = [2, 2];
+  endif
+
+  if (nargin > 0)
+    if (length (imageSize) < 2)
+      error ("Octave:invalid-input-arg", ...
+      "ImageSize must have at least two elements");
+    endif
+
+    validateattributes (imageSize, {"numeric"}, ...
+    {"positive", "integer", "vector"}, "imref2d", "imageSize");
+  endif
+
+  imageSize = imageSize(1:2);
+  m = imageSize(1);
+  n = imageSize(2);
+
+  if (numel (varargin) == 0)
+    xWorldLimits = [0.5, n + 0.5];
+    yWorldLimits = [0.5, m + 0.5];
+  elseif (numel (varargin) == 2)
+    if (isscalar (varargin{1}))
+      validateattributes (varargin{1}, {"numeric"}, ...
+      {"real", "positive", "scalar"}, "imref2d", "pixelExtentInWorldX");
+      validateattributes (varargin{2}, {"numeric"}, ...
+      {"real", "positive", "scalar"}, "imref2d", "pixelExtentInWorldY");
+      pixelExtentInWorldX = varargin{1};
+      pixelExtentInWorldY = varargin{2};
+    else
+      validateattributes (varargin{1}, {"numeric"}, ...
+      {"real", "vector", "size", [1, 2]}, "imref2d", "xWorldLimits");
+      validateattributes (varargin{2}, {"numeric"}, ...
+      {"real", "vector", "size", [1, 2]}, "imref2d", "yWorldLimits");
+      xWorldLimits = varargin{1};
+      yWorldLimits = varargin{2};
+    endif
+  endif
+
+  if (exist ("pixelExtentInWorldX") && exist ("pixelExtentInWorldY"))
+    imageExtentInWorldX = pixelExtentInWorldX * m;
+    imageExtentInWorldY = pixelExtentInWorldY * n;
+    xWorldLimits = [pixelExtentInWorldX / 2, imageExtentInWorldX + ...
+    pixelExtentInWorldX / 2];
+    yWorldLimits = [pixelExtentInWorldY / 2, imageExtentInWorldY + ...
+    pixelExtentInWorldY / 2];
+  elseif (exist ("xWorldLimits") && exist ("yWorldLimits"))
+    imageExtentInWorldX = xWorldLimits(2) - xWorldLimits(1);
+    imageExtentInWorldY = yWorldLimits(2) - yWorldLimits(1);
+    pixelExtentInWorldX = imageExtentInWorldX / n;
+    pixelExtentInWorldY = imageExtentInWorldY / m;
+  endif
+
+  xIntrinsicLimits = [0.5, n + 0.5];
+  yIntrinsicLimits = [0.5, m + 0.5];
+
+  r.ImageSize = imageSize;
+  r.XWorldLimits = xWorldLimits;
+  r.YWorldLimits = yWorldLimits;
+  r.PixelExtentInWorldX = pixelExtentInWorldX;
+  r.PixelExtentInWorldY = pixelExtentInWorldY;
+  r.ImageExtentInWorldX = imageExtentInWorldX;
+  r.ImageExtentInWorldY = imageExtentInWorldY;
+  r.XIntrinsicLimits = xIntrinsicLimits;
+  r.YIntrinsicLimits = yIntrinsicLimits;
+  r = class (r, "imref2d");
+endfunction
+
+%!error id=Octave:invalid-fun-call imref2d (1, 2, 3, 4)
+%!error id=Octave:invalid-input-arg imref2d (42)
+%!error id=Octave:invalid-input-arg imref2d ([42])
+%!error id=Octave:expected-integer imref2d ([4.2, 42])
+%!error id=Octave:expected-positive imref2d ([0, 0])
+%!error id=Octave:expected-positive imref2d ([-4, 2])
+%!error id=Octave:expected-positive imref2d ([4, 2], 0, 2)
+%!error id=Octave:expected-positive imref2d ([4, 2], 2, 0)
+%!error id=Octave:expected-real imref2d ([4, 2], j, 2)
+%!error id=Octave:expected-real imref2d ([4, 2], 2, j)
+%!error id=Octave:expected-real imref2d ([4, 2], [j, 2], [3, 4])
+%!error id=Octave:expected-real imref2d ([4, 2], [1, 2], [j, 4])
+%!error id=Octave:expected-vector imref2d ([4, 2], [], [])
+%!error id=Octave:expected-vector imref2d ([4, 2], [], [1])
+%!error id=Octave:expected-scalar imref2d ([4, 2], [1], [])
+%!error id=Octave:incorrect-size imref2d ([4, 2], [1, 2], [0])
+%!error id=Octave:incorrect-size imref2d ([4, 2], [1, 2], [1, 2, 3])
+%!error id=Octave:incorrect-size imref2d ([4, 2], [1, 2, 3], [1, 2])
+%!error id=Octave:incorrect-size imref2d ([4, 2], [1; 2], [1, 2])
+%!error id=Octave:incorrect-size imref2d ([4, 2], [1, 2], [1; 2])
+%!error id=Octave:invalid-indexing imref2d().InvalidProperty
+
+%!test
+%! r = imref2d;
+%! assert (r.XWorldLimits, [0.5, 2.5])
+%! assert (r.YWorldLimits, [0.5, 2.5])
+%! assert (r.ImageSize, [2, 2])
+%! assert (r.PixelExtentInWorldX, 1)
+%! assert (r.PixelExtentInWorldY, 1)
+%! assert (r.ImageExtentInWorldX, 2)
+%! assert (r.ImageExtentInWorldY, 2)
+%! assert (r.XIntrinsicLimits, [0.5, 2.5])
+%! assert (r.YIntrinsicLimits, [0.5, 2.5])
+
+%!test
+%! r = imref2d ([100, 200]);
+%! assert (r.XWorldLimits, [0.5, 200.5])
+%! assert (r.YWorldLimits, [0.5, 100.5])
+%! assert (r.ImageSize, [100, 200])
+%! assert (r.PixelExtentInWorldX, 1)
+%! assert (r.PixelExtentInWorldY, 1)
+%! assert (r.ImageExtentInWorldX, 200)
+%! assert (r.ImageExtentInWorldY, 100)
+%! assert (r.XIntrinsicLimits, [0.5, 200.5])
+%! assert (r.YIntrinsicLimits, [0.5, 100.5])
+
+%!test
+%! xWorldLimits = [2, 5];
+%! yWorldLimits = [3, 6];
+%! r = imref2d ([291, 240], xWorldLimits, yWorldLimits);
+%! assert (r.XWorldLimits, [2, 5])
+%! assert (r.YWorldLimits, [3, 6])
+%! assert (r.ImageSize, [291, 240])
+%! assert (r.PixelExtentInWorldX, 0.0125)
+%! assert (r.PixelExtentInWorldY, 0.0103, 1e-3)
+%! assert (r.ImageExtentInWorldX, 3)
+%! assert (r.ImageExtentInWorldY, 3)
+%! assert (r.XIntrinsicLimits, [0.5, 240.5])
+%! assert (r.YIntrinsicLimits, [0.5, 291.5])
+
+%!test
+%! pixelExtentInWorldX = 0.3125;
+%! pixelExtentInWorldY = 0.3125;
+%! r = imref2d ([512, 512], pixelExtentInWorldX, pixelExtentInWorldY);
+%! assert (r.XWorldLimits, [0.15625, 160.1562], 1e-4)
+%! assert (r.YWorldLimits, [0.15625, 160.1562], 1e-4)
+%! assert (r.ImageSize, [512, 512])
+%! assert (r.PixelExtentInWorldX, 0.3125)
+%! assert (r.PixelExtentInWorldY, 0.3125)
+%! assert (r.ImageExtentInWorldX, 160)
+%! assert (r.ImageExtentInWorldY, 160)
+%! assert (r.XIntrinsicLimits, [0.5, 512.5])
+%! assert (r.YIntrinsicLimits, [0.5, 512.5])
+
+## changing ImageSize
+%!test
+%! r = imref2d;
+%! assert (r.XWorldLimits, [0.5, 2.5])
+%! assert (r.YWorldLimits, [0.5, 2.5])
+%! assert (r.ImageSize, [2, 2])
+%! assert (r.PixelExtentInWorldX, 1)
+%! assert (r.PixelExtentInWorldY, 1)
+%! assert (r.ImageExtentInWorldX, 2)
+%! assert (r.ImageExtentInWorldY, 2)
+%! assert (r.XIntrinsicLimits, [0.5, 2.5])
+%! assert (r.YIntrinsicLimits, [0.5, 2.5])
+%! r.ImageSize = [800, 600];
+%! assert (r.XWorldLimits, [0.5, 2.5])
+%! assert (r.YWorldLimits, [0.5, 2.5])
+%! assert (r.ImageSize, [800, 600])
+%! assert (r.PixelExtentInWorldX, 0.003333, 1e-5)
+%! assert (r.PixelExtentInWorldY, 0.0025)
+%! assert (r.ImageExtentInWorldX, 2)
+%! assert (r.ImageExtentInWorldY, 2)
+%! assert (r.XIntrinsicLimits, [0.5, 600.5])
+%! assert (r.YIntrinsicLimits, [0.5, 800.5])
+
+## changing XWorldLimits and YWorldLimits
+%!test
+%! r = imref2d;
+%! assert (r.XWorldLimits, [0.5, 2.5])
+%! assert (r.YWorldLimits, [0.5, 2.5])
+%! assert (r.ImageSize, [2, 2])
+%! assert (r.PixelExtentInWorldX, 1)
+%! assert (r.PixelExtentInWorldY, 1)
+%! assert (r.ImageExtentInWorldX, 2)
+%! assert (r.ImageExtentInWorldY, 2)
+%! assert (r.XIntrinsicLimits, [0.5, 2.5])
+%! assert (r.YIntrinsicLimits, [0.5, 2.5])
+%! r.XWorldLimits = [-60, 13.33];
+%! r.YWorldLimits = [-900.8, -560.26];
+%! assert (r.XWorldLimits, [-60, 13.33])
+%! assert (r.YWorldLimits, [-900.8, -560.26])
+%! assert (r.PixelExtentInWorldX, 36.6650)
+%! assert (r.PixelExtentInWorldY, 170.27, 1e-5)
+%! assert (r.ImageExtentInWorldX, 73.33, 1e-5)
+%! assert (r.ImageExtentInWorldY, 340.54, 1e-5)
+%! assert (r.XIntrinsicLimits, [0.5, 2.5])
+%! assert (r.YIntrinsicLimits, [0.5, 2.5])
+
+%!test
+%! r = imref2d;
+%! fail ("r.XWorldLimits = []", "")
+%! fail ("r.XWorldLimits = [1]", "")
+%! fail ("r.XWorldLimits = [j]", "")
+%! fail ("r.XWorldLimits = [1; 2]", "")
+%! fail ("r.YWorldLimits = []", "")
+%! fail ("r.YWorldLimits = [1]", "")
+%! fail ("r.YWorldLimits = [j]", "")
+%! fail ("r.YWorldLimits = [1; 2]", "")
+
+%!assert (imref2d ([4, 2, 3]).ImageSize, [4, 2]);
\ No newline at end of file
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref2d/intrinsicToWorld.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref2d/intrinsicToWorld.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,78 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} {[@var{xWorld}, @var{yWorld}] =} intrinsicToWorld (@var{r}, @var{xIntrinsic}, @var{yIntrinsic})
+## Convert from intrinsic to world coordinates.
+##
+## Converts intrinsic coordinates of the image associated with the spatial
+## referencing object @var{r} to world coordinates @var{xWorld}
+## and @var{yWorld}.  If a point (@var{xIntrinsic}(i), @var{yIntrinsic}(i))
+## falls outside the intrinsic bounds of the image, the world coordinates are
+## extrapolated, possibly resulting in negative values.
+##
+## @seealso{imref2d, imref3d, worldToIntrinsic}
+## @end deftypefn
+
+function [xWorld, yWorld] = intrinsicToWorld (r, xIntrinsic, yIntrinsic)
+  if (nargin != 3)
+    print_usage ();
+  endif
+  
+  validateattributes (xIntrinsic, {"numeric"}, ...
+  {"real"}, "imref2d", "xIntrinsic");
+  validateattributes (yIntrinsic, {"numeric"}, ...
+  {"real"}, "imref2d", "yIntrinsic");
+  
+  if (! all (size (xIntrinsic) == size (yIntrinsic)))
+    error ("Octave:invalid-input-arg", ...
+    "xIntrinsic and yIntrinsic must be of the same size");
+  endif
+  
+  xIntrinsicLimits = r.XIntrinsicLimits;
+  yIntrinsicLimits = r.YIntrinsicLimits;
+  xWorldLimits = r.XWorldLimits;
+  yWorldLimits = r.YWorldLimits;
+  xWorld = xWorldLimits(1) + r.PixelExtentInWorldX * ...
+  (xIntrinsic - xIntrinsicLimits(1));
+  yWorld = yWorldLimits(1) + r.PixelExtentInWorldY * ...
+  (yIntrinsic - yIntrinsicLimits(1));
+endfunction
+
+%!error id=Octave:invalid-fun-call intrinsicToWorld (imref2d)
+%!error id=Octave:invalid-fun-call intrinsicToWorld (imref2d, 1, 2, 3)
+%!error id=Octave:expected-real intrinsicToWorld (imref2d, 1j, 2)
+%!error id=Octave:expected-real intrinsicToWorld (imref2d, 1, 2j)
+%!error id=Octave:invalid-input-arg intrinsicToWorld (imref2d, [1, 2], 3)
+%!error id=Octave:invalid-input-arg intrinsicToWorld (imref2d, [1], [2, 3])
+
+%!test
+%! r = imref2d ([512, 512], 0.3125, 0.3125);
+%! xIntrinsic = [34, 442];
+%! yIntrinsic = [172, 172];
+%! [xWorld, yWorld] = intrinsicToWorld (r, xIntrinsic, yIntrinsic);
+%! assert (xWorld, [10.625, 138.125])
+%! assert (yWorld, [53.75, 53.75])
+
+%!test
+%! [xWorld, yWorld] = intrinsicToWorld (imref2d, -5.3, -2.8);
+%! assert (xWorld, -5.3)
+%! assert (yWorld, -2.8)
+
+%!test
+%! [xW, yW] = intrinsicToWorld (imref2d, [1, 2; 3, 4],  [2, 3; 5, 9]);
+%! assert (xW, [1, 2; 3, 4])
+%! assert (yW, [2, 3; 5, 9])
\ No newline at end of file
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref2d/sizesMatch.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref2d/sizesMatch.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,65 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} {@var{TF} =} sizesMatch (@var{r}, @var{A})
+## Determine if object and image are size-compatible.
+##
+## Outputs logical 1 (true) if the first two dimensions of an n-dimensional
+## image @var{A} match the image size of the spatial referencing object @var{r},
+## otherwise outputs zero (false).
+##
+## @seealso{imref2d, imref3d}
+## @end deftypefn
+
+function TF = sizesMatch (r, A)
+  if (nargin != 2)
+    print_usage();
+  endif
+
+  sizeA = size(A);
+  TF = all(sizeA(1:2) == r.ImageSize);
+endfunction
+
+%!error id=Octave:invalid-fun-call sizesMatch (imref2d)
+
+## example from MATLAB documentation
+%!test
+%! I = zeros (256, 256);
+%! r = imref2d ([256, 256]);
+%! assert (sizesMatch (r, I), true)
+%! I2 = zeros (246, 300);
+%! assert (sizesMatch (r, I2), false)
+
+## accepts empty image
+%!test
+%! r = imref2d ([256, 256]);
+%! assert (sizesMatch (r, []), false)
+
+## accepts 1-D image
+%!test
+%! r = imref2d ([256, 256]);
+%! assert (sizesMatch (r, 42), false)
+
+## accepts N-D image
+%!test
+%! r = imref2d ([256, 256]);
+%! assert (sizesMatch (r, zeros (256, 256, 3, 2)), true)
+
+%!test
+%! I = zeros (384, 512, 3);
+%! r = imref2d (size (I));
+%! assert (sizesMatch (r, I), true)
\ No newline at end of file
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref2d/subsasgn.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref2d/subsasgn.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,84 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} {@var{rout} =} subsasgn (@var{r}, @var{index}, @var{val})
+##
+## @seealso{}
+## @end deftypefn
+
+function rout = subsasgn (r, index, val)
+  switch (index.type)
+    case "."
+      fld = index.subs;
+      switch (fld)
+        case "ImageSize"
+          imageSize = val;
+
+          if (length (imageSize) < 2)
+            error ("Octave:invalid-input-arg", ...
+            "ImageSize must have at least two elements");
+          endif
+
+          validateattributes (imageSize, {"numeric"}, ...
+          {"positive", "integer", "vector"}, "imref2d", "imageSize");
+          
+          m = imageSize(1);
+          n = imageSize(2);
+          
+          rout = r;
+          rout.ImageSize = imageSize;
+          rout.PixelExtentInWorldX = r.ImageExtentInWorldX / n;
+          rout.PixelExtentInWorldY = r.ImageExtentInWorldY / m;
+          rout.XIntrinsicLimits = [0.5, n + 0.5];
+          rout.YIntrinsicLimits = [0.5, m + 0.5];
+        case "XWorldLimits"
+          xWorldLimits = val;
+          
+          validateattributes (xWorldLimits, {"numeric"}, ...
+          {"increasing", "real", "vector", "size", [1, 2]}, ...
+          "imref2d", "xWorldLimits");
+          
+          imageSize = r.ImageSize;
+          imageExtentInWorldX = xWorldLimits(2) - xWorldLimits(1);
+          
+          rout = r;
+          rout.XWorldLimits = val;
+          rout.ImageExtentInWorldX = imageExtentInWorldX;
+          rout.PixelExtentInWorldX = imageExtentInWorldX / imageSize(2);
+        case "YWorldLimits"
+          yWorldLimits = val;
+          
+          validateattributes (yWorldLimits, {"numeric"}, ...
+          {"increasing", "real", "vector", "size", [1, 2]}, ...
+          "imref2d", "yWorldLimits");
+          
+          imageSize = r.ImageSize;
+          imageExtentInWorldY = yWorldLimits(2) - yWorldLimits(1);
+          
+          rout = r;
+          rout.YWorldLimits = val;
+          rout.ImageExtentInWorldY = imageExtentInWorldY;
+          rout.PixelExtentInWorldY = imageExtentInWorldY / imageSize(1);
+        otherwise
+          error ("Octave:invalid-indexing", ...
+          "@imref2d/subsasgn: invalid property '%s'", fld);
+      endswitch
+
+    otherwise
+      error ("Octave:invalid-indexing", "@imref2d/subsasgn: invalid index type")
+  endswitch
+endfunction
\ No newline at end of file
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref2d/subsref.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref2d/subsref.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,49 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} {@var{r} =} subref (@var{val}, @var{idx})
+##
+## @seealso{}
+## @end deftypefn
+
+function r = subsref (val, idx)
+  if (strcmp (idx.type, "."))
+    switch (idx.subs)
+      case "ImageSize"
+        r = val.ImageSize;
+      case "XWorldLimits"
+        r = val.XWorldLimits;
+      case "YWorldLimits"
+        r = val.YWorldLimits;
+      case "PixelExtentInWorldX"
+        r = val.PixelExtentInWorldX;
+      case "PixelExtentInWorldY"
+        r = val.PixelExtentInWorldY;
+      case "ImageExtentInWorldX"
+        r = val.ImageExtentInWorldX;
+      case "ImageExtentInWorldY"
+        r = val.ImageExtentInWorldY;
+      case "XIntrinsicLimits"
+        r = val.XIntrinsicLimits;
+      case "YIntrinsicLimits"
+        r = val.YIntrinsicLimits;
+      otherwise
+        error ("Octave:invalid-indexing", ...
+        strcat ("unknown property '", idx.subs, "' for class imref2d"));
+    endswitch
+  endif
+endfunction
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref2d/worldToIntrinsic.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref2d/worldToIntrinsic.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,69 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} {[@var{xIntrinsic}, @var{yIntrinsic}] =} worldToIntrinsic (@var{r}, @var{xWorld}, @var{yWorld})
+## Convert from world to intrinsic coordinates.
+##
+## Converts world coordinates @var{xWorld} and @var{yWorld} to intrinsic
+## coordinates @var{xIntrinsic} and @var{yIntrinsic} of an image associated
+## with the spatial referencing object @var{r}.  If a point
+## (@var{xWorld}(i), @var{yWorld}(i)) falls outside the bounds of the image,
+## its intrinsic coordinates are extrapolated, possibly resulting in
+## negative values.
+##
+## @seealso{imref2d, imref3d, intrinsicToWorld}
+## @end deftypefn
+
+function [xIntrinsic, yIntrinsic] = worldToIntrinsic (r, xWorld, yWorld)
+  if (nargin != 3)
+    print_usage ();
+  endif
+  
+  validateattributes (xWorld, {"numeric"}, ...
+  {"real"}, "imref2d", "xWorld");
+  validateattributes (yWorld, {"numeric"}, ...
+  {"real"}, "imref2d", "yWorld");
+  
+  if (! all (size (xWorld) == size (yWorld)))
+    error ("Octave:invalid-input-arg", ...
+    "xWorld and yWorld must be of the same size");
+  endif
+  
+  xIntrinsicLimits = r.XIntrinsicLimits;
+  yIntrinsicLimits = r.YIntrinsicLimits;
+  xWorldLimits = r.XWorldLimits;
+  yWorldLimits = r.YWorldLimits;
+  xIntrinsic = xIntrinsicLimits(1) + (xWorld  - xWorldLimits(1)) ...
+  / r.PixelExtentInWorldX;
+  yIntrinsic = yIntrinsicLimits(1) + (yWorld  - yWorldLimits(1)) ...
+  / r.PixelExtentInWorldY;
+endfunction
+
+%!error id=Octave:invalid-fun-call worldToIntrinsic (imref2d)
+%!error id=Octave:invalid-fun-call worldToIntrinsic (imref2d, 1, 2, 3)
+%!error id=Octave:expected-real worldToIntrinsic (imref2d, 1j, 2)
+%!error id=Octave:expected-real worldToIntrinsic (imref2d, 1, 2j)
+%!error id=Octave:invalid-input-arg worldToIntrinsic (imref2d, [1, 2], 3)
+%!error id=Octave:invalid-input-arg worldToIntrinsic (imref2d, [1], [2, 3])
+
+%!test
+%! r = imref2d ([512, 512], 0.3125, 0.3125);
+%! xW = [38.44, 39.44, 38.44, -0.2];
+%! yW = [68.75, 68.75, 75.75, -1];
+%! [xI, yI] = worldToIntrinsic (r, xW, yW);
+%! assert (xI, [123.008, 126.208, 123.008, -0.64], 1e-6)
+%! assert (yI, [220, 220, 242.4, -3.2], 1e-6)
\ No newline at end of file
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref2d/worldToSubscript.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref2d/worldToSubscript.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,69 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} {[@var{i}, @var{j}] =} worldToSubscript (@var{r}, @var{xWorld}, @var{yWorld})
+## Convert world coordinates to row and column subscripts.
+##
+## Converts world coordinates to row and column subscripts of an image
+## associated with the spatial referencing object @var{r}.  A point located at
+## (@var{xWorld}(i), @var{yWorld}(i)) world coordinates maps to row and column
+## subscripts @var{i}(i) and @var{j}(i) respectively.  Note the reversed order
+## of the dimensions.  If the point falls outside the bounds of the image, both
+## its row and column subscripts are NaN.
+##
+## @seealso{imref2d, imref3d, worldToIntrinsic}
+## @end deftypefn
+
+function [i, j] = worldToSubscript (r, xWorld, yWorld)
+  if (nargin != 3)
+    print_usage ();
+  endif
+
+  validateattributes (xWorld, {"numeric"}, ...
+  {"real"}, "imref2d", "xWorld");
+  validateattributes (yWorld, {"numeric"}, ...
+  {"real"}, "imref2d", "yWorld");
+  
+  if (! all (size (xWorld) == size (yWorld)))
+    error ("Octave:invalid-input-arg", ...
+    "xWorld and yWorld must be of the same size");
+  endif
+
+  [xIntrinsic, yIntrinsic] = worldToIntrinsic (r, xWorld, yWorld);
+  xIntrinsicLimits = r.XIntrinsicLimits;
+  yIntrinsicLimits = r.YIntrinsicLimits;
+  inImage = contains (r, xWorld, yWorld);
+  xIntrinsic(! inImage) = NaN;
+  yIntrinsic(! inImage) = NaN;
+  i = round (yIntrinsic);
+  j = round (xIntrinsic);
+endfunction
+
+%!error id=Octave:invalid-fun-call worldToSubscript (imref2d)
+%!error id=Octave:invalid-fun-call worldToSubscript (imref2d, 1, 2, 3)
+%!error id=Octave:expected-real worldToSubscript (imref2d, 1j, 2)
+%!error id=Octave:expected-real worldToSubscript (imref2d, 1, 2j)
+%!error id=Octave:invalid-input-arg worldToSubscript (imref2d, [1, 2], 3)
+%!error id=Octave:invalid-input-arg worldToSubscript (imref2d, [1], [2, 3])
+
+%!test
+%! r = imref2d ([512, 512], 0.3125, 0.3125);
+%! xW = [38.44, 39.44, 38.44, -0.2];
+%! yW = [68.75, 68.75, 75.75, -1];
+%! [rS, cS] = worldToSubscript (r, xW, yW);
+%! assert (rS, [220, 220, 242, NaN])
+%! assert (cS, [123, 126, 123, NaN])
\ No newline at end of file
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref3d/contains.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref3d/contains.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,68 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} {@var{tf} =} contains (@var{r}, @var{xWorld}, @var{yWorld}, @var{zWorld})
+## Determine if image contains points in world coordinate system.
+##
+## Outputs a logical array @var{tf}, where the i-th nonzero value means the
+## point (@var{xWorld}(i), @var{yWorld}(i), @var{zWorld}(i)) lies within the
+## bounds of an image associated with a spatial referencing object @var{r}.
+##
+## @seealso{imref2d, imref3d}
+## @end deftypefn
+
+function tf = contains (r, xWorld, yWorld, zWorld)
+  if (nargin != 4)
+    print_usage();
+  endif
+
+  validateattributes (xWorld, {"numeric"}, ...
+  {"real"}, "imref3d", "xWorld");
+  validateattributes (yWorld, {"numeric"}, ...
+  {"real"}, "imref3d", "yWorld");
+  validateattributes (zWorld, {"numeric"}, ...
+  {"real"}, "imref3d", "zWorld");
+
+  if (! all (size (xWorld) == size (yWorld)) ...
+    || ! all (size (xWorld) == size (zWorld)))
+    error ("Octave:invalid-input-arg", ...
+    "imref3d/contains: xWorld, yWorld and zWorld must be of the same size");
+  endif
+  
+  xWorldLimits = r.XWorldLimits;
+  yWorldLimits = r.YWorldLimits;
+  zWorldLimits = r.ZWorldLimits;
+  containsX = xWorld >= xWorldLimits(1) & xWorld <= xWorldLimits(2);
+  containsY = yWorld >= yWorldLimits(1) & yWorld <= yWorldLimits(2);
+  containsZ = zWorld >= zWorldLimits(1) & zWorld <= zWorldLimits(2);
+  tf = containsX & containsY & containsZ;
+endfunction
+
+%!error id=Octave:invalid-fun-call contains (imref3d)
+%!error id=Octave:invalid-fun-call contains (imref3d, 1)
+%!error id=Octave:invalid-fun-call contains (imref3d, 1, 2)
+%!error id=Octave:invalid-fun-call contains (imref3d, 1, 2, 3, 4)
+%!error id=Octave:invalid-input-arg contains (imref3d, [1, 2], 3, 4)
+%!error id=Octave:invalid-input-arg contains (imref3d, 1, [2, 3], 4)
+%!error id=Octave:invalid-input-arg contains (imref3d, 1, 2, [3, 4])
+%!error id=Octave:expected-real contains (imref3d, 1j, 2, 3)
+%!error id=Octave:expected-real contains (imref3d, 1, 2j, 3)
+%!error id=Octave:expected-real contains (imref3d, 1, 2, 3j)
+
+%!test
+%! r = imref3d ([128, 128, 27]);
+%! assert (contains (r, [5, 6, 6, 8], [5, 10, 10, 257], [1, 27.5, 28, 1]), logical ([1, 1, 0, 0]))
\ No newline at end of file
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref3d/disp.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref3d/disp.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,39 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} disp (@var{r})
+##
+## @seealso{}
+## @end deftypefn
+
+function disp (r)
+  printf("%s with properties:\n", class (r));
+  printf("\n");
+  printf("         XWorldLimits: [%d %d]\n", r.XWorldLimits);
+  printf("         YWorldLimits: [%d %d]\n", r.YWorldLimits);
+  printf("         ZWorldLimits: [%d %d]\n", r.ZWorldLimits);
+  printf("            ImageSize: [%d %d %d]\n", r.ImageSize);
+  printf("  PixelExtentInWorldX: %d\n", r.PixelExtentInWorldX);
+  printf("  PixelExtentInWorldY: %d\n", r.PixelExtentInWorldY);
+  printf("  PixelExtentInWorldZ: %d\n", r.PixelExtentInWorldZ);
+  printf("  ImageExtentInWorldX: %d\n", r.ImageExtentInWorldX);
+  printf("  ImageExtentInWorldY: %d\n", r.ImageExtentInWorldY);
+  printf("  ImageExtentInWorldZ: %d\n", r.ImageExtentInWorldZ);
+  printf("     XIntrinsicLimits: [%d %d]\n", r.XIntrinsicLimits);
+  printf("     YIntrinsicLimits: [%d %d]\n", r.YIntrinsicLimits);
+  printf("     ZIntrinsicLimits: [%d %d]\n", r.ZIntrinsicLimits);
+endfunction
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref3d/imref3d.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref3d/imref3d.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,324 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} {@var{r} =} imref3d
+## @deftypefnx {} {@var{r} =} imref3d (@var{imageSize})
+## @deftypefnx {} {@var{r} =} imref3d (@var{imageSize}, @var{pixelExtentInWorldX}, @var{pixelExtentInWorldY}, @var{pixelExtentInWorldZ})
+## @deftypefnx {} {@var{r} =} imref3d (@var{imageSize}, @var{xWorldLimits}, @var{yWorldLimits}, @var{zWorldLimits})
+## Reference 3-D image to world coordinates.
+##
+## Creates an imref3d object referencing a 3-D m-by-n-by-p image with the size
+## @var{imageSize} to world coordinates.  The world extent is either given
+## by @var{xWorldLimits}, @var{yWorldLimits} and @var{xWorldLimits} or computed
+## from @var{pixelExtentInWorldX}, @var{pixelExtentInWorldY} and
+## @var{pixelExtentInWorldZ}.  @var{imageSize} is [2, 2, 2] by default.
+##
+## Intrinsic coordinates are x = 1.0, y = 1.0, z = 1.0 in the center of the
+## top left pixel in the first plane and x = n, y = m, z = p in the center
+## of the bottom right pixel in the last plane.  Spatial resolution in each
+## dimension can be different.
+##
+## imref3d object has the following properties:
+##
+## ImageSize - two element integer vector with image height and width
+## in pixels.
+##
+## XWorldLimits - limits of the image along the x-axis in world units
+## specified as a two element real vector.
+##
+## YWorldLimits - limits of the image along the y-axis in world units
+## specified as a two element real vector.
+##
+## ZWorldLimits - limits of the image along the z-axis in world units
+## specified as a two element real vector.
+##
+## PixelExtentInWorldX - pixel extent along the x-axis in world units
+## specified as a real scalar.
+##
+## PixelExtentInWorldY - pixel extent along the y-axis in world units
+## specified as a real scalar.
+##
+## PixelExtentInWorldZ - pixel extent along the z-axis in world units
+## specified as a real scalar.
+##
+## ImageExtentInWorldX - image extent along the x-axis in world units
+## specified as a real scalar.
+##
+## ImageExtentInWorldY - image extent along the y-axis in world units
+## specified as a real scalar.
+##
+## ImageExtentInWorldZ - image extent along the z-axis in world units
+## specified as a real scalar.
+##
+## XIntrinsicLimits - limits of the image along the x-axis in intrinsic
+## units, equals to @code{[n - 0.5, n + 0.5]}.
+##
+## YIntrinsicLimits - limits of the image along the y-axis in intrinsic
+## units, equals to @code{[m - 0.5, m + 0.5]}.
+##
+## ZIntrinsicLimits - limits of the image along the z-axis in intrinsic
+## units, equals to @code{[p - 0.5, p + 0.5]}.
+##
+## @seealso{imref2d}
+## @end deftypefn
+
+function r = imref3d (imageSize, varargin)
+  if (nargin > 4)
+    print_usage ();
+  endif
+
+  if (nargin == 0)
+    imageSize = [2, 2, 2];
+  elseif (nargin > 0)
+    validateattributes (imageSize, {"numeric"}, ...
+    {"positive", "integer", "vector", "size", [1, 3]}, "imref3d", "imageSize");
+    
+    imageSize = imageSize(1:3);
+  endif
+  
+  m = imageSize(1);
+  n = imageSize(2);
+  p = imageSize(3);
+  
+  if (numel (varargin) == 0)
+    xWorldLimits = [0.5, n + 0.5];
+    yWorldLimits = [0.5, m + 0.5];
+    zWorldLimits = [0.5, p + 0.5];
+    r2 = @imref2d (imageSize);
+  elseif (numel (varargin) == 3)
+    if (isscalar (varargin{1}))
+      validateattributes (varargin{1}, {"numeric"}, ...
+      {"real", "positive", "scalar"}, "imref3d", "pixelExtentInWorldX");
+      validateattributes (varargin{2}, {"numeric"}, ...
+      {"real", "positive", "scalar"}, "imref3d", "pixelExtentInWorldY");
+      validateattributes (varargin{3}, {"numeric"}, ...
+      {"real", "positive", "scalar"}, "imref3d", "pixelExtentInWorldZ");
+      pixelExtentInWorldX = varargin{1};
+      pixelExtentInWorldY = varargin{2};
+      pixelExtentInWorldZ = varargin{3};
+    else
+      validateattributes (varargin{1}, {"numeric"}, ...
+      {"real", "vector", "size", [1, 2]}, "imref3d", "xWorldLimits");
+      validateattributes (varargin{2}, {"numeric"}, ...
+      {"real", "vector", "size", [1, 2]}, "imref3d", "yWorldLimits");
+      validateattributes (varargin{3}, {"numeric"}, ...
+      {"real", "vector", "size", [1, 2]}, "imref3d", "zWorldLimits");
+      xWorldLimits = varargin{1};
+      yWorldLimits = varargin{2};
+      zWorldLimits = varargin{3};
+    endif
+  endif
+
+  if (exist ("pixelExtentInWorldX") && exist ("pixelExtentInWorldY") ...
+    && exist ("pixelExtentInWorldZ"))
+    imageExtentInWorldX = pixelExtentInWorldX * m;
+    imageExtentInWorldY = pixelExtentInWorldY * n;
+    imageExtentInWorldZ = pixelExtentInWorldZ * p;
+    xWorldLimits = [pixelExtentInWorldX / 2, imageExtentInWorldX + ...
+    pixelExtentInWorldX / 2];
+    yWorldLimits = [pixelExtentInWorldY / 2, imageExtentInWorldY + ...
+    pixelExtentInWorldY / 2];
+    zWorldLimits = [pixelExtentInWorldZ / 2, imageExtentInWorldZ + ...
+    pixelExtentInWorldZ / 2];
+  elseif (exist ("xWorldLimits") && exist ("yWorldLimits") ...
+    && exist ("zWorldLimits"))
+    imageExtentInWorldX = xWorldLimits(2) - xWorldLimits(1);
+    imageExtentInWorldY = yWorldLimits(2) - yWorldLimits(1);
+    imageExtentInWorldZ = zWorldLimits(2) - zWorldLimits(1);
+    pixelExtentInWorldX = imageExtentInWorldX / n;
+    pixelExtentInWorldY = imageExtentInWorldY / m;
+    pixelExtentInWorldZ = imageExtentInWorldZ / p;
+  endif
+
+  xIntrinsicLimits = [0.5, n + 0.5];
+  yIntrinsicLimits = [0.5, m + 0.5];
+  zIntrinsicLimits = [0.5, p + 0.5];
+
+  r.ImageSize = imageSize;
+  r.XWorldLimits = xWorldLimits;
+  r.YWorldLimits = yWorldLimits;
+  r.ZWorldLimits = zWorldLimits;
+  r.PixelExtentInWorldX = pixelExtentInWorldX;
+  r.PixelExtentInWorldY = pixelExtentInWorldY;
+  r.PixelExtentInWorldZ = pixelExtentInWorldZ;
+  r.ImageExtentInWorldX = imageExtentInWorldX;
+  r.ImageExtentInWorldY = imageExtentInWorldY;
+  r.ImageExtentInWorldZ = imageExtentInWorldZ;
+  r.XIntrinsicLimits = xIntrinsicLimits;
+  r.YIntrinsicLimits = yIntrinsicLimits;
+  r.ZIntrinsicLimits = zIntrinsicLimits;
+  
+  ## in MATLAB imref3d isa imref2d
+  r = class (r, "imref3d", @imref2d());
+endfunction
+
+%!error id=Octave:invalid-fun-call imref3d (1, 2, 3, 4, 5)
+%!error id=Octave:incorrect-size imref3d (42)
+%!error id=Octave:incorrect-size imref3d ([42])
+%!error id=Octave:incorrect-size imref3d ([4, 2])
+%!error id=Octave:incorrect-size imref3d ([4, 2, 3, 3])
+%!error id=Octave:expected-integer imref3d ([4.2, 42])
+%!error id=Octave:expected-positive imref3d ([0, 0])
+%!error id=Octave:expected-positive imref3d ([-4, 2])
+%!error id=Octave:expected-positive imref3d ([4, 2, 3], 0, 1, 2)
+%!error id=Octave:expected-positive imref3d ([4, 2, 3], 1, 0, 2)
+%!error id=Octave:expected-positive imref3d ([4, 2, 3], 1, 2, 0)
+%!error id=Octave:expected-real imref3d ([4, 2, 3], j, 1, 2)
+%!error id=Octave:expected-real imref3d ([4, 2, 3], 1, j, 2)
+%!error id=Octave:expected-real imref3d ([4, 2, 3], 1, 2, j)
+%!error id=Octave:expected-real imref3d ([4, 2, 3], [j, 2], [3, 4], [5, 6])
+%!error id=Octave:expected-real imref3d ([4, 2, 3], [1, 2], [j, 4], [5, 6])
+%!error id=Octave:expected-real imref3d ([4, 2, 3], [1, 2], [3, 4], [5, j])
+%!error id=Octave:expected-vector imref3d ([4, 2, 3], [], [], [])
+%!error id=Octave:expected-vector imref3d ([4, 2, 3], [], [1], [2])
+%!error id=Octave:expected-scalar imref3d ([4, 2, 3], [1], [], [])
+%!error id=Octave:incorrect-size imref3d ([4, 2, 3], [1, 2], [3, 4], [0])
+%!error id=Octave:incorrect-size imref3d ([4, 2, 3], [1, 2], [3, 4, 5], [6, 7])
+%!error id=Octave:incorrect-size imref3d ([4, 2, 3], [1, 2], [3, 4], [5, 6, 7])
+%!error id=Octave:incorrect-size imref3d ([4, 2, 3], [1; 2], [3, 4], [5, 6])
+%!error id=Octave:incorrect-size imref3d ([4, 2, 3], [1, 2], [3; 4], [5, 6])
+%!error id=Octave:incorrect-size imref3d ([4, 2, 3], [1, 2], [3, 4], [5; 6])
+%!error id=Octave:invalid-indexing imref3d().InvalidProperty
+
+%!assert (imref3d ([4, 2, 3]).ImageSize, [4, 2, 3])
+
+%!test
+%! r = imref3d;
+%! assert (r.XWorldLimits, [0.5, 2.5])
+%! assert (r.YWorldLimits, [0.5, 2.5])
+%! assert (r.ZWorldLimits, [0.5, 2.5])
+%! assert (r.ImageSize, [2, 2, 2])
+%! assert (r.PixelExtentInWorldX, 1)
+%! assert (r.PixelExtentInWorldY, 1)
+%! assert (r.PixelExtentInWorldZ, 1)
+%! assert (r.ImageExtentInWorldX, 2)
+%! assert (r.ImageExtentInWorldY, 2)
+%! assert (r.ImageExtentInWorldZ, 2)
+%! assert (r.XIntrinsicLimits, [0.5, 2.5])
+%! assert (r.YIntrinsicLimits, [0.5, 2.5])
+%! assert (r.ZIntrinsicLimits, [0.5, 2.5])
+
+%!test
+%! r = imref3d ([128, 128, 27]);
+%! assert (r.XWorldLimits, [0.5, 128.5])
+%! assert (r.YWorldLimits, [0.5, 128.5])
+%! assert (r.ZWorldLimits, [0.5, 27.5])
+%! assert (r.ImageSize, [128, 128, 27])
+%! assert (r.PixelExtentInWorldX, 1)
+%! assert (r.PixelExtentInWorldY, 1)
+%! assert (r.PixelExtentInWorldZ, 1)
+%! assert (r.ImageExtentInWorldX, 128)
+%! assert (r.ImageExtentInWorldY, 128)
+%! assert (r.ImageExtentInWorldZ, 27)
+%! assert (r.XIntrinsicLimits, [0.5, 128.5])
+%! assert (r.YIntrinsicLimits, [0.5, 128.5])
+%! assert (r.ZIntrinsicLimits, [0.5, 27.5])
+
+%!test
+%! r = imref3d ([128, 128, 27], 2, 2, 4);
+%! assert (r.XWorldLimits, [1, 257])
+%! assert (r.YWorldLimits, [1, 257])
+%! assert (r.ZWorldLimits, [2, 110])
+%! assert (r.ImageSize, [128, 128, 27])
+%! assert (r.PixelExtentInWorldX, 2)
+%! assert (r.PixelExtentInWorldY, 2)
+%! assert (r.PixelExtentInWorldZ, 4)
+%! assert (r.ImageExtentInWorldX, 256)
+%! assert (r.ImageExtentInWorldY, 256)
+%! assert (r.ImageExtentInWorldZ, 108)
+%! assert (r.XIntrinsicLimits, [0.5, 128.5])
+%! assert (r.YIntrinsicLimits, [0.5, 128.5])
+%! assert (r.ZIntrinsicLimits, [0.5, 27.5])
+
+## changing ImageSize
+%!test
+%! r = imref3d;
+%! assert (r.XWorldLimits, [0.5, 2.5])
+%! assert (r.YWorldLimits, [0.5, 2.5])
+%! assert (r.ZWorldLimits, [0.5, 2.5])
+%! assert (r.ImageSize, [2, 2, 2])
+%! assert (r.PixelExtentInWorldX, 1)
+%! assert (r.PixelExtentInWorldY, 1)
+%! assert (r.PixelExtentInWorldZ, 1)
+%! assert (r.ImageExtentInWorldX, 2)
+%! assert (r.ImageExtentInWorldY, 2)
+%! assert (r.ImageExtentInWorldZ, 2)
+%! assert (r.XIntrinsicLimits, [0.5, 2.5])
+%! assert (r.YIntrinsicLimits, [0.5, 2.5])
+%! assert (r.ZIntrinsicLimits, [0.5, 2.5])
+%! r.ImageSize = [128, 128, 27];
+%! assert (r.XWorldLimits, [0.5, 2.5])
+%! assert (r.YWorldLimits, [0.5, 2.5])
+%! assert (r.ZWorldLimits, [0.5, 2.5])
+%! assert (r.ImageSize, [128, 128, 27])
+%! assert (r.PixelExtentInWorldX, 0.015625, 1e-6)
+%! assert (r.PixelExtentInWorldY, 0.015625, 1e-6)
+%! assert (r.PixelExtentInWorldZ, 0.074074, 1e-6)
+%! assert (r.ImageExtentInWorldX, 2)
+%! assert (r.ImageExtentInWorldY, 2)
+%! assert (r.ImageExtentInWorldZ, 2)
+%! assert (r.XIntrinsicLimits, [0.5, 128.5])
+%! assert (r.YIntrinsicLimits, [0.5, 128.5])
+%! assert (r.ZIntrinsicLimits, [0.5, 27.5])
+
+## changing XWorldLimits, YWorldLimits and ZWorldLimits
+%!test
+%! r = imref3d;
+%! assert (r.XWorldLimits, [0.5, 2.5])
+%! assert (r.YWorldLimits, [0.5, 2.5])
+%! assert (r.ZWorldLimits, [0.5, 2.5])
+%! assert (r.ImageSize, [2, 2, 2])
+%! assert (r.PixelExtentInWorldX, 1)
+%! assert (r.PixelExtentInWorldY, 1)
+%! assert (r.PixelExtentInWorldZ, 1)
+%! assert (r.ImageExtentInWorldX, 2)
+%! assert (r.ImageExtentInWorldY, 2)
+%! assert (r.ImageExtentInWorldZ, 2)
+%! assert (r.XIntrinsicLimits, [0.5, 2.5])
+%! assert (r.YIntrinsicLimits, [0.5, 2.5])
+%! assert (r.ZIntrinsicLimits, [0.5, 2.5])
+%! r.XWorldLimits = [-60, 13.33];
+%! r.YWorldLimits = [-900.8, -560.26];
+%! r.ZWorldLimits = [-302.48, 1500.333];
+%! assert (r.XWorldLimits, [-60, 13.33])
+%! assert (r.YWorldLimits, [-900.8, -560.26])
+%! assert (r.ZWorldLimits, [-302.48, 1500.333])
+%! assert (r.ImageSize, [2, 2, 2])
+%! assert (r.PixelExtentInWorldX, 36.6650)
+%! assert (r.PixelExtentInWorldY, 170.27, 1e-5)
+%! assert (r.PixelExtentInWorldZ, 901.4065)
+%! assert (r.ImageExtentInWorldX, 73.33, 1e-5)
+%! assert (r.ImageExtentInWorldY, 340.54, 1e-5)
+%! assert (r.ImageExtentInWorldZ, 1802.813, 1e-5)
+%! assert (r.XIntrinsicLimits, [0.5, 2.5])
+%! assert (r.YIntrinsicLimits, [0.5, 2.5])
+%! assert (r.ZIntrinsicLimits, [0.5, 2.5])
+
+%!test
+%! r = imref3d;
+%! fail ("r.XWorldLimits = []", "")
+%! fail ("r.XWorldLimits = [1]", "")
+%! fail ("r.XWorldLimits = [j]", "")
+%! fail ("r.XWorldLimits = [1; 2]", "")
+%! fail ("r.YWorldLimits = []", "")
+%! fail ("r.YWorldLimits = [1]", "")
+%! fail ("r.YWorldLimits = [j]", "")
+%! fail ("r.YWorldLimits = [1; 2]", "")
+%! fail ("r.ZWorldLimits = []", "")
+%! fail ("r.ZWorldLimits = [1]", "")
+%! fail ("r.ZWorldLimits = [j]", "")
+%! fail ("r.ZWorldLimits = [1; 2]", "")
\ No newline at end of file
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref3d/intrinsicToWorld.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref3d/intrinsicToWorld.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,96 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} {[@var{xWorld}, @var{yWorld}, @var{zWorld}] =} intrinsicToWorld (@var{r}, @var{xIntrinsic}, @var{yIntrinsic}, @var{zIntrinsic})
+## Convert from intrinsic to world coordinates.
+##
+## Converts intrinsic coordinates of the image associated with the spatial
+## referencing object @var{r} to world coordinates @var{xWorld}, @var{yWorld}
+## and @var{zWorld}.  If a point
+## (@var{xIntrinsic}(i), @var{yIntrinsic}(i), @var{zIntrinsic}(i))
+## falls outside the intrinsic bounds of the image, the world coordinates are
+## extrapolated, possibly resulting in negative values.
+##
+## @seealso{imref2d, imref3d, worldToIntrinsic}
+## @end deftypefn
+
+function [xWorld, yWorld, zWorld] = intrinsicToWorld (r, xIntrinsic, yIntrinsic, zIntrinsic)
+  if (nargin != 4)
+    print_usage ();
+  endif
+  
+  validateattributes (xIntrinsic, {"numeric"}, ...
+  {"real"}, "imref3d", "xIntrinsic");
+  validateattributes (yIntrinsic, {"numeric"}, ...
+  {"real"}, "imref3d", "yIntrinsic");
+  validateattributes (zIntrinsic, {"numeric"}, ...
+  {"real"}, "imref3d", "zIntrinsic");
+  
+  if (! all (size (xIntrinsic) == size (yIntrinsic)) ...
+    || ! all (size (xIntrinsic) == size (zIntrinsic)))
+    error ("Octave:invalid-input-arg", ...
+    "xIntrinsic, yIntrinsic and zIntrinsic must be of the same size");
+  endif
+  
+  xIntrinsicLimits = r.XIntrinsicLimits;
+  yIntrinsicLimits = r.YIntrinsicLimits;
+  zIntrinsicLimits = r.ZIntrinsicLimits;
+  xWorldLimits = r.XWorldLimits;
+  yWorldLimits = r.YWorldLimits;
+  zWorldLimits = r.ZWorldLimits;
+  xWorld = xWorldLimits(1) + r.PixelExtentInWorldX * ...
+  (xIntrinsic - xIntrinsicLimits(1));
+  yWorld = yWorldLimits(1) + r.PixelExtentInWorldY * ...
+  (yIntrinsic - yIntrinsicLimits(1));
+  zWorld = zWorldLimits(1) + r.PixelExtentInWorldZ * ...
+  (zIntrinsic - zIntrinsicLimits(1));
+endfunction
+
+%!error id=Octave:invalid-fun-call intrinsicToWorld (imref3d)
+%!error id=Octave:invalid-fun-call intrinsicToWorld (imref3d, 1)
+%!error id=Octave:invalid-fun-call intrinsicToWorld (imref3d, 1, 2)
+%!error id=Octave:invalid-fun-call intrinsicToWorld (imref3d, 1, 2, 3, 4)
+%!error id=Octave:expected-real intrinsicToWorld (imref3d, 1j, 2, 3)
+%!error id=Octave:expected-real intrinsicToWorld (imref3d, 1, 2j, 3)
+%!error id=Octave:expected-real intrinsicToWorld (imref3d, 1, j, 3j)
+%!error id=Octave:invalid-input-arg intrinsicToWorld (imref3d, [1, 2], 3, 4)
+%!error id=Octave:invalid-input-arg intrinsicToWorld (imref3d, 1, [2, 3], 4)
+%!error id=Octave:invalid-input-arg intrinsicToWorld (imref3d, 1, 2, [3, 4])
+
+%!test
+%! r = imref3d ([128, 128, 27], 2, 2, 4);
+%! xI = [54, 71, 57, 70];
+%! yI = [46, 48, 79, 80];
+%! zI = [13, 13, 13, 13];
+%! [xW, yW, zW] = intrinsicToWorld (r, xI, yI, zI);
+%! assert (xW, [108, 142, 114, 140])
+%! assert (yW, [92, 96, 158, 160])
+%! assert (zW, [52, 52, 52, 52])
+
+%!test
+%! [xW, yW, zW] = intrinsicToWorld (imref3d, -5.3, -2.8, -15.88);
+%! assert (xW, -5.3)
+%! assert (yW, -2.8)
+%! assert (zW, -15.88, 1e-6)
+
+%!test
+%! [xW, yW, zW] = intrinsicToWorld (imref3d, [1, 2; 3, 4],
+%!                                           [2, 3; 5, 9],
+%!                                           [-5, 8; 19, 42.8]);
+%! assert (xW, [1, 2; 3, 4])
+%! assert (yW, [2, 3; 5, 9])
+%! assert (zW, [-5, 8; 19, 42.8])
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref3d/subsasgn.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref3d/subsasgn.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,101 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} {@var{rout} =} subsasgn (@var{r}, @var{index}, @var{val})
+##
+## @seealso{}
+## @end deftypefn
+
+function rout = subsasgn (r, index, val)
+  switch (index.type)
+    case "."
+      fld = index.subs;
+      switch (fld)
+        case "ImageSize"
+          imageSize = val;
+
+          if (length (imageSize) < 3)
+            error ("Octave:invalid-input-arg", ...
+            "ImageSize must have three elements");
+          endif
+
+          validateattributes (imageSize, {"numeric"}, ...
+          {"positive", "integer", "vector"}, "imref3d", "imageSize");
+          
+          m = imageSize(1);
+          n = imageSize(2);
+          p = imageSize(3);
+          
+          rout = r;
+          rout.ImageSize = imageSize;
+          rout.PixelExtentInWorldX = r.ImageExtentInWorldX / n;
+          rout.PixelExtentInWorldY = r.ImageExtentInWorldY / m;
+          rout.PixelExtentInWorldZ = r.ImageExtentInWorldZ / p;
+          rout.XIntrinsicLimits = [0.5, n + 0.5];
+          rout.YIntrinsicLimits = [0.5, m + 0.5];
+          rout.ZIntrinsicLimits = [0.5, p + 0.5];
+        case "XWorldLimits"
+          xWorldLimits = val;
+          
+          validateattributes (xWorldLimits, {"numeric"}, ...
+          {"increasing", "real", "vector", "size", [1, 2]}, ...
+          "imref3d", "xWorldLimits");
+          
+          imageSize = r.ImageSize;
+          imageExtentInWorldX = xWorldLimits(2) - xWorldLimits(1);
+          
+          rout = r;
+          rout.XWorldLimits = val;
+          rout.ImageExtentInWorldX = imageExtentInWorldX;
+          rout.PixelExtentInWorldX = imageExtentInWorldX / imageSize(2);
+        case "YWorldLimits"
+          yWorldLimits = val;
+          
+          validateattributes (yWorldLimits, {"numeric"}, ...
+          {"increasing", "real", "vector", "size", [1, 2]}, ...
+          "imref3d", "yWorldLimits");
+          
+          imageSize = r.ImageSize;
+          imageExtentInWorldY = yWorldLimits(2) - yWorldLimits(1);
+          
+          rout = r;
+          rout.YWorldLimits = val;
+          rout.ImageExtentInWorldY = imageExtentInWorldY;
+          rout.PixelExtentInWorldY = imageExtentInWorldY / imageSize(1);
+        case "ZWorldLimits"
+          zWorldLimits = val;
+          
+          validateattributes (zWorldLimits, {"numeric"}, ...
+          {"increasing", "real", "vector", "size", [1, 2]}, ...
+          "imref3d", "zWorldLimits");
+          
+          imageSize = r.ImageSize;
+          imageExtentInWorldZ = zWorldLimits(2) - zWorldLimits(1);
+          
+          rout = r;
+          rout.ZWorldLimits = val;
+          rout.ImageExtentInWorldZ = imageExtentInWorldZ;
+          rout.PixelExtentInWorldZ = imageExtentInWorldZ / imageSize(3);
+        otherwise
+          error ("Octave:invalid-indexing", ...
+          "@imref3d/subsasgn: invalid property '%s'", fld);
+      endswitch
+
+    otherwise
+      error ("Octave:invalid-indexing", "@imref3d/subsasgn: invalid index type")
+  endswitch
+endfunction
\ No newline at end of file
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref3d/subsref.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref3d/subsref.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,57 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} {@var{r} =} subref (@var{val}, @var{idx})
+##
+## @seealso{}
+## @end deftypefn
+
+function r = subsref (val, idx)
+  if (strcmp (idx.type, "."))
+    switch (idx.subs)
+      case "ImageSize"
+        r = val.ImageSize;
+      case "XWorldLimits"
+        r = val.XWorldLimits;
+      case "YWorldLimits"
+        r = val.YWorldLimits;
+      case "ZWorldLimits"
+        r = val.ZWorldLimits;
+      case "PixelExtentInWorldX"
+        r = val.PixelExtentInWorldX;
+      case "PixelExtentInWorldY"
+        r = val.PixelExtentInWorldY;
+      case "PixelExtentInWorldZ"
+        r = val.PixelExtentInWorldZ;
+      case "ImageExtentInWorldX"
+        r = val.ImageExtentInWorldX;
+      case "ImageExtentInWorldY"
+        r = val.ImageExtentInWorldY;
+      case "ImageExtentInWorldZ"
+        r = val.ImageExtentInWorldZ;
+      case "XIntrinsicLimits"
+        r = val.XIntrinsicLimits;
+      case "YIntrinsicLimits"
+        r = val.YIntrinsicLimits;
+      case "ZIntrinsicLimits"
+        r = val.ZIntrinsicLimits;
+      otherwise
+        error ("Octave:invalid-indexing", ...
+        strcat ("unknown property '", idx.subs, "' for class imref3d"));
+    endswitch
+  endif
+endfunction
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref3d/worldToIntrinsic.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref3d/worldToIntrinsic.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,81 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} {[@var{xIntrinsic}, @var{yIntrinsic}, @var{zIntrinsic}] =} worldToIntrinsic (@var{r}, @var{xWorld}, @var{yWorld}, @var{zWorld})
+## Convert from world to intrinsic coordinates.
+##
+## Converts world coordinates @var{xWorld}, @var{yWorld} and @var{zWorld} to
+## intrinsic coordinates @var{xIntrinsic}, @var{yIntrinsic} and @var{zIntrinsic}
+## of an image associated with the spatial referencing object @var{r}.  If a
+## point (@var{xWorld}(i), @var{yWorld}(i), @var{zWorld}(i)) falls outside
+## the bounds of the image, its intrinsic coordinates are extrapolated,
+## possibly resulting in negative values.
+##
+## @seealso{imref2d, imref3d, intrinsicToWorld}
+## @end deftypefn
+
+function [xIntrinsic, yIntrinsic, zIntrinsic] = worldToIntrinsic (r, xWorld, yWorld, zWorld)
+  if (nargin != 4)
+    print_usage ();
+  endif
+  
+  validateattributes (xWorld, {"numeric"}, ...
+  {"real"}, "imref3d", "xWorld");
+  validateattributes (yWorld, {"numeric"}, ...
+  {"real"}, "imref3d", "yWorld");
+  validateattributes (zWorld, {"numeric"}, ...
+  {"real"}, "imref3d", "zWorld");
+  
+  if (! all (size (xWorld) == size (yWorld)) ...
+    || ! all (size (xWorld) == size (zWorld)))
+    error ("Octave:invalid-input-arg", ...
+    "xWorld, yWorld and zWorld must be of the same size");
+  endif
+  
+  xIntrinsicLimits = r.XIntrinsicLimits;
+  yIntrinsicLimits = r.YIntrinsicLimits;
+  zIntrinsicLimits = r.ZIntrinsicLimits;
+  xWorldLimits = r.XWorldLimits;
+  yWorldLimits = r.YWorldLimits;
+  zWorldLimits = r.ZWorldLimits;
+  xIntrinsic = xIntrinsicLimits(1) + (xWorld  - xWorldLimits(1)) ...
+  / r.PixelExtentInWorldX;
+  yIntrinsic = yIntrinsicLimits(1) + (yWorld  - yWorldLimits(1)) ...
+  / r.PixelExtentInWorldY;
+  zIntrinsic = zIntrinsicLimits(1) + (zWorld  - zWorldLimits(1)) ...
+  / r.PixelExtentInWorldZ;
+endfunction
+
+%!error id=Octave:invalid-fun-call worldToIntrinsic (imref3d)
+%!error id=Octave:invalid-fun-call worldToIntrinsic (imref3d, 1, 2)
+%!error id=Octave:invalid-fun-call worldToIntrinsic (imref3d, 1, 2, 3, 4)
+%!error id=Octave:expected-real worldToIntrinsic (imref3d, 1j, 2, 3)
+%!error id=Octave:expected-real worldToIntrinsic (imref3d, 1, 2j, 3)
+%!error id=Octave:expected-real worldToIntrinsic (imref3d, 1, 2, 3j)
+%!error id=Octave:invalid-input-arg worldToIntrinsic (imref3d, [1, 2], 3, 4)
+%!error id=Octave:invalid-input-arg worldToIntrinsic (imref3d, 1, [2, 3], 4)
+%!error id=Octave:invalid-input-arg worldToIntrinsic (imref3d, 1, 2, [3, 4])
+
+%!test
+%! r = imref3d ([128, 128, 27], 2, 2, 4);
+%! xW = [108, 108, 108.2, 2];
+%! yW = [92, 92, 92, -1];
+%! zW = [52, 55, 52, 0.33];
+%! [xI, yI, zI] = worldToIntrinsic (r, xW, yW, zW);
+%! assert (xI, [54, 54, 54.1, 1], 1e-6)
+%! assert (yI, [46, 46, 46, -0.5], 1e-6)
+%! assert (zI, [13, 13.75, 13, 0.0825], 1e-6)
\ No newline at end of file
diff -r 3b1f25568ccd -r 76d1155e62fd inst/@imref3d/worldToSubscript.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/@imref3d/worldToSubscript.m	Thu Nov 15 01:00:09 2018 +0100
@@ -0,0 +1,82 @@
+## Copyright (C) 2018 Martin Janda <janda.martin1@gmail.com>
+## 
+## This program is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program.  If not, see
+## <https://www.gnu.org/licenses/>.
+
+## -*- texinfo -*- 
+## @deftypefn {} {[@var{i}, @var{j}, @var{k}] =} worldToSubscript (@var{r}, @var{xWorld}, @var{yWorld}, @var{zWorld})
+## Convert world coordinates to row, column and plane subscripts.
+##
+## Converts world coordinates to row, column and plane subscripts of an image
+## associated with the spatial referencing object @var{r}.  A point located at
+## (@var{xWorld}(i), @var{yWorld}(i), @var{zWorld}(i)) world coordinates maps
+## to row, column and plane subscripts @var{i}(i), @var{j}(i), @var{k}(i)
+## respectively.  Note the reversed order of the first two dimensions.  If the
+## point falls outside the bounds of the image, all of its subscripts are NaN.
+##
+## @seealso{imref2d, imref3d, worldToIntrinsic}
+## @end deftypefn
+
+function [i, j, k] = worldToSubscript (r, xWorld, yWorld, zWorld)
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  validateattributes (xWorld, {"numeric"}, ...
+  {"real"}, "imref3d", "xWorld");
+  validateattributes (yWorld, {"numeric"}, ...
+  {"real"}, "imref3d", "yWorld");
+  validateattributes (zWorld, {"numeric"}, ...
+  {"real"}, "imref3d", "zWorld");
+  
+  if (! all (size (xWorld) == size (yWorld)) ...
+    || ! all (size (xWorld) == size (zWorld)))
+    error ("Octave:invalid-input-arg", ...
+    "xWorld, yWorld and zWorld must be of the same size");
+  endif
+
+  [xIntrinsic, yIntrinsic, zIntrinsic] ...
+  = worldToIntrinsic (r, xWorld, yWorld, zWorld);
+  xIntrinsicLimits = r.XIntrinsicLimits;
+  yIntrinsicLimits = r.YIntrinsicLimits;
+  zIntrinsicLimits = r.ZIntrinsicLimits;
+  inImage = contains (r, xWorld, yWorld, zWorld);
+  xIntrinsic(! inImage) = NaN;
+  yIntrinsic(! inImage) = NaN;
+  zIntrinsic(! inImage) = NaN;
+  i = round (yIntrinsic);
+  j = round (xIntrinsic);
+  k = round (zIntrinsic);
+endfunction
+
+%!error id=Octave:invalid-fun-call worldToSubscript (imref3d)
+%!error id=Octave:invalid-fun-call worldToSubscript (imref3d, 1)
+%!error id=Octave:invalid-fun-call worldToSubscript (imref3d, 1, 2)
+%!error id=Octave:invalid-fun-call worldToSubscript (imref3d, 1, 2, 3, 4)
+%!error id=Octave:expected-real worldToSubscript (imref3d, 1j, 2, 3)
+%!error id=Octave:expected-real worldToSubscript (imref3d, 1, 2j, 3)
+%!error id=Octave:expected-real worldToSubscript (imref3d, 1, 2, 3j)
+%!error id=Octave:invalid-input-arg worldToSubscript (imref3d, [1, 2], 3, 4)
+%!error id=Octave:invalid-input-arg worldToSubscript (imref3d, 1, [2, 3], 4)
+%!error id=Octave:invalid-input-arg worldToSubscript (imref3d, 1, 2, [3, 4])
+
+%!test
+%! r = imref3d ([128, 128, 27], 2, 2, 4);
+%! xW = [108, 108, 113.2, 2];
+%! yW = [92, 92, 92, -1];
+%! zW = [52, 55, 52, 0.33];
+%! [rS, cS, pS] = worldToSubscript (r, xW, yW, zW);
+%! assert (rS, [46, 46, 46, NaN])
+%! assert (cS, [54, 54, 57, NaN])
+%! assert (pS, [13, 14, 13, NaN])
\ No newline at end of file
