image morphing in Halide

morph.cpp

#include <Halide.h>
using namespace Halide;

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

#define PI 3.14159

Expr width;
Expr height;

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

    Expr QX0 = x-segAfter(2, i);
    Expr QX1 = y-segAfter(3, i);

    answer[0] = segBefore(0, i) + u*PQ20 + v*PQ2perp0/length2;
    answer[1] = segBefore(1, i) + v*PQ2perp1/length2 + u*PQ21;
    answer[2] = select(cast<float>(u)<0.0f, sqrt(PX0*PX0+PX1*PX1), \
                           select(cast<float>(u)>1.0f, sqrt(QX0*QX0+QX1*QX1), \
                           select(cast<float>(v)<0.0f, -v, v)));
    answer[3] = length;

}

Func warpBy1(Func f, UniformImage segBefore, UniformImage segAfter){
    Var x,y;
    Func out;
    Expr la[4];
    twoSeg(segBefore, segAfter, 0, x, y, la);
    out(x,y) =f(la[0], la[1]);
    return out;
    }

// bi-linear interpolation
Expr interpolate(Func f, Expr x, Expr y){
    Expr newX = cast<int>(x);
    Expr newY = cast<int>(y);

    Expr weight1 = newX+1 - x;
    Expr weight2 = newY+1 - y;
    Expr inter = (f(newX, newY)*weight1 + f(newX+1, newY)*(1-weight1))*weight2 +\
        (f(newX, newY+1)*weight1 + f(newX+1, newY+1)*(1-weight1))*(1-weight2);
    return inter;
}

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


    Func X2, dist, len;
    X2(x,y,i) = (segBefore(0, i) + u*PQ20 + v*PQ2perp0/length2, \
                 segBefore(1, i) + v*PQ2perp1/length2 + u*PQ21);

    Expr distance = select(cast<float>(u)<0.0f, sqrt(PX0*PX0+PX1*PX1), \
                           select(cast<float>(u)>1.0f, sqrt(QX0*QX0+QX1*QX1), \
                           select(cast<float>(v)<0.0f, -v, v)));
    dist(x,y,i) = distance;
    len(x,y,i) = length;

    // go through all the segments
    RVar index(0, no_segments);
    Func weightsum;
    weightsum(x,y) = sum(len(x,y,index)/cast<float>(10+dist(x,y,index)));

    Func DSUM;
    DSUM(x,y) = (sum((X2(x,y,index,0)-x)*(len(x,y,index)/cast<float>(10+dist(x,y,index)))),\
                 sum((X2(x,y,index,1)-y)*(len(x,y,index)/cast<float>(10+dist(x,y,index)))));
    //calculate new x and y coordinates and interpolate
    Expr newX = x+DSUM(x,y,0)/weightsum(x,y);
    Expr newY = y+DSUM(x,y,1)/weightsum(x,y);
    out(x,y) = interpolate(f,newX, newY);
    return out;

}

Func morph(Func im1, Func im2, Func segBefore, Func segAfter, Expr no_segments, Uniform<float> time){
    Var x,y,c;
    Func intermediateSeg, out1, out2, out;
    intermediateSeg(x,y) = segBefore(x,y)*(1.0f-time) + segAfter(x,y)*time;
    out1 = warp(im1, segBefore, intermediateSeg, no_segments);
    out2 = warp(im2, segAfter, intermediateSeg, no_segments);
    out(x,y,c) = time*out2(x,y,c) + (1.0f-time)*out1(x,y,c);
    return out;
}


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

    Var x("x"),y("y"),c("c"), i, j;
    UniformImage inputStart(UInt(16),2);
    UniformImage inputEnd(UInt(16),2);
    UniformImage segBefore(Float(32), 2);
    UniformImage segAfter(Float(32), 2);
    Uniform<float> time;
    Uniform<int> no_segments;

    width = inputStart.width();
    height = inputStart.height();
    Func output("output");

    // Convert to floating point
    Func floatingStart("floatingStart");
    floatingStart(x, y, c) = cast<float>(inputStart(x, y, c)) / 65535.0f;
    Func floatingEnd("floatingEnd");
    floatingEnd(x, y, c) = cast<float>(inputEnd(x, y, c)) / 65535.0f;

    // Convert to functions
    Func segmentsBefore;
    segmentsBefore(x,y) = segBefore(x,y);
    Func segmentsAfter;
    segmentsAfter(x,y) = segAfter(x,y);


    // Set a boundary condition
    Func clampedStart;
    clampedStart(x, y, c) = floatingStart(clamp(x, 0, width-1), clamp(y, 0, height-1), c);
    Func clampedEnd;
    clampedEnd(x, y, c) = floatingEnd(clamp(x, 0, width-1), clamp(y, 0, height-1), c);


    //Func warped = warp(clampedStart, segmentsBefore, segmentsAfter, 5);
    Func warped = morph(clampedStart,clampedEnd, segmentsBefore, segmentsAfter, no_segments, time);

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

    output.compileToFile("morph", {time, no_segments, inputStart, inputEnd, segBefore, segAfter});


    return 0;
}

