## スパース行列(SparseMat)の使い方2：超解像

22 7月 2010 Under: opencv2.x-samples

SparseMatをつかった超解像処理を行います．このデモンストレーションは下記論文の実装になっています．詳細は以下論文やビデオをご覧ください．また，このサンプルコードはOpenMPを有効化するとコアを有効に使います．メモリは１Gほど使うため実行時には注意してください．

Farsiu, S.,Robinson, D., Elad, M., Milanfar, P.”Fast and robust multiframe super resolution,” IEEETrans.ImageProcessing 13 (2004)1327?1344.

# C++

//Super resolution with Bilateral Total Variation
//Implimentation of a paper;
//Farsiu, S.,Robinson, D., Elad, M., Milanfar, P."Fast and robust multiframe super resolution," IEEETrans.ImageProcessing 13 (2004)1327?1344.

#include <cv.h>
#include <highgui.h>
#ifdef _OPENMP
#include <omp.h>
#endif
using namespace cv;

#include <iostream>
using namespace std;

//sparse matrix and matrix function
void mulSparseMat32f(SparseMat& smat, Mat& src, Mat& dest, bool isTranspose = false);
Mat visualizeSparseMat(SparseMat& smat, Size out_imsize);
void subtract_sign(Mat& src1, Mat&src2, Mat& dest);//dest = sign(src1-src2);

//image processing
double getPSNR(Mat& src1, Mat& src2, int bb);
void addspikenoise(Mat& src, Mat& dest, int val);
void addgaussnoise(Mat& src, Mat& dest, double sigma);

//for super resolution
SparseMat createDownsampledMotionandBlurCCDSparseMat32f(Mat& src, int amp, Point2d move);
SparseMat createDegradedImageandSparseMat32F(Mat& src, Mat& dest, Point2d move, int amp);
void btvregularization(Mat& srcVec, Size kernel, float alpha, Mat& dstVec, Size size);

#define sign_float(a,b) (a>b)?1.0f:(a<b)?-1.0f:0.0f

enum
{
SR_DATA_L1=0,
SR_DATA_L2
};

//Bilateral Total Variation based Super resolution with
void superresolutionSparseMat32f(Mat src[], Mat& dest, SparseMat DHF[], const int numofview,int iteration, float beta, float lambda, float alpha, Size reg_window,int method, Mat& ideal);

int main(int argc, char** argv)
{
Mat image = imread("lenna.png");//input image
if(image.empty())
{
cout<<"invalid image name"<<endl;
return -1;
}
Mat dest = Mat(image.size(),CV_8UC3);

const int image_count = 16;//number of input images for super resolution
const int rfactor = 4;//magnification factor

Point2d move[image_count];//memory of motion of input images
SparseMat A[image_count];//degrading matrices
Mat degimage[image_count];//degraded images

RNG rnd;//random number generator of image shift/movement

//(1) generate degraded images and degrading matrices for super resolution
for(int i=0;i<image_count;i++)
{
cout<<i<<endl;
move[i].x=rnd.uniform(0.0,(double)rfactor);//randam shift for x
move[i].y=rnd.uniform(0.0,(double)rfactor);//randam shift for y
if(i==0)// fix first image
{
move[i].x=0;
move[i].y=0;
}
Mat imtemp(image.rows/rfactor,image.cols/rfactor,CV_8UC3);
degimage[i].create(image.rows/rfactor,image.cols/rfactor,CV_8UC3);

A[i]=createDegradedImageandSparseMat32F(image, imtemp,move[i],rfactor);
addgaussnoise(imtemp,degimage[i],10.0);//add gaussian noise
addspikenoise(degimage[i],degimage[i],500);//add spike noise

imshow("im",degimage[i]);
char name[32];
sprintf(name,"input%03d.png",i);
imwrite(name,degimage[i]);
waitKey(30);
}

//(2) super resolution
//beta: asymptotic value of steepest descent method
//lambda: weight parameter to balance data term and smoothness term
////alpha: parameter of spacial distribution in btv
//btv kernel size: kernel size of btv filter
superresolutionSparseMat32f(degimage,dest,A,image_count,
180,//number of iteration
1.3f,//beta
0.03f,//lambda
0.7f,//alpha
Size(7,7),//btv kernel size
SR_DATA_L1,//L1 norm minimization for data term
image);//ideal image for evaluation
return 0;
}

