/*
Copyright (c) 2019- Ferhat Kurtulmuş
Boost Software License - Version 1.0 - August 17th, 2003
Permission is hereby granted, free of charge, to any person or organization
obtaining a copy of the software and accompanying documentation covered by
this license (the "Software") to use, reproduce, display, distribute,
execute, and transmit the Software, and to prepare derivative works of the
Software, and to permit third-parties to whom the Software is furnished to
do so, all subject to the following:
The copyright notices in the Software and this entire statement, including
the above license grant, this restriction and the following disclaimer,
must be included in all copies of the Software, in whole or in part, and
all derivative works of the Software, unless such copies or derivative
works are solely in the form of machine-executable object code generated by
a source language processor.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.
*/

module measure.types;

import std.math;
import std.algorithm;
import std.algorithm.searching;

struct XYList {
    // this is useful for many situation instead of using Point[]
    int[] xs;
    int[] ys;
}

struct Rectangle{
    int x;
    int y;
    int width;
    int height;
}

struct Ellipse{
    double angle;
    double center_x; 
    double center_y;
    double maj;
    double min;
}

struct Point {
    int x;
    int y;
 
    int opCmp(Point rhs) {
        if (x < rhs.x) return -1;
        if (rhs.x < x) return 1;
        return 0;
    }
 
    void toString(scope void delegate(const(char)[]) sink) const {
        import std.format;
        sink("(");
        formattedWrite(sink, "%d", x);
        sink(",");
        formattedWrite(sink, "%d", y);
        sink(")");
    }
}

class Region{
    Mat2D!ubyte image;
    
    // spatial raw moments
    double m00, m10, m01, m20, m11, m02, m30, m21, m12, m03;
    
    ulong area;
    double areaFromContour;
    double perimeter;
    Point centroid;
    double aspect_Ratio;
    Rectangle bBox;
    XYList convexHull;
    double convexArea;
    Ellipse ellipse;
    double extent;
    double solidity;
    double majorAxisLength;
    double minorAxisLength;
    double orientation;
    double eccentricity;
    double equivalentDiameter;
    XYList contourPixelList; // chain sorted!
    XYList pixelList;
    
    this(){}
}

struct Mat2D(T)
{
    T[] data;
    size_t rows;
    size_t cols;
    
    alias height = rows;
    alias width = cols;
    
    this(T)(T[] array, size_t r, size_t c)
    {
        data = array;
        rows = r;
        cols = c;
        
        assert (data.length == rows* cols,
            "array length does not match with rows * cols!");
    }
    
    this(size_t r, size_t c)
    {
        data = new T[r*c];
        rows = r;
        cols = c;
        data[0..$] = 0;
    }
    
    this(size_t n){
        data = new T[n*n];
        rows = n;
        cols = n;
        data[0..$] = 0;
    }
    
    T opIndex( size_t i, size_t j)
    {
        assert((i < rows) && (j < cols), "index out of bounds!");
        return data[i * cols + j];
    }
    
    void opIndexAssign(T val, size_t i, size_t j)
    {
        assert((i < rows) && (j < cols), "index out of bounds!");
        data[i * cols + j] = val;
    }
    
    bool opEquals(Mat2D!T other)
    {
        return data == other.data;
    }
    
    T max()
    {
        return data.maxElement;
    }
    
    T min()
    {
        return data.minElement;
    }
    
    int[2] argMax()
    {
        const size_t flatInd = maxIndex(data);
        int yy = cast(int)(flatInd % cols);
        int xx = cast(int)(flatInd / cols);
        return [xx, yy];
    }
    int[2] argMin()
    {
        const size_t flatInd = minIndex(data);
        int yy = cast(int)(flatInd % cols);
        int xx = cast(int)(flatInd / cols);
        return [xx, yy];
    }
    
    Mat2D!T opMul(S)(Mat2D!S b)
    {
        const size_t r1 = rows;
        const size_t c1 = cols;

        const size_t r2 = b.rows;
        const size_t c2 = b.cols;
        
        assert ((c1 == r2), "column of first matrix in not equal to row of second matrix"); 
        
        Mat2D!T mult = Mat2D!T(r1, c2);
        
        foreach(i; 0..r1)
            foreach(j; 0..c2){
                T tempSum = 0;
                foreach(k; 0..c1)
                    tempSum += this[i, k] * b[k, j];
                mult[i, j] = tempSum;
            }
                
        
        return mult;
    }
    
    Mat2D!T opMul(T scalar)
    {
        Mat2D!T mult = Mat2D!T(rows, cols);
        foreach(i, val; data) mult.data[i] = cast(T)(val * scalar);
        return mult;
    }
    
    Mat2D!T opMulAssign (Mat2D!T mat)
    {
        this = this * mat;
        return this;
    }
    
