Terrible Seam Carving

#define cimg_display 0
#include "lib/CImg.h"
#include "lib/slVector.H"

#include <math.h>   //For pow, abs, maybe a few other things
using namespace cimg_library;

//Get the Lab values of a pixel. For points not on the 2d array, "wrap around". The algorithm below should document itself really.
SlVector3 getPix ( int row, int col, int wid, int hei, int actWid, SlVector3 *image ) {

    //SWITCH TO RETURN 0.0 AT BOUNDARY IF RUSULTS UNFAVORABLE
    //Get wrapped pixel value if it's not withing image bounds, for energy computation.
    //Handle values above or to left of image
    //while ( row < 0 ) row += hei;
    //while ( col < 0 ) col += wid;
    //Handle values below or to right of image
    //row %= hei; col %= wid;

    //If beyond boundary, just return a 0 vector.
    if ( row < 0 || row >= hei || col < 0 || col >= wid )

        return SlVector3 ( 0, 0, 0 );
    return image[row*actWid + col];
}

//Get bounded array value, that is, get valid array value for point in array, get twice highest possible value for point outside array.
double getArBd ( int row, int col, int wid, int hei, double *en, double high ) {

    //Return something greater than any possible energy at boundaries, so that seams do not go across boundaries
    if ( row < 0 || row >= hei || col < 0 || col >= wid )

        return high * 2;
    return en[row*wid + col];
}

//Shift a row/column in array by 1 starting from some point. This is the double array variant.
void shift ( int rowStart, int colStart, int wid, int hei, int actWid, bool vert, double *array ) {

    for ( int ind = (vert?rowStart:colStart); ind < (vert?hei:wid) - 1; ind ++ )

        array[(vert?ind:rowStart)*actWid + (vert?colStart:ind)] = array[(vert?ind+1:rowStart)*actWid + (vert?colStart:ind+1)];
}

//Shift a row/column in array by 1 starting from some point. This is the SlVector3 array variant.
void shift ( int rowStart, int colStart, int wid, int hei, int actWid, bool vert, SlVector3 *array ) {

    for ( int ind = (vert?rowStart:colStart); ind < (vert?hei:wid) - 1; ind ++ )

        array[(vert?ind:rowStart)*actWid + (vert?colStart:ind)] = array[(vert?ind+1:rowStart)*actWid + (vert?colStart:ind+1)];
}

//Compute the gradient at each pixel using Scharr operator
double computeEnergy ( int row, int col, int wid, int hei, int actWid, SlVector3 *image ) {

    //If ranges are invalid, return early.
    if ( row < 0 || row >= hei || col < 0 || col >= wid ) return 0.0;

    //Calculate the horizontal and vertical gradients of each pixel, using common Scharr operator.
    SlVector3 difV = 3 * getPix( row-1, col-1, wid, hei, actWid, image )
                   +10 * getPix( row-1, col  , wid, hei, actWid, image )
                   + 3 * getPix( row-1, col+1, wid, hei, actWid, image )
                   - 3 * getPix( row+1, col-1, wid, hei, actWid, image )
                   -10 * getPix( row+1, col  , wid, hei, actWid, image )
                   - 3 * getPix( row+1, col+1, wid, hei, actWid, image );

    SlVector3 difH = 3 * getPix( row-1, col-1, wid, hei, actWid, image )
                   +10 * getPix( row  , col-1, wid, hei, actWid, image )
                   + 3 * getPix( row+1, col-1, wid, hei, actWid, image )
                   - 3 * getPix( row-1, col+1, wid, hei, actWid, image )
                   -10 * getPix( row  , col+1, wid, hei, actWid, image )
                   - 3 * getPix( row+1, col+1, wid, hei, actWid, image );

    //Return the sum of the magnitudes of the gradients
    return sqrt ( sqrMag ( difV ) + sqrMag ( difH ) );
}