void sum_float_OMP(Mat src[], Mat& dest, int numofview, float beta)
{
for(int n=0;n<numofview;n++)
{
#pragma omp parallel for
for(int j=0;j<dest.rows;j++)
{
dest.ptr<float>(j)[0]-=beta*src[n].ptr<float>(j)[0];
dest.ptr<float>(j)[1]-=beta*src[n].ptr<float>(j)[1];
dest.ptr<float>(j)[2]-=beta*src[n].ptr<float>(j)[2];
}
}
}
void superresolutionSparseMat32f(Mat src[], Mat& dest, SparseMat DHF[], const int numofview,int iteration, float beta, float lambda, float alpha, Size reg_window,int method, Mat& ideal)
{
//(3) create initial image by simple linear interpolation
resize(src[0],dest,dest.size());
cout<<"PSNR"<<getPSNR(dest,ideal,10)<<"dB"<<endl;
imwrite("linear.png",dest);

//(4)convert Mat image structure to 1D vecor structure
Mat dstvec;
dest.reshape(3,dest.cols*dest.rows).convertTo(dstvec,CV_32FC3);

Mat* dstvectemp=new Mat[numofview];
Mat* svec = new Mat[numofview];
Mat* svec2 = new Mat[numofview];

for(int n=0;n<numofview;n++)
{
src[n].reshape(3,src[0].cols*src[0].rows).convertTo(svec[n],CV_32FC3);
src[n].reshape(3,src[0].cols*src[0].rows).convertTo(svec2[n],CV_32FC3);

dstvectemp[n]=dstvec.clone();
}

Mat reg_vec=Mat::zeros(dest.rows*dest.cols,1,CV_32FC3);//regularization vector

//(5)steepest descent method for L1 norm minimization
for(int i=0;i<iteration;i++)
{
cout<<"iteration"<<i<<endl;
int64 t = getTickCount();
Mat diff=Mat::zeros(dstvec.size(),CV_32FC3);

//(5-1)btv
if(lambda>0.0) btvregularization(dstvec,reg_window,alpha,reg_vec,dest.size());

#pragma omp parallel for
for(int n=0;n<numofview;n++)
{
//degrade current estimated image
mulSparseMat32f(DHF[n],dstvec,svec2[n]);

//compere input and degraded image
Mat temp(src[0].cols*src[0].rows,1, CV_32FC3);
if(method==SR_DATA_L1)
{
subtract_sign(svec2[n], svec[n],temp);
}
else
{
subtract(svec2[n],svec[n],temp);
//temp = svec2[n]- svec[n]; //supported in OpenCV2.1
}

//blur the subtructed vector with transposed matrix
mulSparseMat32f(DHF[n],temp,dstvectemp[n],true);
}
//creep ideal image, beta is parameter of the creeping speed.
//add transeposed difference vector. sum_float_OMP is parallelized function of following for loop
/*for(int n=0;n<numofview;n++)
{
addWeighted(dstvec,1.0,dstvectemp[n],-beta,0.0,dstvec);
//dstvec -= (beta*dstvectemp[n]);//supported in OpenCV2.1
}*/
sum_float_OMP(dstvectemp,dstvec,numofview,beta);

//add smoothness term
if(lambda>0.0)
{
addWeighted(dstvec,1.0,reg_vec,-beta*lambda,0.0,dstvec);
//dstvec -=lambda*beta*reg_vec;//supported in OpenCV2.1
}

//show SR imtermediate process information. these processes does not be required at actural implimentation.
dstvec.reshape(3,dest.rows).convertTo(dest,CV_8UC3);
cout<<"PSNR"<<getPSNR(dest,ideal,10)<<"dB"<<endl;

char name[64];
sprintf(name,"%03d: %.1f dB",i,getPSNR(dest,ideal,10));
putText(dest,name,Point(15,50), FONT_HERSHEY_DUPLEX,1.5,CV_RGB(255,255,255),2);

sprintf(name,"iteration%04d.png",i);
imshow("SRimage",dest);waitKey(30);
imwrite(name,dest);
cout<<"time/iteration"<<(getTickCount()-t)*1000.0/getTickFrequency()<<"ms"<<endl;
}

//re-convert  1D vecor structure to Mat image structure
dstvec.reshape(3,dest.rows).convertTo(dest,CV_8UC3);
imwrite("sr.png",dest);
}

