scripts/builtin/img_transform.dml

#-------------------------------------------------------------
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.
#
#-------------------------------------------------------------

# The Image Transform function applies an affine transformation to an image.
# Optionally resizes the image (without scaling).
# Uses nearest neighbor sampling.
#
# .. code-block:: python
#
#   >>> import numpy as np
#   >>> from systemds.context import SystemDSContext
#   >>> from systemds.operator.algorithm import img_transform
#   >>> 
#   >>> with SystemDSContext() as sds:
#   ...     img = sds.from_numpy(
#   ...         np.array([[ 10., 20., 30.],
#   ...                   [ 40., 50., 60.],
#   ...                   [ 70., 80., 90.]], dtype=np.float32)
#   ...     )
#   ...     result_img = img_transform(img, 3, 3, -1., 0., 2., 0., 1., 0., 255.).compute()
#   ...     print(result_img)
#   [[ 20.  10. 255.]
#    [ 50.  40. 255.]
#    [ 80.  70. 255.]]
#
#
# INPUT:
# -------------------------------------------------------------------------------------------
# img_in       Input image as 2D matrix with top left corner at [1, 1]
# out_w        Width of the output image
# out_h        Height of the output image
# a,b,c,d,e,f  The first two rows of the affine matrix in row-major order
# fill_value   The background of the image
# -------------------------------------------------------------------------------------------
#
# OUTPUT:
# ---------------------------------------------------------------------------------------
# img_out  Output image as 2D matrix with top left corner at [1, 1]
# ---------------------------------------------------------------------------------------

m_img_transform = function(Matrix[Double] img_in, Integer out_w, Integer out_h, Double a,
    Double b, Double c, Double d, Double e, Double f, Double fill_value) 
    return (Matrix[Double] img_out) {
  divisor = a * e - b * d
  if(divisor == 0) {
    print("Inverse matrix does not exist! Returning input.")
    img_out = img_in
  }
  else {
    orig_w = ncol(img_in)
    orig_h = nrow(img_in)
    # inverted transformation matrix
    # inversion is necessary because we compute the sampling position of pixels in the output image
    # and not the output coordinates of input pixels
    T_inv = matrix(0, rows=3, cols=3)
    T_inv[1, 1] = e / divisor
    T_inv[1, 2] = -b / divisor
    T_inv[1, 3] = (b * f - c * e) / divisor
    T_inv[2, 1] = -d / divisor
    T_inv[2, 2] = a / divisor
    T_inv[2, 3] = (c * d - a * f) / divisor
    T_inv[3, 3] = 1

    # coordinates of output pixel-centers linearized in row-major order
    coords = matrix(1, rows=3, cols=out_w*out_h)
    coords[1,] = t((seq(0, out_w*out_h-1) %% out_w) + 0.5)
    coords[2,] = t((seq(0, out_w*out_h-1) %/% out_w) + 0.5)

    # compute sampling pixel indices
    coords = floor(T_inv %*% coords) + 1

    img_out = matrix(fill_value, rows=out_h, cols=out_w)
    
    inx = t(coords[1,])
    iny = t(coords[2,])

    # any out-of-range pixels, if present, correspond to an extra pixel with fill_value at the end of the input
    index_vector = (orig_w *(iny-1) + inx) * ((0<inx) & (inx<=orig_w) & (0<iny) & (iny<=orig_h)) 
    index_vector = t(index_vector)
    xs = ((index_vector == 0)*(orig_w*orig_h +1)) + index_vector

    y = matrix(img_in, 1, orig_w*orig_h)
    if(min(index_vector) == 0){
      ys= cbind(y, matrix(fill_value,1, 1))
    }else{
      ys = y
    }

    ind= matrix(seq(1,ncol(xs),1),1,ncol(xs))
    z = table(xs, ind) 
    output = ys%*%z

    img_out = matrix(output, rows=out_h, cols=out_w)
  }
  
}