Spline Catmull-Rom for image zooming using C++ and opencv

Question

I'm trying to implement spline Catmull-Rom for image zooming using C++ and OpenCV. I performed two tests, the first is image zooming (X2), and the second image reconstruction (zooming image decimated). My problem is that in the image interpolated appear some white and black pixel (image1) when I displayed the value of pixels I found that white pixels have a negative value and the black one has a value greater than 255, also the image reconstructed appear blurred (image2 and image3).

 float CalCurveInt(float t, float p0, float p1, float p2, float p3)
    {
    
        float t2 = t * t;
        float t3 = t2 * t;
    
        float x = 0.5f * ((2.0f * p1) +
            (-p0 + p2) * t +
            (2.0f * p0 - 5.0f * p1 + 4 * p2 - p3) * t2 +
            (-p0 + 3.0f * p1 - 3.0f * p2 + p3) * t3);
    return x;
    }

    Mat CalcCatmull(Mat &src, int zoom)
    {
        int v1, v2, v3, v4, Ptr, Xmax, Ymax;
        float Result, t, c1, c2, c3, c4;
        //------------------------------------------------------------
        Xmax = src.cols;
        Ymax = src.rows;
        Size srcSize(zoom*Xmax, Ymax);
        Mat dst(srcSize, CV_8UC1);
    
        for (int j = 0; j < Ymax; j++)
        {
            Ptr = 0;
            for (int i = 0; i < Xmax; i++)
            {
                v1 = i - 1; v2 = i; v3 = i + 1; v4 = i + 2;
    
                if (i - 1 < 0)      v1 = 0;
                if (Xmax <= i + 1)  v3 = Xmax - 1;
                if (Xmax <= i + 2)  v4 = Xmax - 1;
    
                for (double J = 1; J <= zoom; J++)
                {
                    t = J / zoom;
                    Result = 0.0;
    
                    c1 = src.at<uchar>(j, v1);
                    c2 = src.at<uchar>(j, v2);
                    c3 = src.at<uchar>(j, v3);
                    c4 = src.at<uchar>(j, v4);
    
                    Result = CalCurveInt(t, c1, c2, c3, c4);
                    dst.at<uchar>(j, Ptr) = abs(Result);
                    Ptr++;
                }
            }
        }

        //------------------------------------------------
        Xmax = dst.cols;
        Ymax = dst.rows;
    
        Size srcSize1(Xmax, zoom*Ymax);
        Mat dest(srcSize1, CV_8UC1);
    
        for (int i = 0; i < Xmax; i++)
        {
            Ptr = 0;
            for (int j = 0; j < Ymax; j++)
            {
    
                v1 = j - 1; v2 = j; v3 = j + 1; v4 = j + 2;
    
                if (j - 1 < 0)      v1 = 0;
                if (Ymax <= j + 1)  v3 = Ymax - 1;
                if (Ymax <= j + 2)  v4 = Ymax - 1;
    
                for (double J = 1; J <= zoom; J++)
                {
                    t = J / zoom;
                    Result = 0.0;
    
                    c1 = dst.at<uchar>(v1, i);
                    c2 = dst.at<uchar>(v2, i);
                    c3 = dst.at<uchar>(v3, i);
                    c4 = dst.at<uchar>(v4, i);
                    Result = CalCurveInt(t, c1, c2, c3, c4);           
                    dest.at<uchar>(Ptr, i) = Result;
                    Ptr++;
                }
            }
        }
    
        return dest;
    }
    
    float zoom = 2.0;

    
    int main()
    {
        Mat src = imread("fruits.png", CV_LOAD_IMAGE_GRAYSCALE);
        int width = src.cols;
        int hight = src.rows;
    
        /*Image Decimation*/
        Size srcdSize(int(width / zoom), int(hight / zoom));
        Mat srcd;
        pyrDown(src, srcd, srcdSize);
        imshow("decimation", srcd);
    
        Mat dst = CalcCatmull(srcd, zoom);
    
        imshow("Image Source", src);
        imshow("Image dest", dst);
        imwrite("Image dest.png", dst);
    
        waitKey(0);
        return 0;
    }

Thanks in advance.

Answer 1

My old implementation, seems it worked fine.

#include <iostream>
#include <vector>
#include <stdio.h>
#include <stdarg.h>
#include "opencv2/opencv.hpp"
#include "fstream"
#include "iostream"
using namespace std;
using namespace cv;

