README.md

MotionDeBlur

DIP final project . Doing a restore for motion blur

（使用時記得修改圖片讀檔的相對路徑）

---

>>> finalproject_psf <<<

實作PSF 解釋大概的演算法

第一個for迴圈是將圖檔讀進來

第二個for迴圈是建立點分散函式 這邊的位移量L是定30，可以從最上面的define去修改 而角度是45度，要改的話就要改if的y x之間關係 （ tan 三角函數 ，ex tan 45度 = 1 = y/x ）

第三個for迴圈就是實作psf PSF其實就是點分散的矩陣與原來圖檔的矩陣做convolution，就像做filter一樣 但有點不同的是，不能讓每一點相加後不處理， 所以我還多了一個變數number去做正規化的處理。

第四個迴圈將結果輸出成圖檔

final_project_psf_myself.cpp

//
//  main.cpp
//  finalproject-psf
//
//  Created by 陳昱安 on 13/5/25.
//  Copyright (c) 2013年 陳昱安. All rights reserved.
//

#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;

#define L 30
#define IMGL 512

int global_matrix[IMGL][IMGL]={0};
double output_matrix[IMGL][IMGL]={0};
int out_number[IMGL][IMGL] = {0} ;

int main(int argc, const char * argv[])
{
    fstream fin ;
    fin.open("/Users/mikeva/Chen/02NTU/大四/11數位影像處理/final project/lenna512.raw",ios::in);
    ofstream fout("/Users/mikeva/Chen/02NTU/大四/11數位影像處理/final project/result_myself_L_30.raw");
    char y;
    
    for(int i = 0 ; i < IMGL ; i++)
    {
        for(int j = 0 ; j < IMGL ; j++){
            fin.get(y);
            global_matrix[i][j]= ( (((int)y) + 256 )% 256 ) ;
        }
    }
    int ii = 0 ;
    int jj = 0 ;
    for(int t = 0 ; t < L ; t ++)
    {
        for(int i = 0 + t ; i < IMGL ; i ++){
            for(int j = 0 + t ; j < IMGL ; j ++ ){
                output_matrix[i][j] = output_matrix[i][j] + ( (double) global_matrix[ii][jj] ) ;
                out_number[i][j] ++ ;
                jj ++ ;
            }
            ii ++ ;
            jj = 0 ;
        }
        ii = 0 ;
        jj = 0 ;
    }
    for(int i = 0 ; i < IMGL ; i++)
    {
        for(int j = 0 ; j < IMGL ; j ++)
        {
            fout << (unsigned char) ( (output_matrix[i][j]) / out_number[i][j] );
        }
    }
    cout << " >> finish ... " <<endl;
    return 0;
}

finalproject_psf.cpp

//
//  main.cpp
//  finalproject-psf
//
//  Created by 陳昱安 on 13/5/25.
//  Copyright (c) 2013年 陳昱安. All rights reserved.
//

#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;

#define L 40
#define IMGL 512

int global_matrix[IMGL][IMGL]={0};
double h_matrix[IMGL][IMGL]={0};
double output_matrix[IMGL][IMGL]={0};

double PSF(int a[IMGL][IMGL],double b[IMGL][IMGL],int i,int j);

int main(int argc, const char * argv[])
{
    fstream fin ;
    fin.open("/Users/mikeva/Chen/02NTU/大四/11數位影像處理/final project/lenna512.raw",ios::in);
    ofstream fout("/Users/mikeva/Chen/02NTU/大四/11數位影像處理/final project/result_L_40.raw");
    char y;
    
    for(int i = 0 ; i < IMGL ; i++)
    {
        for(int j = 0 ; j < IMGL ; j++){
            fin.get(y);
            global_matrix[i][j]= ( (((int)y) + 256 )% 256 ) ;
        }
    }
    for(int i = 0 ; i < IMGL ; i++)
    {
        for(int j = 0 ; j < IMGL ; j ++)
        {
//            cout << j <<endl;

            if( sqrt( (double)(i*i+ j*j) )  <= L && i==j ){
                h_matrix[i][j] =  (double) 1  ;
                    }
            else{
                h_matrix[i][j] = 0 ;
            }
        }
    }

    for(int i = 0 ; i < IMGL ; i++)
    {
        for(int j = 0 ; j < IMGL ; j ++)
        {
            output_matrix[i][j] = PSF(global_matrix, h_matrix,i,j);
           
        }
         cout << i <<endl;
    }
    
    for(int i = 0 ; i < IMGL ; i++)
    {
        for(int j = 0 ; j < IMGL ; j ++)
        {
            fout << (unsigned char) (output_matrix[i][j]);
//            cout <<i << "<< "<< j << " : " << output_matrix[i][j] ;
        }
    }
    
    
    return 0;
}


double PSF(int a[IMGL][IMGL],double b[IMGL][IMGL],int i,int j){
    
    double result = 0.0;
    int h_i = 0 ;
    int h_j = 0 ;
    int number = 0 ;
    for(int ii = i- (L/2) ; ii < i + (L/2) ; ii ++ ){
        for(int jj = j - (L/2) ; jj < j + (L/2)  ; jj ++ ){
            if(jj >= 0  && jj < IMGL && ii >= 0 && ii < IMGL)
            {
                result = result + (double) a[ii][jj]* b[h_i][h_j]  ;
                if( b[h_i][h_j] > 0 ){
                    number ++ ;
                }
                
            }
//            cout << jj <<endl;
            h_j ++ ;
        }
        h_i ++ ;
    }
//    cout << "nu "<< number << endl;
    return ( result / number );
}