process.cpp

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <string>
#include <sstream>
extern "C" {
#include "morph.h"
}
#include "../Util.h"
#include "../png.h"
#include <sys/time.h>

//helper function that converts integer to string
std::string IntToString ( int number )
{
  std::ostringstream oss;

  // Works just like cout
  oss<< number;

  // Return the underlying string
  return oss.str();
}

int main(int argc, char **argv) {
    if (argc < 4) {
        printf("Usage: ./process input.png levels alpha beta output.png\n"
               "e.g.: ./process input.png 8 1 1 output.png\n");
        return 0;
    }
    //two images that we want to morph from one to the other
    Image<uint16_t> inputStart = load<uint16_t>(argv[1]);
    Image<uint16_t> inputEnd = load<uint16_t>(argv[2]);
    // number of morphing steps in between
    int N = atoi(argv[3]);


    // Image<float> segBefore = {{85, 131, 112, 128},{141, 130, 164, 131},{103, 199, 116, 216}, \
    //                           {163, 169, 146, 210}, {36, 176, 58, 229}};
    // Image<float> segAfter = {{81, 113, 109, 106}, {140, 104, 157, 102}, {101, 176, 125, 196}, \
    //                          {160, 131, 147, 176}, {64, 156, 85, 194}};

    // user-specified segments that depend on input images
    Image<float> segBefore = {{71, 93, 89, 93}, {113, 92, 132, 90}, {87, 132, 123, 128}, \
                              {64, 132, 104, 172}, {145, 129, 112, 171}, {70, 41, 121, 40}, {102, 94, 102, 119}};
    Image<float> segAfter = {{68, 100, 89, 98}, {115, 97, 134, 96}, {88, 142, 120, 141}, {56, 135, 105, 187},\
                             {146, 133, 115, 185}, {62, 49, 131, 47}, {100, 95, 102, 125}};


    float time;
    save(inputStart, "morph0.png");
    save(inputEnd, "morph"+ IntToString(N+1) + ".png");
    for (int i = 0; i < N; i++){
        Image<uint16_t> output(inputStart.width(), inputStart.height(), 3);
        time = (i+1)*(1.0f/float(N+1));
        morph(time,segBefore.height(), inputStart, inputEnd,segBefore,segAfter,output);
        save(output, "morph"+ IntToString(i+1) + ".png");
    }

    //some other segments picked by user for other images

 // Image<float> segBefore = {{85, 131, 112, 128},{141, 130, 164, 131},{103, 199, 116, 216}, \
    //                           {163, 169, 146, 210}, {36, 176, 58, 229}};
    // Image<float> segAfter = {{81, 113, 109, 106}, {140, 104, 157, 102}, {101, 176, 125, 196}, \
    //                          {160, 131, 147, 176}, {64, 156, 85, 194}};

// Image<float> segBefore = {{71.0f, 93.0f, 89.0f, 93.0f},
//                               {113.0f, 92.0f, 132.0f, 90.0f},
//                               {87, 132, 123, 128},
//                               {64, 132, 104, 172},
//                           {145, 129, 112, 171},
//                           {70, 41, 121, 40},
//                           {102, 94, 102, 119}};
//     Image<float> segAfter = {{68.0f, 100.0f, 89.0f, 98.0f},
//                              {115, 97, 134, 96},
//                              {88, 142, 120, 141},
//                              {56, 135, 105, 187},
//                               {146, 133, 115, 185},
//                               {2, 49, 131, 47},
//                               {100, 95, 102, 125}};


  // Image<float> segBefore = {{71.0f, 93.0f, 89.0f, 93.0f},
  //                             {113.0f, 92.0f, 132.0f, 90.0f},
  //                             {103.0f, 199.0f, 116.0f, 216.0f},
  //                             {163.0f, 169.0f, 146.0f, 210.0f},
  //                             {36.0f, 176.0f, 58.0f, 229.0f}};
  //   Image<float> segAfter = {{68.0f, 100.0f, 89.0f, 98.0f},
  //                            {115.0f, 97.0f, 134.0f, 96.0f},
  //                            {101.0f, 176.0f, 125.0f, 196.0f},
  //                            {160.0f, 131.0f, 147.0f, 176.0f},
  //                            {64.0f, 156.0f, 85.0f, 194.0f}};
    return 0;
}