void btvregularization(Mat& srcVec, Size kernel, float alpha, Mat& dstVec, Size size)
{
Mat src;
srcVec.reshape(3,size.height).convertTo(src,CV_32FC3);
Mat dest(size.height,size.width,CV_32FC3);

const int kw = (kernel.width-1)/2;
const int kh = (kernel.height-1)/2;

float* weight = new float[kernel.width*kernel.height];
for(int m=0,count=0;m<=kh;m++)
{
for(int l=kw;l+m>=0;l--)
{
weight[count]=pow(alpha,abs(m)+abs(l));
count++;
}
}
//a part of under term of Eq (22) lambda*\sum\sum ...
#pragma omp parallel for
for(int j=kh;j<src.rows-kh;j++)
{
float* s = src.ptr<float>(j);
float* d = dest.ptr<float>(j);
for(int i=kw;i<src.cols-kw;i++)
{
float r=0.0;
float g=0.0;
float b=0.0;

const float sr=s[3*(i)+0];
const float sg=s[3*(i)+1];
const float sb=s[3*(i)+2];
for(int m=0,count=0;m<=kh;m++)
{
float* s2 = src.ptr<float>(j-m);
float* ss = src.ptr<float>(j+m);
for(int l=kw;l+m>=0;l--)
{
r+=weight[count]*(sign_float(sr,ss[3*(i+l)+0]) -sign_float(s2[3*(i-l)+0],sr));
g+=weight[count]*(sign_float(sg,ss[3*(i+l)+1]) -sign_float(s2[3*(i-l)+1],sg));
b+=weight[count]*(sign_float(sb,ss[3*(i+l)+2]) -sign_float(s2[3*(i-l)+2],sb));
count++;
}
}
d[3*i+0]=r;
d[3*i+1]=g;
d[3*i+2]=b;
}
}
dest.reshape(3,size.height*size.width).convertTo(dstVec,CV_32FC3);
delete[] weight;
}