//Get the minimum of three numbers.
inline double minOfThree ( double a, double b, double c ) {

    return ( a < b ? ( a < c ? a : c ) : ( b < c ? b : c ) );
}

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

    //Debug variables
    bool debugWriteEnergy = false;
    if ( argc == 6 )

        debugWriteEnergy = atoi( argv[5] );

    //Open the image proper and store the input and output dimensions and their differences
    CImg<double> input(argv[1]);
    CImg<double> lab = input.RGBtoLab();
    int inHei = input.height();
    int inWid = input.width();
    int outHei = atoi( argv[4] );
    int outWid = atoi( argv[3] );
    int difHei = inHei - outHei;
    int difWid = inWid - outWid;

    //Map all of lab's pixels to an SlVector3 array for easy pixel manipulation
    SlVector3 *image = new SlVector3 [inWid * inHei];
    for ( int row = 0; row < inHei; row++ )

        for ( int col = 0; col < inWid; col++ )

            for ( int iter = 0; iter < 3; iter ++ )

                image[row*inWid + col][iter] = lab( col, row, iter );

    //Get the energy and max energy
    double maxEnergy = 0.0;
    double *energy = new double [inWid*inHei];
    for ( int row = 0; row < inHei; row++ ) {

        for ( int col = 0; col < inWid; col++ ) {

            double en = computeEnergy ( row, col, inWid, inHei, inWid, image );
            energy[row*inWid + col] = en;
            maxEnergy = ( maxEnergy < en ? en : maxEnergy );
        }
    }

    if ( debugWriteEnergy ) {

        //Set all the image pixels based on their energy
        for ( int row = 0; row < inHei; row++ ) {

            for ( int col = 0; col < inWid; col++ ) {

                double set = pow ( energy[row*inWid + col] / maxEnergy, 0.25 );
                image[row*inWid + col][0] = (int) ( 100 * set );
                image[row*inWid + col][1] = 0;
                image[row*inWid + col][2] = 0;
            }
        }
    } else {

        //Select a certain number of horizontal seams and remove them.
        for ( int hSeams = 0; hSeams < difHei; hSeams++ ) {

            //Create and calculate the seam values
            //Set the first column of seam values equal to the energies at those points.
            double *seamMap = new double [inWid*(inHei-hSeams)];
            for ( int row = 0; row < inHei-hSeams; row++ )

                seamMap[row*inWid] = energy[row*inWid];

            //For every subsequent column, set each pixel's seam value equal to the sum of its energy and the least of all the left neighbor's energies.
            for ( int col = 1; col < inWid ; col++ )

                for ( int row = 0; row < inHei-hSeams; row++ )

                    seamMap[row*inWid + col] = energy[row*inWid + col]
                                 + minOfThree ( getArBd(row-1,col-1,inWid,inHei-hSeams,seamMap,maxEnergy*inWid),
                                        getArBd(row  ,col-1,inWid,inHei-hSeams,seamMap,maxEnergy*inWid),
                                        getArBd(row+1,col-1,inWid,inHei-hSeams,seamMap,maxEnergy*inWid));

            //Store each row index for every column, starting from last and working backward. This way we can easily get each pixel of the seam.
            int *rowIndicies = new int [inWid];

            //Find lowest value seam ending by looking at ending columns and finding lowest cumulative seam.
            rowIndicies[inWid-1] = 0;
            double minSeam = getArBd ( 0, inWid-1, inWid, inHei-hSeams, seamMap, maxEnergy*inWid );
            for ( int row = 1; row < inHei - hSeams; row++ ) {

                double current = getArBd ( row, inWid-1, inWid, inHei-hSeams, seamMap, maxEnergy*inWid );
                if ( current < minSeam ) { minSeam = current; rowIndicies[inWid-1] = row; }
            }

            //Work from ending to recover entire least important seam.
            for ( int col = inWid - 2; col >= 0; col -- ) {

                double a = getArBd ( rowIndicies[col+1]-1, col, inWid, inHei-hSeams, seamMap, maxEnergy*inWid );
                double b = getArBd ( rowIndicies[col+1]  , col, inWid, inHei-hSeams, seamMap, maxEnergy*inWid );
                double c = getArBd ( rowIndicies[col+1]+1, col, inWid, inHei-hSeams, seamMap, maxEnergy*inWid );
                double min = minOfThree ( a, b, c );
                if ( min == a )         rowIndicies[col] = rowIndicies[col+1]-1;
                else if ( min == b )    rowIndicies[col] = rowIndicies[col+1]  ;
                else                    rowIndicies[col] = rowIndicies[col+1]+1;
            }

            //Call shift function on the seam point of each column of image and energy; get rid of that seam in the image and in energy map
            for ( int col = 0; col < inWid; col ++ ) {

                shift ( rowIndicies[col], col, inWid, inHei - hSeams, inWid, true, image );
                shift ( rowIndicies[col], col, inWid, inHei - hSeams, inWid, true, energy );
            }

            //Recompute the energy at the locations where the seam was removed, and at the locations exactly below that
            for ( int col = 0; col < inWid; col ++ ) {

                for ( int dif = -2; dif <= 2; dif ++ ) {

                    int index = rowIndicies[col] + dif;
                    if ( index >= 0 && index < (inHei-hSeams) )

                        energy[(index)*inWid + col] = computeEnergy ( index, col, inWid, inHei-hSeams-1, inWid, image );
                }
            }
            delete [] rowIndicies;
            delete [] seamMap;
        }

        //Select a certain number of horizontal seams and remove them...
        for ( int vSeams = 0; vSeams < difWid; vSeams++ ) {

            //Create and calculate the seam values
            //Set the first row of seam values equal to the energies at those points.
            double *seamMap = new double [(inWid-vSeams)*outHei];
            for ( int col = 0; col < inWid-vSeams; col++ )

                seamMap[col] = energy[col];

            //For every subsequent row, set each pixel's seam value equal to the sum of its energy and the least of all the top neighbor's energies.
            for ( int row = 1; row < outHei ; row++ )

                for ( int col = 0; col < inWid-vSeams; col++ )

                    seamMap[row*(inWid-vSeams) + col] = energy[row*inWid + col]
                                                      + minOfThree ( getArBd(row-1,col-1,inWid-vSeams,outHei,seamMap,maxEnergy*outHei),
                                                                     getArBd(row-1,col  ,inWid-vSeams,outHei,seamMap,maxEnergy*outHei),
                                                                     getArBd(row-1,col+1,inWid-vSeams,outHei,seamMap,maxEnergy*outHei));

            //Store each column index for every row, starting from last and working backward. This way we can easily get each pixel of the seam.
            int *colIndicies = new int [outHei];

            //Find lowest value seam ending by looking at ending columns and finding lowest cumulative seam.
            colIndicies[outHei-1] = 0;
            double minSeam = getArBd ( outHei-1, 0, inWid-vSeams, outHei, seamMap, maxEnergy*outHei );
            for ( int col = 1; col < inWid-vSeams; col++ ) {

                double current = getArBd ( outHei-1, col, inWid-vSeams, outHei, seamMap, maxEnergy*inWid );
                if ( current < minSeam ) { minSeam = current; colIndicies[outHei-1] = col; }
            }

            //Work from ending to recover entire least important seam.
            for ( int row = outHei - 2; row >= 0; row -- ) {

                double a = getArBd ( row, colIndicies[row+1]-1, inWid-vSeams, outHei, seamMap, maxEnergy*outHei );
                double b = getArBd ( row, colIndicies[row+1]  , inWid-vSeams, outHei, seamMap, maxEnergy*outHei );
                double c = getArBd ( row, colIndicies[row+1]+1, inWid-vSeams, outHei, seamMap, maxEnergy*outHei );
                double min = minOfThree ( a, b, c );
                if ( min == a )         colIndicies[row] = colIndicies[row+1]-1;
                else if ( min == b )    colIndicies[row] = colIndicies[row+1]  ;
                else                    colIndicies[row] = colIndicies[row+1]+1;
            }

            //Call shift function on the seam point of each column of image and energy; get rid of that seam in the image and in energy map
            for ( int row = 0; row < outHei; row ++ ) {

                shift ( row, colIndicies[row], inWid-vSeams, outHei, inWid, false, image );
                shift ( row, colIndicies[row], inWid-vSeams, outHei, inWid, false, energy );
            }

            //Recompute the energy at the locations where the seam was removed, and at the locations exactly below that
            for ( int row = 0; row < outHei; row ++ ) {

                for ( int dif = -2; dif <= 2; dif ++ ) {

                    int index = colIndicies[row] + dif;
                    if ( index >= 0 && index < (inWid-vSeams) )

                        energy[row*inWid + index] = computeEnergy ( row, index, inWid-vSeams-1, outHei, inWid, image );
                }
            }
            delete [] colIndicies;
            delete [] seamMap;
        }
    }

    delete [] energy;

    //Open the output image and place all map all the processing array points to it.
    CImg<double> output( outWid, outHei, input.depth(), input.spectrum(), 0 );
    for ( int row = 0; row < outHei; row++ )

        for ( int col = 0; col < outWid; col++ )

            for ( int iter = 0; iter < 3; iter ++ )

                output( col, row, iter ) = image[row*inWid + col][iter];

    //Save the image as the appropriate format, deallocate processing array and exit.
    CImg<double> rgb = output.LabtoRGB();
    if ( strstr( argv[2], "png" ) )         rgb.save_png( argv[2] );
    else if ( strstr( argv[2], "jpg" ) )    rgb.save_jpeg( argv[2] );

    delete [] image;
    return 0;
}