#-------------------------------------------------------------
#
# 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.
#
#-------------------------------------------------------------

# Image Transform
# Applies an affine transformation to an image.
# Optionally resizes the image (without scaling).
# Uses nearest neighbor sampling.
#
# INPUT PARAMETERS:
# ---------------------------------------------------------------------------------------------
# NAME       TYPE   DEFAULT MEANING
# ---------------------------------------------------------------------------------------------
# img_in     Matrix ---     Input image as 2D matrix with top left corner at [1, 1]
# out_w      Int    ---     Width of the output image
# out_h      Int    ---     Height of the output image
# abcdef     Int    ---     The first two rows of the affine matrix in row-major order
# fill_value Double ---     The background of the image

# RETURN VALUES
# ---------------------------------------------------------------------------------------------
# img_out    Matrix ---     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

    # fill output image
    img_out = matrix(fill_value, rows=out_w*out_h, cols=1)
    parfor (cell in 1:(out_w*out_h)) {
      inx = as.scalar(coords[1, cell])
      iny = as.scalar(coords[2, cell])
      if ((0 < inx) & (inx <= orig_w) & (0 < iny) & (iny <= orig_h))
        img_out[cell] = img_in[iny, inx]
    }

    # TODO replace above loop with following vectorized code
    # but additional size mismatch / fill handling necessary
    # ---
    # P = order(target=t(coords), by=matrix("1 2",1,2), index.return=TRUE);
    # inx = t(coords[1,]);
    # iny = t(coords[2,]);
    # vals = P %*% matrix(img_in, out_w*out_h, 1);
    # img_out = ((0<inx) & (inx<=orig_w) & (0<iny) & (iny<=orig_h)) * vals; 

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