//-----------------------------------------------------------------------------------------------------
// Take 2 points, compute values between p1 and p2, p0 and p3 need for tangents computation
// on the bouunds. Parameter t - changes in range 0 to 1 (0 - we are in p1, 1 - we are in p2)
//-----------------------------------------------------------------------------------------------------
void PointOnCurve(Point2f &out, float t, Point2f p0, Point2f p1, Point2f p2, Point2f p3)
{
    float t2 = t * t;
    float t3 = t2 * t;
    out.x = 0.5f * ( ( 2.0f * p1.x ) + ( -p0.x + p2.x ) * t +
        ( 2.0f * p0.x - 5.0f * p1.x + 4 * p2.x - p3.x ) * t2 +
        ( -p0.x + 3.0f * p1.x - 3.0f * p2.x + p3.x ) * t3 );
    out.y = 0.5f * ( ( 2.0f * p1.y ) +  ( -p0.y + p2.y ) * t +
        ( 2.0f * p0.y - 5.0f * p1.y + 4 * p2.y - p3.y ) * t2 +
        ( -p0.y + 3.0f * p1.y - 3.0f * p2.y + p3.y ) * t3 );
}
//-----------------------------------------------------------------------------------------------------
// interpolation of 4х4 patch
// 
//          S * S * S * S 
//          * * * * * * * 
//          S * S * S * S 
//          * * * * * * * 
//          S * S * S * S 
//          * * * * * * * 
//          S * S * S * S 
//  
//          S- pixels of source imgage
//
//          sequentially take 2 middle columns and computte D.
//
//          S * 1 * 2 * S
//          * * * * * * *
//          S * 1 * 2 * S
//          * * D * D * *
//          S * 1 * 2 * S
//          * * * * * * *
//          S * 1 * 2 * S
//
//          same for rows and we will have F
//
//          S * S * S * S
//          * * * * * * *
//          3 * 3 F 3 * 3
//          * * D * D * *
//          4 * 4 F 4 * 4
//          * * * * * * *
//          S * S * S * S
//
//          then compute diagonals and after averafing with neihbours will find С
//
//          1 * S * S * 2
//          * * * * * * *
//          S * 1 F 2 * S
//          * * D C D * *
//          S * 2 F 1 * S
//          * * * * * * *
//          2 * S * S * 1
//-----------------------------------------------------------------------------------------------------
void PointOnSurface(Mat& src,Mat& dst)
{
float t=0.5;
Point2f out;
dst=Mat(3,3,CV_32FC1);
// Угловые точки результата совпадают с точками центральной ячейки исходного патча
dst.at<float>(0,0)=src.at<float>(1,1);
dst.at<float>(2,0)=src.at<float>(2,1);

dst.at<float>(0,2)=src.at<float>(1,2);
dst.at<float>(2,2)=src.at<float>(2,2);

Point2f p0;
Point2f p1;
Point2f p2;
Point2f p3;

p0.x=0;p0.y=src.at<float>(0,1);
p1.x=1;p1.y=src.at<float>(1,1);
p2.x=2;p2.y=src.at<float>(2,1);
p3.x=3;p3.y=src.at<float>(3,1);

PointOnCurve(out,t,p0,p1,p2,p3);
dst.at<float>(1,0)=out.y;

p0.x=0;p0.y=src.at<float>(0,2);
p1.x=1;p1.y=src.at<float>(1,2);
p2.x=2;p2.y=src.at<float>(2,2);
p3.x=3;p3.y=src.at<float>(3,2);

PointOnCurve(out,t,p0,p1,p2,p3);
dst.at<float>(1,2)=out.y;

p0.x=0;p0.y=src.at<float>(1,0);
p1.x=1;p1.y=src.at<float>(1,1);
p2.x=2;p2.y=src.at<float>(1,2);
p3.x=3;p3.y=src.at<float>(1,3);


PointOnCurve(out,t,p0,p1,p2,p3);
dst.at<float>(0,1)=out.y;

p0.x=0;p0.y=src.at<float>(2,0);
p1.x=1;p1.y=src.at<float>(2,1);
p2.x=2;p2.y=src.at<float>(2,2);
p3.x=3;p3.y=src.at<float>(2,3);

PointOnCurve(out,t,p0,p1,p2,p3);
dst.at<float>(2,1)=out.y;

// diagonals

// 1
p0.x=0;p0.y=src.at<float>(0,0);
p1.x=1;p1.y=src.at<float>(1,1);
p2.x=2;p2.y=src.at<float>(2,2);
p3.x=3;p3.y=src.at<float>(3,3);


PointOnCurve(out,t,p0,p1,p2,p3);
float d1=out.y;

// 2
p0.x=0;p0.y=src.at<float>(3,0);
p1.x=1;p1.y=src.at<float>(2,1);
p2.x=2;p2.y=src.at<float>(1,2);
p3.x=3;p3.y=src.at<float>(0,3);

PointOnCurve(out,t,p0,p1,p2,p3);
float d2=out.y;

// averaging
dst.at<float>(1,1)=1.0/6.0*(d1+d2+dst.at<float>(0,1)+dst.at<float>(1,0)+dst.at<float>(1,2)+dst.at<float>(2,1));
}
//-----------------------------------------------------------------------------------------------------
//
//-----------------------------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------------------------
void Scale2Times(Mat& src_img,Mat &dstImg)
{
    Mat imgf,img;
    Mat dst;
    Mat src;
    src_img.convertTo(imgf,CV_32FC1,1.0/255.0);

    cv::copyMakeBorder(imgf,img,1,1,1,1,cv::BORDER_REFLECT);

    dstImg=Mat(src_img.rows*2,src_img.cols*2,CV_32FC1);

    for(int i=0;i<img.rows-4;i++)
    {
        for(int j=0;j<img.cols-4;j++)
        {
            img(Rect(j,i,4,4)).copyTo(src);
            PointOnSurface(src,dst);
            dst.copyTo(dstImg(Rect(2*j+1,2*i+1,3,3)));
        }
    }
    dstImg=dstImg(Rect(0,0,dstImg.cols-2,dstImg.rows-2)).clone();
}
//-----------------------------------------------------------------------------------------------------
//
//-----------------------------------------------------------------------------------------------------
int main( int argc, char** argv )
{
    namedWindow("Src");
    namedWindow("cvResize");
    namedWindow("Catmul-Rom");
    Mat Img=imread("C:\\ImagesForTest\\1.tiff",0);
    imshow("Src",Img);
    Mat dstImg;
    Scale2Times(Img,dstImg);
    imshow("Catmul-Rom",dstImg);

    Mat ImgLin(Img.rows*2,Img.cols*2,CV_8UC1);
    cv::resize(Img,ImgLin,Size(Img.cols*2,Img.rows*2),INTER_CUBIC);
    imshow("cvResize",ImgLin);

    waitKey(0);
    //getchar();

    return 0;
}