--- a/INDEX
+++ b/INDEX
@@ -24,6 +24,7 @@
  imadd
  imcomplement
  imdivide
+ imlincomb
  immultiply
  imsubtract
 Black and white image functions

--- a/NEWS
+++ b/NEWS
@@ -10,6 +10,7 @@
       imadd
       imbothat
       imdivide
+      imlincomb
       immultiply
       imsubtract
       iptcheckconn

--- a/inst/imabsdiff.m
+++ b/inst/imabsdiff.m
@@ -31,7 +31,7 @@
 ##
 ## @emph{Note 2}: the values are truncated to the mininum value of the output
 ## class.
-## @seealso{imadd, imcomplement, imdivide, immultiply, imsubtract}
+## @seealso{imadd, imcomplement, imdivide, imlincomb, immultiply, imsubtract}
 ## @end deftypefn
 
 function img = imabsdiff (img, val, out_class = class (img))

--- a/inst/imadd.m
+++ b/inst/imadd.m
@@ -32,7 +32,7 @@
 ##
 ## @emph{Note 2}: the values are truncated to the maximum value of the output
 ## class.
-## @seealso{imabsdiff, imcomplement, imdivide, immultiply, imsubtract}
+## @seealso{imabsdiff, imcomplement, imdivide, imlincomb, immultiply, imsubtract}
 ## @end deftypefn
 
 function img = imadd (img, val, out_class = class (img))

--- a/inst/imcomplement.m
+++ b/inst/imcomplement.m
@@ -21,7 +21,7 @@
 ## For binary images, the complement is computed as @code{!@var{A}}, for floating
 ## point images it is computed as @code{1 - @var{A}}, and for integer images as
 ## @code{intmax(class(@var{A})) - @var{A}}.
-## @seealso{imadd, imdivide, immultiply, imsubtract}
+## @seealso{imadd, imdivide, imlincomb, immultiply, imsubtract}
 ## @end deftypefn
 
 function B = imcomplement(A)

--- a/inst/imdivide.m
+++ b/inst/imdivide.m
@@ -28,7 +28,7 @@
 ##
 ## @emph{Note}: the values are truncated to the mininum value of the output
 ## class.
-## @seealso{imabsdiff, imadd, imcomplement, immultiply, imsubtract}
+## @seealso{imabsdiff, imadd, imcomplement, immultiply, imlincomb, imsubtract}
 ## @end deftypefn
 
 function img = imdivide (img, val, out_class = class (img))

new file mode 100644
--- /dev/null
+++ b/inst/imlincomb.m
@@ -0,0 +1,90 @@
+## Copyright (C) 2011 Carnë Draug <carandraug+dev@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 <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{out} =} imlincomb (@var{fac}, @var{img})
+## @deftypefnx {Function File} {@var{out} =} imlincomb (@var{fac1}, @var{img1}, @var{fac2}, @var{img2}, @dots{})
+## @deftypefnx {Function File} {@var{out} =} imlincomb (@dots{}, @var{class})
+## Combine images linearly.
+##
+## Returns the computed image as per:
+##
+## @var{out} = @var{fac}1*@var{img}1 + @var{fac}2*@var{img}2 + @dots{} + @var{fac}n*@var{img}n
+##
+## The images @var{img}1..n must all be of same size and class. The factors @var{fac}1..n
+## must all be floating-point scalars.
+##
+## The class of @var{out} will be the same as @var{img}s unless @var{img}s are logical
+## in which case @var{out} will be double. Alternatively, it can be specified
+## with @var{class}.
+##
+## If applying several arithmetic operations on images, @code{imlincomb} is more
+## precise since calculations are performed at double precision.
+##
+## @emph{Note 1}: you can force output class to be logical by specifying
+## @var{class} though it possibly doesn't make a lot of sense.
+##
+## @seealso{imadd, imcomplement, imdivide, immultiply, imsubtract}
+## @end deftypefn
+
+function out = imlincomb (varargin)
+
+  if (nargin < 2)
+    print_usage;
+  else
+    if (rem (nargin, 2) == 0)
+      ## use default for output class; the class of first image (second argument)
+      out_class = class (varargin{2});
+      def_class = true;
+    else
+      ## last argument is requested output class
+      out_class = varargin{end};
+      def_class = false;
+    endif
+  endif
+
+  facI = 1:2:nargin-1;        # index for factors
+  imgI = 2:2:nargin;          # index for images
+  imgC = class (varargin{2}); # class of the first image
+
+  out = zeros (size (varargin{2}));
+  for i = 1:numel (imgI)
+    ## we keep index the images froom varargin rather than copying to new variable to
+    ## avoid taking up a lot of memory
+    if (!isreal (varargin{facI(i)}) || !isscalar (varargin{facI(i)}) || !isfloat (varargin{facI(i)}))
+      error ("factor to multiply each image must be a real, floating-point, scalar.");
+    elseif ((!isnumeric (varargin{imgI(i)}) && !islogical (varargin{imgI(i)})) ...
+                || isempty (varargin{imgI(i)}) || issparse (varargin{imgI(i)}) ...
+                || !isreal (varargin{imgI(i)}) || !isa (varargin{imgI(i)}, imgC))
+      error ("images must be a numeric or logical, non-empty, non-sparse real matrix of the same class.");
+    endif
+    img  = double(varargin{imgI(i)}) .* double(varargin{facI(i)});
+    out += img;
+  endfor
+
+  ## this is probably matlab imcompatible since by their documentation, they don't even
+  ## support logical matrix. If specified by user, respect and return a logical
+  ## matrix. Otherwise, return a double, even if images were all logical
+  if (strcmpi (out_class, "logical") && def_class)
+    cnv = @double;
+  else
+    cnv = str2func (tolower (out_class));
+  endif
+  out = cnv (out);
+
+endfunction
+
+%!assert (imlincomb (0.5, uint8 ([255 10]), 0.5, uint8 ([50 20])),           uint8  ([153 15])); # default to first class and truncate
+%!assert (imlincomb (0.5, uint8 ([255 10]), 0.5, uint8 ([50 20]), "uint16"), uint16 ([153 15])); # defining output class works

--- a/inst/immultiply.m
+++ b/inst/immultiply.m
@@ -28,7 +28,7 @@
 ##
 ## @emph{Note}: the values are truncated to the mininum value of the output
 ## class.
-## @seealso{imabsdiff, imadd, imcomplement, imdivide, imsubtract}
+## @seealso{imabsdiff, imadd, imcomplement, imdivide, imlincomb, imsubtract}
 ## @end deftypefn
 
 function img = immultiply (img, val, out_class = class (img))

--- a/inst/imsubtract.m
+++ b/inst/imsubtract.m
@@ -46,7 +46,7 @@
 ##
 ## Because both 190 and 200 were truncated to 127 before subtraction, their difference
 ## is zero.
-## @seealso{imabsdiff, imadd, imcomplement, imdivide, immultiply}
+## @seealso{imabsdiff, imadd, imcomplement, imdivide, imlincomb, immultiply}
 ## @end deftypefn
 
 function img = imsubtract (img, val, out_class = class (img))