finalproject_restore.cpp

//
//  main.cpp
//  finalproject_restore
//
//  Created by 陳昱安 on 13/5/25.
//  Copyright (c) 2013年 陳昱安. All rights reserved.
//
#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;

#define L 40
#define IMGL 512

int global_matrix[IMGL][IMGL]={0};
double f_matrix[IMGL][IMGL]={0};
double f_matrix_n[IMGL][IMGL]={0};
double h_matrix[IMGL][IMGL]={0};
double h_matrix_n[IMGL][IMGL]={0};
double h_adjoint[IMGL][IMGL]={0};
double execute_temp[IMGL][IMGL]={0};


double adjoint(double m,int i , int j );
double execute(int g[IMGL][IMGL] ,double h[IMGL][IMGL],double f[IMGL][IMGL],double h_a[IMGL][IMGL],int i,int j );
double PSF(double a[IMGL][IMGL],double b[IMGL][IMGL],int i,int j);
//double execute1(int g[IMGL][IMGL] ,double h[IMGL][IMGL],double f[IMGL][IMGL],int i , int j );

int main(int argc, const char * argv[])
{
    fstream fin ;
    fin.open("/Users/mikeva/Chen/02NTU/大四/11數位影像處理/final project/result_L_40.raw",ios::in);
    ofstream fout("/Users/mikeva/Chen/02NTU/大四/11數位影像處理/final project/result_restore_L_40.raw");
    char y;
    cout<<" >> initialize ... "<<endl;
    // initialize 
    for(int i = 0 ; i < IMGL ; i++)
    {
        for(int j = 0 ; j < IMGL ; j++){
            fin.get(y);
            global_matrix[i][j]= ( (((int)y) + 256 )% 256 ) ;
        }
    }

    for(int i = 0 ; i < IMGL ; i++)
    {
        for(int j = 0 ; j < IMGL ; j ++)
        {
            //            cout << j <<endl;
            
            if( sqrt( (double)(i*i+ j*j) )  <= L && i==j ){
                h_matrix[i][j] =  (double) 1  ;
            }
            else{
                h_matrix[i][j] = 0 ;
            }
        }
    }
    for(int i = 0 ; i < IMGL ; i++)
    {
        for(int j = 0 ; j < IMGL ; j++){
            h_adjoint[i][j] = adjoint(h_matrix[j][i], (i+1), (j+1));
            cout <<i<<"  "<<j <<":"<<h_adjoint[i][j]<<endl;
        }
    }
    
    for(int i = 0 ; i < IMGL ; i++)
    {
        for(int j = 0 ; j < IMGL ; j++){
            f_matrix[i][j] =(double) global_matrix[i][j] ;
        }
    }

    cout<<" >> executing ... "<<endl;
    //  開始實作
    for(int i = 0 ; i < IMGL ; i++)
    {
        for(int j = 0 ; j < IMGL ; j++){
            execute_temp[i][j] = execute(global_matrix, h_matrix, f_matrix,h_adjoint, i, j);
//            cout << execute_temp[i][j]<<endl;
        }
    }
    
    for(int i = 0 ; i < IMGL ; i++)
    {
        for(int j = 0 ; j < IMGL ; j++){
            f_matrix_n[i][j] = PSF(execute_temp, h_adjoint, i, j);
//            cout << f_matrix_n[i][j]<<endl;
            f_matrix_n[i][j] = PSF(f_matrix,f_matrix_n,i,j)  ;
           
        }
    }
    for(int i = 0 ; i < IMGL ; i++)
    {
        for(int j = 0 ; j < IMGL ; j++){
            fout << (unsigned char) f_matrix_n[i][j];
        }
    }
    
    cout<<" >> finishing ... "<<endl;
    return 0;
}

// function area
double adjoint(double m,int i , int j ){
    double result;
    
    result = pow( (-1),(i+j) ) * m ;
    
    return result;
}


double execute(int g[IMGL][IMGL] ,double h[IMGL][IMGL],double f[IMGL][IMGL],double h_a[IMGL][IMGL],int i,int j ){

    double result = 0.0 ;
    double temp = 0.0 ;
//    double temp2 = 0.0;
    temp = (double) g[i][j] / PSF(f,h,i,j);
//    temp2 = PSF(f, h_a, i, j);
    
//    result = temp * temp2 ;
    result = temp;
    return result;
}


//double execute1(int g[IMGL][IMGL] ,double h[IMGL][IMGL],double f[IMGL][IMGL],double h_a[IMGL][IMGL],int i,int j ){
//    
//    double result = 0.0 ;
//    double temp = 0.0 ;
//    temp = (double) g[i][j] / PSF(f,h,i,j);
//     
//    result = temp ;
//    return result;
//}




double PSF(double a[IMGL][IMGL],double b[IMGL][IMGL],int i,int j){
    
    double result = 0.0;
    int h_i = 0 ;
    int h_j = 0 ;
    int number = 0 ;
    for(int ii = i- (L/2) ; ii < i + (L/2) ; ii ++ ){
        for(int jj = j - (L/2) ; jj < j + (L/2)  ; jj ++ ){
            if(jj >= 0  && jj < IMGL && ii >= 0 && ii < IMGL)
            {
                result = result + (double) a[ii][jj]* b[h_i][h_j]  ;
                if( b[h_i][h_j] > 0 ){
                    number ++ ;
                }
                
            }
            //            cout << jj <<endl;
            h_j ++ ;
        }
        h_i ++ ;
    }
    //    cout << "nu "<< number << endl;
    return ( result / number );
}