SparseMat createDownsampledMotionandBlurCCDSparseMat32f(Mat& src, int amp, Point2d move)
{
//(1)'
//D down sampling matrix
//H blur matrix, in this case, we use only ccd sampling blur.
//F motion matrix, in this case, threr is only global shift motion.

float div = 1.0f/((float)(amp*amp));
int x1 = (int)(move.x+1);
int x0 = (int)(move.x);
float a1 = (float)(move.x-x0);
float a0 = (float)(1.0-a1);

int y1 = (int)(move.y+1);
int y0 = (int)(move.y);
float b1 = (float)(move.y-y0);
float b0 = (float)(1.0-b1);

int bsx =x1;
int bsy =y1;

int matsize =src.cols*src.rows ;
int matsize2 =src.cols*src.rows/(amp*amp);
int size2[2]={matsize,matsize2};
SparseMat DHF(2,size2,CV_32FC1);

const int step = src.cols/amp;
for(int j=0;j<src.rows;j+=amp)
{
for(int i=0;i<src.cols;i+=amp)
{
int y = src.cols*j+i;
int s = step*j/amp+i/amp;

if(i>=bsx &&i<src.cols-bsx-amp &&j>=bsy &&j<src.rows-bsy-amp)
{
for(int l=0;l<amp;l++)
{
for(int k=0;k<amp;k++)
{
DHF.ref<float>(y+src.cols*(y0+l)+x0+k,s)+=a0*b0*div;
DHF.ref<float>(y+src.cols*(y0+l)+x1+k,s)+=a1*b0*div;
DHF.ref<float>(y+src.cols*(y1+l)+x0+k,s)+=a0*b1*div;
DHF.ref<float>(y+src.cols*(y1+l)+x1+k,s)+=a1*b1*div;
}
}
}
}
}
return DHF;
}
SparseMat createDegradedImageandSparseMat32F(Mat& src, Mat& dest, Point2d move, int amp)
{
SparseMat DHF=createDownsampledMotionandBlurCCDSparseMat32f(src,amp,move);

int matsize =src.cols*src.rows ;
int matsize2 =src.cols*src.rows/(amp*amp);

Mat svec;
src.reshape(3,matsize).convertTo(svec,CV_32FC3);
Mat dvec(matsize2,1,CV_32FC3);

mulSparseMat32f(DHF,svec,dvec);

imshow("smat",visualizeSparseMat(DHF,Size(512,512/amp/amp)));waitKey(30);
//re-convert  1D vecor structure to Mat image structure
dvec.reshape(3,dest.rows).convertTo(dest,CV_8UC3);

return DHF;
}

void subtract_sign(Mat& src1, Mat&src2, Mat& dest)
{
for(int j=0;j<src1.rows;j++)
{
float* s1 = src1.ptr<float>(j);
float* s2 = src2.ptr<float>(j);
float* d = dest.ptr<float>(j);
for(int i=0;i<src1.cols;i++)
{
d[3*i]=sign_float(s1[3*i],s2[3*i]);
d[3*i+1]=sign_float(s1[3*i+1],s2[3*i+1]);
d[3*i+2]=sign_float(s1[3*i+2],s2[3*i+2]);
}
}
}

Mat visualizeSparseMat(SparseMat& smat, Size out_imsize)
{
Mat data = Mat::zeros(out_imsize,CV_8U);
double inv_size0 = 1.0/smat.size(0);
double inv_size1 = 1.0/smat.size(1);

SparseMatIterator it = smat.begin(),it_end = smat.end();
for(;it!=it_end;++it)
{
int j = (int)(((double)it.node()->idx[0]*inv_size0*out_imsize.height));
int i = (int)(((double)it.node()->idx[1]*inv_size1*out_imsize.width));
data.at<uchar>(j,i)=255;
}

Mat zeromat = Mat::zeros(out_imsize,CV_8U);
vector<Mat> image;
image.push_back(zeromat);
image.push_back(data);
image.push_back(zeromat);

Mat ret;
cv::merge(image,ret);

cout<<"number of non zero elements: "<<smat.nzcount()<<endl;
return ret;
}

void addgaussnoise(Mat& src, Mat& dest, double sigma)
{
Mat noise(src.rows,src.cols,CV_32FC1);
Mat src_f;
vector<Mat> images;
split(src,images);
for(int i=0;i<src.channels();i++)
{
images[i].convertTo(src_f,CV_32FC1);
randn(noise,Scalar(0.0),Scalar(sigma));
Mat temp = noise+src_f;
temp.convertTo(images[i],CV_8UC1);
}
merge(images,dest);
}

void addspikenoise(Mat& src, Mat& dest, int val)
{
src.copyTo(dest);
#pragma omp parallel for
for(int j=0;j<src.rows;j++)
{
RNG rnd(getTickCount());
uchar* d = dest.ptr<uchar>(j);
for(int i=0;i<src.cols;i++)
{
if(rnd.uniform(0,val)<1)
{
d[3*i]=255;
d[3*i+1]=255;
d[3*i+2]=255;
}
}
}
}