    Mat2D!T opMulAssign(T scalar)
    {
        this = this * scalar;
        return this;
    }
    
    static Mat2D!T diagFromArray(T)(T[] v){
        size_t len = v.length;
        Mat2D!T d = Mat2D!T(len, len);
        
        foreach(y; 0..len)
            foreach(x; 0..len){
                    d[y, x] = 0;
                    if (y == x){
                        d[y, x] = v[x];
                    }
                    else
                        d[y, x] = 0;
            }
        
        return d;
    }
    
    T[] diag()
    {
        assert(rows == cols, "The matrix is not square! Diagonal cannot be created!");
        
        T[] d; d.length = cols;
        
        foreach(y; 0..rows)
            foreach(x; 0..rows){
                if (y == x) d[x] = this[y, x];
            }
        return d;
    }
    
    Mat2D!T transpose()
    {

        Mat2D!T transp = Mat2D!T(cols, rows);

        foreach(i; 0..rows)
            foreach(j; 0..cols)
                transp[j, i] = this[i, j];
        return transp;
    }
    
    Mat2D!double inverse() 
    { 
        // https://www.geeksforgeeks.org/adjoint-inverse-matrix/
        
        assert(rows == cols, "The matrix is not square, cannot inverse!");
        
        size_t N = rows;
        Mat2D!double inv = Mat2D!double(rows, rows);

        const double det = determinant(this, N);
        assert (det != 0); //  Singular matrix, can't find its inverse
      
        // Find adjoint 
        Mat2D!double adj = Mat2D!double(rows, rows);
        adjoint(this, adj);
        
        for (int i=0; i<N; i++) 
            for (int j=0; j<N; j++) 
                inv[i, j] = adj[i, j]/cast(double)det; 
      
        return inv; 
    }
}

/+ planned for performance, but it badly affects the speed for some reason?
struct ROIViewOfMat2D(T)
{
    T* data;
    size_t rows;
    size_t cols;
    
    alias height = rows;
    alias width = cols;
    
    Rectangle roi;
    
    this(T)(Mat2D!T parentMat, Rectangle _roi)
    {
        data = parentMat.data.ptr;
        rows = parentMat.rows;
        cols = parentMat.cols;
        roi = _roi;
    }
    
    T opIndex( size_t _i, size_t _j)
    {
        size_t i = _i + roi.y;
        size_t j = _j + roi.x;
        
        assert((i < rows) && (j < cols), "index out of bounds!");
        return data[i * cols + j];
    }
}
+/

void getCofactor(T)(Mat2D!T A, Mat2D!double temp, int p, int q, ulong n) { 
    int i = 0, j = 0; 
  
    for (int row = 0; row < n; row++) 
    { 
        for (int col = 0; col < n; col++) 
        {  
            if (row != p && col != q) 
            { 
                temp[i, j++] = A[row, col]; 
  
                if (j == n - 1) 
                { 
                    j = 0; 
                    i++; 
                } 
            } 
        } 
    } 
} 

double determinant(T)(Mat2D!T A, ulong n)
{ 
    size_t N = A.rows;
    double D = 0;
    
    if (n == 1) 
        return A[0, 0]; 
  
    Mat2D!double temp = Mat2D!double(N, N);
  
    int sign = 1;
    
    for (int f = 0; f < n; f++) 
    {
        getCofactor(A, temp, 0, f, n); 
        D += sign * A[0, f] * determinant(temp, n - 1); 
        sign = -sign; 
    } 
  
    return D; 
}

void adjoint(T)(Mat2D!T A, Mat2D!double adj) 
{ 
    size_t N = A.rows;
    
    if (N == 1) 
    { 
        adj[0, 0] = 1; 
        return; 
    } 
  
    int sign = 1;
    Mat2D!double temp = Mat2D!double(N, N); 
  
    for (int i=0; i<N; i++) 
    { 
        for (int j=0; j<N; j++) 
        {
            getCofactor(A, temp, i, j, N); 
            sign = ((i+j)%2==0)? 1: -1; 
            adj[j, i] = (sign)*(determinant(temp, N-1)); 
        } 
    } 
}

unittest
{
    double[] aa = [1,2,3,4,5,6];
    Mat2D!double maa = Mat2D!double(aa, 2, 3);
    
    double[] bb = [7,8,9,10,11,12];
    Mat2D!double mbb = Mat2D!double(bb, 3, 2);
    
    auto mul = maa*mbb;
    assert(maa*mbb == Mat2D!double([58.0, 64.0, 139.0, 154.0], 2, 2));
    
    Mat2D!int cc = Mat2D!int([1,2,3,0,1,5,5,6,0], 3, 3);
    assert(cc.inverse().data == [-6.0, 3.6, 1.4, 5.0, -3.0, -1.0, -1.0, 0.8, 0.2] );
}