morph.cpp

#include <Halide.h>
using namespace Halide;

#include "../png.h"
#include <typeinfo>
#include <iostream>
#include <string>
#include <math.h>

#define PI 3.141592653589

Expr width;
Expr height;
Var x("x"),y("y"),c("c");

Expr twoSegX(UniformImage segBefore, UniformImage segAfter,  Var x, Var y){
    //RVar i(0,1);
    Expr i =0;
    Expr PX0 = x-segAfter(i,0);
    Expr PX1 = y-segAfter(i,1);
    Expr PQ0 = segAfter(i, 2)-segAfter(i,0);
    Expr PQ1 = segAfter(i,3)-segAfter(i,1);
    Expr length =sqrt(PQ0*PQ0+PQ1*PQ1);
    Expr u = (PX0*PQ0+PX1*PQ1)/cast<float>(PQ0*PQ0+PQ1*PQ1);
    Expr PQperp0 = -segAfter(i,3)+segAfter(i,1);
    Expr PQperp1 = segAfter(i,2) - segAfter(i,0);
    Expr v = (PX0*PQperp0+PX1*PQperp1)/length;
    Expr PQ20 = segBefore(i,2)-segBefore(i,0);
    Expr PQ21 = segBefore(i,3)-segBefore(i,1);
    Expr length2 =sqrt(PQ20*PQ20+PQ21*PQ21);
    Expr PQ2perp0 = -segBefore(i,3)+segBefore(i,1);
    Expr PQ2perp1 = segBefore(i,2) - segBefore(i,0);

    Expr newx = segBefore(i,0) + u*PQ20 + v*PQ2perp0/length2;
    return newx;
}

Expr twoSegY(Image<float> segBefore, Image<float> segAfter, Expr index, Var x, Var y){
    Func PX, PQ, PQperp, PQ2perp, PQ2;
    PX = (x-segAfter(index,0),y-segAfter(index,1));
    PQ = (segAfter(index, 2)-segAfter(index,0), segAfter(index,3)-segAfter(index,1));
    PQperp = (-segAfter(index,3)+segAfter(index,1), segAfter(index,2) - segAfter(index,0));
    PQ2 = (segBefore(index,2)-segBefore(index,0), segBefore(index,3)-segBefore(index,1));
    PQ2perp = (-segBefore(index,3)+segBefore(index,1), segBefore(index,2) - segBefore(index,0));
    Expr bla = PQ(0)*PQ(0)+PQ(1)*PQ(1);
    Expr length =sqrt(bla);
    Expr u = (PX(0)*PQ(0)+PX(1)*PQ(1))/cast<float>(bla);
    Expr v = (PX(0)*PQperp(0)+PX(1)*PQperp(1))/length;
    Expr length2 =sqrt(PQ2(0)*PQ2(0)+PQ2(1)*PQ2(1));

    Expr newy = segBefore(index,1) + v*PQ2perp(1)/length2 + u*PQ2(1);
    return cast<int>(newy);
}

Func warpBy1(Func f, Image<float> segBefore, Image<float> segAfter){
    Func X2;
    Func out;
    X2(x,y) = (twoSegY(segBefore,segAfter,0,x,y),twoSegY(segBefore,segAfter,0,x,y));
    out(x,y) =f(X2(x,y,0), X2(x,y,1));
    return out;
    }

Func warp(Func f, Func segBefore, Func segAfter){
    RVar i(0,2),j(0,1);
    //Expr segBefore = listOfSegBefore[i];
    //Expr segAfter = listOfSegAfter[i];
    //float* segAfter = listOfSegAfter[0];
    //float *segBefore = listOfSegBefore[0];
    // Func lala;
    // // lala(x) = segAfter[cast<int>(x)];
    Expr PX0 = x-segAfter(0,0);
    Expr PX1 = y-segAfter(0,1);
    Expr PQ0 = segAfter(0,2)-segAfter(0,0);
    Expr PQ1 = segAfter(0,3)-segAfter(0,1);
    Expr length =sqrt(PQ0*PQ0+PQ1*PQ1);
    Expr u = (PX0*PQ0+PX1*PQ1)/cast<float>(PQ0*PQ0+PQ1*PQ1);
    Expr PQperp0 = -segAfter(0,3)+segAfter(0,1);
    // Expr PQperp1 = segAfter[0][2] - segAfter[0][0];
    // Expr v = (PX0*PQperp0+PX1*PQperp1)/length;
    // Expr PQ20 = segBefore[0][2]-segBefore[0][0];
    // Expr PQ21 = segBefore[0][3]-segBefore[0][1];
    // Expr length2 =sqrt(PQ20*PQ20+PQ21*PQ21);
    // Expr PQ2perp0 = -segBefore[0][3]+segBefore[0][1];
    // Expr PQ2perp1 = segBefore[0][2] - segBefore[0][0];

    // Expr newx = segBefore[0][0] + u*PQ20 + v*PQ2perp0/length2;
    // Expr newy = segBefore[0][1] + u*PQ21 + v*PQ2perp1/length2;

    // Expr bla = i+segBefore[0][0];
    // Func out;
    // out(x,y) = f(x+sum(bla), cast<int>(newy));

    //krenemo od pozadi
    //out(x,y) = f(x+sum(/weightsum
    return f;


}


int main(int argc, char **argv) {

    /* THE ALGORITHM */

    UniformImage input(UInt(16),2);
    UniformImage input2(UInt(16),2);
    UniformImage sobel(Float(32), 2);
    Image<float> segBefore = {{71, 93, 89, 93},
                            {113, 92, 132, 90}};
    Image<float> segAfter = {{68, 100, 89, 98},
                           {115, 97, 134, 96}};
    Func sobel_float;
    sobel_float(x,y,c) = cast<float>(sobel(x,y,c));
    UniformImage segvBefore(UInt(16), 2);
    UniformImage segvAfter(UInt(16), 2);

    width = input.width();
    height = input.height();
    Func output("output");
    Var i,j;

    // Convert to floating point
    Func floating("floating");
    floating(x, y, c) = cast<float>(input(x, y, c)) / 65535.0f;
    Func floating2("floating2");
    floating2(x, y, c) = cast<float>(input2(x, y, c)) / 65535.0f;
    floating.root();

    // Set a boundary condition
    Func clamped;
    clamped(x, y, c) = floating(clamp(x, 0, width-1), clamp(y, 0, height-1), c);
    Func clamped2;
    clamped2(x, y, c) = floating(clamp(x, 0, width-1), clamp(y, 0, height-1), c);

    //test twosegments
    //float seg1[2][4] = {{71, 93, 89, 93},{71, 93, 89, 93}};
    //float seg2[4] = {68, 100, 89, 98};
    //Func warped = warpBy1(clamped, segmentsBefore, segmentsAfter, sobel_float);
    //Func warped_clamped;
    //warped_clamped(x,y,c) = warped(clamp(x, 0, width-1), clamp(y, 0, height-1), c);
    clamped.root();

    Func warped = warpBy1(clamped,segBefore, segAfter);
    warped.root();
    output(x,y) = cast<uint16_t>(clamp(warped(x,y), 0.0f, 1.0f) * 65535.0f);

    output.compileToFile("local_laplacian", {input, input2, segvBefore, segvAfter});


    return 0;
}