void mulSparseMat32f(SparseMat& smat, Mat& src, Mat& dest, bool isTranspose)
{
dest.setTo(0);
SparseMatIterator it = smat.begin(),it_end = smat.end();
if(!isTranspose)
{
for(;it!=it_end;++it)
{
int i=it.node()->idx[0];
int j=it.node()->idx[1];
float* d = dest.ptr<float>(j);
float* s = src.ptr<float>(i);
for(int i=0;i<3;i++)
{
d[i]+= it.value<float>() * s[i];
}
}
}
else // for transpose matrix
{
for(;it!=it_end;++it)
{
int i=it.node()->idx[1];
int j=it.node()->idx[0];
float* d = dest.ptr<float>(j);
float* s = src.ptr<float>(i);
for(int i=0;i<3;i++)
{
d[i]+= it.value<float>() * s[i];
}
}
}
}

double getPSNR(Mat& src1, Mat& src2, int bb)
{
int i,j;
double sse,mse,psnr;
sse = 0.0;

Mat s1,s2;
cvtColor(src1,s1,CV_BGR2GRAY);
cvtColor(src2,s2,CV_BGR2GRAY);

int count=0;
for(j=bb;j<s1.rows-bb;j++)
{
uchar* d=s1.ptr(j);
uchar* s=s2.ptr(j);

for(i=bb;i<s1.cols-bb;i++)
{
sse += ((d[i] - s[i])*(d[i] - s[i]));
count++;
}
}
if(sse == 0.0 || count==0)
{
return 0;
}
else
{
mse =sse /(double)(count);
psnr = 10.0*log10((255*255)/mse);
return psnr;
}
}



## （１）劣化画像と劣化行列（疎行列）を生成します．

Y=DHFX

で表すことができます．行列DHFををまとめて書けば，

Y=AX

となります．最後にガウシアンノイズとスパイクノイズを劣化画像に追加します．これを図に示すと以下のようになります．

## (2)　劣化の逆変換を最急降下法で行うことで超解像処理を行います．

ぼけ除去のサンプルに加えて，Bilateral Total Variation(BTV)による正則化が追加されており，ノイズにロバストになっています．ここでは，この式の計算を繰り返し処理により行います．（L1ノルム最小化の場合

Xn+1Xn – β{Σk=1NAkTsign(AkXn-Y)+λΣl=-ppΣm=pp[I-Sx-lSy-m]sign(Xn-SxlSymXn)}

βは収束のステップを，λは画像の滑らかさの拘束の強さを，αはBTVの距離関する減衰のパラメータとなっています．また，SR_DATA_L2，SR_DATA_L１で最小化するノルムを選択可能になっています．数式中，SxlSymは画像をｘ,y方向にl,mピクセル平行移動する行列を表しており，Iは何もしない単位行列となっています．また関数signは符号関数で，正の場合1を負の場合-1を返します．

## （３）まずはじめに初期値として，リサイズ（線形補間）した画像を入力します．

ここから超解像処理の詳細に入ります．ここでは，初期画像のPSNR（画質の評価関数）も測定しています．

## （５）最急降下法により，少しずつ近づけていきます．

1.まず，BTVによるペナルティを計算します．このペナルティは画像が急激に変化するところで大きくなります．

2.複数の入力画像と劣化行列をかけた現推定画像と比較（L1かL2ノルムで）し，残差を記憶します．

3.それらを足した後，正解画像との比較を行ってどれくらいの精度が出ているか確認します．実際の使用時には，正解画像はないため，ここではただ評価のためだけに正解画像を使用しています．

## 結果

16枚の解像度128×128の入力画像から，線形補間で拡大した画像と超解像処理を行ったものを以下に示します．線形補間では失われている高周波情報が復元されていることや，スパイクノイズ，ガウシアンノイズが超解像後には消えているのがわかります．

16枚の低解像度画像からの超解像処理結果

リサイズ（線形補間）

### “スパース行列(SparseMat)の使い方2：超解像” への17件のコメント

1. Wang Yongwen より:

Thank you for your sharing! I try to run your code on VS2010 and Open CV but there are some errors existing such as”Unhandled exception at 0x7603c6e3 in SuperResolution.exe: Microsoft C++ exception: cv::Exception at memory location 0x001be79c..”, which I don’t know how to debug it.Could you give me some suggestions about how to run this code successfully to achieve super_resolution?
Thank you very much!

2. Nguyen Le より:

Hi Norishige FUKUSHIMA,

If I already had my own sample images which were shifted a little in comparison with the first image, how can I modify your code to achieve the matrix A as well as the final super resolution image?
I already tried to use the function “calcOpticalFlowPyrLK” to calculate the displacement between the first image and the others based on the iterative Lucas-Kanade method with pyramids. Then, I returned the result of that function to the array move[i].x and move[i].y in order to calculate the matrix A. However, my solution did not seem to work correctly at all!!
Please let me know if you have any suggestion!!
Thanks in advance!!

• fukushima より:

Please modify the loop of (1) and the function of “createDegradedImageandSparseMat32″. This loop generates randomly shifted-multiple images from single image. If you want to use pre-captured image, please modify this part.
BRs,

3. Le Nguyen より:

Hi Norishige FUKUSHIMA,

If I already had my own sample images which were shifted a little in comparison with the first image, how can I modify your code to achieve the matrix A as well as the final super resolution image?
I already tried to use the function “calcOpticalFlowPyrLK” to calculate the displacement between the first image and the others based on the iterative Lucas-Kanade method with pyramids. Then, I returned the result of that function to the array move[i].x and move[i].y in order to calculate the matrix A. However, my solution did not seem to work correctly at all!!
Please let me know if you have any suggestion!!
Thanks in advance!!

4. I have a problem with this code.
What is the allowed maximum size of the input image?

• fukushima より:

I cannot know that because the size depends on OpenCV core function’s version deeply. In my computer, I cannot set large image size (256×256, 512×512).

5. iOn_azuma より:

haartrainingを使用しています。-nstages 21にて機械学習を行っていると処理に時間が掛かっているstage（だいたい17か18）の途中でtraing_data.xmlが作成されます。
CPU負荷と計算中の％表示から一見計算は続いている様に見えますが、この現象は各パラメータの設定が不適切で機械学習が最後まで実施できていない際の処理と考えています。
この場合の処置として正しいのは、
１．処理を途中で打ち切って再度パラメータを変更して作業を行う。
できたxmlは使えない。
２．まだ計算中なので継続して処理を行う。
最終的に　xmlが上書きされる。
作業としてどちらが正解か、ご存知の方いないでしょうか？
あれこれ試しておりますが、顔認識用サンプルのxmlのようにファイルができません。
よろしくお願いいたします。

6. Utkarsh より:

Awesome – thanks for the code. Do you know if any newer, better techniques have come up?

7. David より:

Dear Fukushima:

Sorry if this is not the appropiate place for a question, let me know
where to post if not here.

I tried to use your code in a set of images not generated by the progam itself.
These images are aligned; so no movement is taken in account. They are noisy already.
I can’t figure how to calculate the A matrix as createDegradedImageandSparseMat32F( ) does;
because createDownsampledMotionandBlurCCDSparseMat32f( ) expects a source of the same size as the output.
I tried resizing the input images to destination size (rfactor times), and using move parameter 0,0.
The noise was removed but no superresolution achieved.
Can you explain how to work with preloaded images or how to calculate the A matrix in the situation
of aligned pictures ?

Best regards, David

• fukushima より:

Dear David

Please use following function for computing sparse A matrix.
createDownsampledMotionandBlurCCDSparseMat32f(src,amp,move);
src: image which has desired (super resesoluved)resolution, which has not contain image data, bacause only image size is required.
amp:(high resolution)/(low resolution)
move:global motion per image.

>>The noise was removed but no superresolution achieved.
This implementation is based on multi-frame super resolution[1].
Thus, this method has weak capability for single image SR.

[1]Farsiu, S.,Robinson, D., Elad, M., Milanfar, P.”Fast and robust multiframe super resolution,” IEEETrans.ImageProcessing 13 (2004)1327-1344.

Regards,

Norishige FUKUSHIMA
in the main function (at line 74-76)
A[i]=createDegradedImageandSparseMat32F(image, imtemp,move[i],rfactor);
addgaussnoise(imtemp,degimage[i],10.0);//add gaussian noise
addspikenoise(degimage[i],degimage[i],500);//add spike noise

• David より:

Thank you !
David

8. 長澤亨 より:

すみません、お教え下さい。
プログラム中のimages1が未定義なのではありませんか？
コンパイルエラーがでますが、、、

void addgaussnoise(Mat& src, Mat& dest, double sigma) {

images1.convertTo(src_f,CV_32FC1); // <– images1が未定義？

}

• fukushima より:

コードをアップするときにソースが特殊文字扱いになってしまいimages1…と表示されている部分ですが，本来はimage[ c ]　という文字列が入っています．
コードにこの文字化け？変更が行われないようにアップし直しましたので確認してみてください．

9. PVDHP より:

All it’s cool! But where is dowload link and how to compile it?

10. fukushima より:

福嶋です．
このコードでは，スパース行列の入力部分の例外処理（294行目）で負の値が入った場合にうまく機能しないようになってしまっています。ここを修正すればプログラムが落ちることはなくなると思います．

しかし，行列を作る部分（300～303行目）では，第一象限（x,yともに正）の場合のみうまく行くように線形補間での行列アクセスを設計しています．

例えば平行移動が，x,y = (1.4,1.8)の場合，この画素のこの線形補間は，
0.6*0.2*[1,1] + 0.4*0.2*[2,1] + 0.6*0.8*[1,2]+ 0.4*0.8*[2,2]
で表されます．ここで，[m,n]は(m,n)pixelの画素の値です．
しかしここにマイナスが入って来た場合x,y = (-1.4,-1.8)
このプログラムでは，
[-1,-1],[-1,0],[0,-1],[0,0]
の４ペアを使って線形補間しようとします．
正しい最近傍は
[-2,-2],[-2,-1],[-1,-2],[-1,-1]
のペアとなるため，プログラムを修正する必要があると思います．

対策として，もっとも小さな（負の方向に絶対値が大きい）ものを基準画像に設定してしまうかこのプログラムを修正してみてください．

• r_takuma より:

何度もご親切なアドバイスありがとうございました．
福嶋さんのおかげで無事，負の値をとる場合でも超解像処理を行うことができました．

今回だけでなく，無知な私に優しく親切なアドバイスをしていただき，本当にありがとうございました．
私は，学生なのですが，福嶋さんのおかげで研究(といっても，それほど大したものでもありませんが．．．)を進めることができます．

お忙しいところ，本当にありがとうございました．

11. r_takuma より:

はじめまして．こちら(http://oshiete.goo.ne.jp/qa/6345132.html)でも何度か質問させていただいた者です．

今，こちらのプログラムを参考にさせていただいて，複数の低解像度画像から高解像度画像を生成するプログラムを作っています．
そこで，どうしても分からないところがあり，こちらに質問させていただきました．
コメント欄にこのような質問を投稿してしまう無礼をお許しください．

実行しようとしていることの概要としましては，
１，複数の低解像度画像からの高解像度画像生成
２，複数の低解像度画像は，平行移動のみと仮定(アフィン変換等はしていない)
３，別のプログラムで，移動量は推定済み
４，その他，ブラーカーネル等のパラメータはも推定済み

このような条件のもと，プログラムしたところ，移動量が正の値をとる場合の超解像処理は無事行うことができました．
しかし，移動量が負の値をとる場合，処理が最後まで実行できないバグが発生してしまいました．
訂正すべき箇所は，260-310行目のcreateDownsampledMotionandBlurCCDSparseMat32f内であることは理解できたのですが，この関数がどのような処理を行っているかが理解しかねているため，どのような訂正を行えば良いかが分かりません．

よろしければ，アドバイス等いただけたらと思います．
長文失礼しました．