平滑化・エッジ抽出フィルタの高速化

[loss]

画像の平滑化およびエッジ抽出のプログラムを書きました。
この処理内容を変えずに、消費時間の短縮(高速化)を行うためにはどのような改善が必要でしょうか。
改善内容とその理由も明記していただけると助かります。
複数パターンの改善法やその短縮度合いも比較してみたいため、様々な方法でのお力添えをいただきたいです。

コード:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

typedef unsigned char UCHAR;

				/* 画像サイズは変更不可 */
#define DIM1  960
#define DIM2 1280
#define DIM3    3

#define N_REPEAT 100		/* 繰返し回数．消費時間測定時は必ず100とする． */

void write_ppm_cip(UCHAR [][DIM2][DIM3],char *,int ,int );
void read_ppm_cip(UCHAR [][DIM2][DIM3], char *);
void filtering(UCHAR [][DIM2][DIM3], UCHAR [][DIM2][DIM3], double [], int, int);

UCHAR org[DIM1][DIM2][DIM3];	/* 原画像(入力画像) */
UCHAR res1[DIM1][DIM2][DIM3];	/* 結果画像１(出力画像) */
UCHAR res2[DIM1][DIM2][DIM3];	/* 結果画像２(出力画像) */

int main(int argc, char *argv[])
{
  int i;
  char in_fname[200];
  char out_fname[200];
  double kernel[][9]={
    {1/9.0,1/9.0,1/9.0,1/9.0,1/9.0,1/9.0,1/9.0,1/9.0,1/9.0}, /* 平滑化のカーネル（係数） */
    {-1,0,1,-2,0,2,-1,0,1}	/* 横エッジ抽出(sobelフィルタ)のカーネル（係数） */
  };
  
  if(argc < 2)
  {
    printf("usage: a.out filename\n");
    exit(1);
  }

  strcpy(in_fname,argv[argc-1]); /* 入力画像のファイル名 */

  read_ppm_cip(org, in_fname);	/* 画像をファイルより読み込み */

				/* フィルタリング(時間測定のため必ず100回行う．) */
  for(i=0;i<N_REPEAT;i++)
  {
    filtering(res1,org,kernel[0],DIM1,DIM2); /* 平滑化 */
    filtering(res2,res1,kernel[1],DIM1,DIM2); /* 横エッジ抽出 */
  }  
				/* 出力ファイル名設定 */
  in_fname[strlen(in_fname)-4]='\0'; /* .ppmを除く */
  sprintf(out_fname,"%s_hei.ppm",in_fname);	/* 平滑化画像 */
  write_ppm_cip(res1,out_fname,DIM2,DIM1); /* 画像をファイルへ出力 */

  sprintf(out_fname,"%s_edge.ppm",in_fname);	/* 横エッジ抽出画像 */
  write_ppm_cip(res2,out_fname,DIM2,DIM1); /* 画像をファイルへ出力 */

  return 0;
}


/* フィルタリング処理 */
/* res1, res2は静的変数のため初期値は０．従って画像の最外周の画素値は０． */
void filtering(UCHAR ans[][DIM2][DIM3], UCHAR data[][DIM2][DIM3], double a[], int n_gyou, int n_retu)
{
  int g,r,k;
  int n;
  int i,j;
  double tmp;
  
  for(g=1;g<n_gyou-1;g++)
    for(r=1;r<n_retu-1;r++)
      for(k=0;k<DIM3;k++)
      {
	tmp=0.0;
	n=0;
	for(i=-1;i<=1;i++)
	  for(j=-1;j<=1;j++)
	  {
	    tmp+=data[g+i][r+j][k]*a[n];
	    n++;
	  }
	if(tmp < 0)		/* 絶対値をとる */
	  tmp=-tmp;
	if(tmp > 255)
	  ans[g][r][k]=255;
	else
	  ans[g][r][k]=(UCHAR)(tmp+0.5);
      }
}


#define W_BYTE DIM1*DIM2*DIM3

/* CIP形式の画像データよりPPMファイル作成 */
void write_ppm_cip(UCHAR data_buf[][DIM2][DIM3],char *fname,int width,int height)
{
  FILE *fp;
  
				/* ファイルを開く */
  if((fp = fopen(fname, "wb")) == NULL) {
    fprintf(stderr, "file(%s) can't open\n", fname) ;
    exit(1) ;
  }

  fprintf(fp, "P6\n") ;		/* カラー画像かつバイナリーデータの記号 */
  fprintf(fp, "%d %d\n", width, height) ; /* 画像の幅(列数)と高さ(行数) */
  fprintf(fp, "255\n") ;	/* 最大値 */

  {				/* 画像データをrepeat(+1)個に分割して書き込み */
    int i;
    int repeat=(DIM1*DIM2*DIM3)/(W_BYTE); /* 分割数 */
    int rest=DIM1*DIM2*DIM3-repeat*(W_BYTE); /* 余りデータ量 */
    UCHAR *pt=(UCHAR *)&data_buf[0][0][0]; /* 書き込むデータの位置を持つポインタ */

    for(i=0;i<repeat;i++)
    {
      fwrite(pt, sizeof(UCHAR), W_BYTE, fp);
      pt += W_BYTE;
    }
    if(rest > 0)
      fwrite(pt, sizeof(UCHAR), rest, fp);
  }
  
  fclose(fp);			/* ファイルを閉じる */    
}


#define R_BYTE DIM1*DIM2*DIM3

/* ppmフォーマットの画像ファイルをCIP形式の画像メモリに読み込み */
void read_ppm_cip(UCHAR data_buf[][DIM2][DIM3], char *fname)
{
  FILE	*fp ;
  char	str_buf[1024] ;
  char	magic_num[8] ;		/* マジックナンバー */
  int	max_val ;		/* 画素値の最大値 */    
  int width,height;
  int count_limit;

				/* ファイルを開く */    
  if((fp = fopen(fname, "r")) == NULL) {
    fprintf(stderr, "file(%s) can't open.\n", fname) ;
    exit(1) ;
  }
				/* マジックナンバー読み込み */
  count_limit=0;
  sprintf(str_buf,"#");
  while((str_buf[0]=='#') || (str_buf[0]=='\n') || (str_buf[0]=='\t') || (str_buf[0]==' '))
  {
    fgets(str_buf,1024,fp);
    count_limit++;
    if(count_limit > 1000)
    {
      fprintf(stderr,"ERROR: Irregal file format.\n");
      exit(1);
    }
  }
  strcpy(magic_num, str_buf);
  magic_num[strlen(magic_num)-1]='\0'; /* \nを除く */
  if(strcmp(magic_num, "P6") != 0) 
  {
    fprintf(stderr, "ERROR: magic number(%s) not match.\n", magic_num) ;
    exit(1) ;
  }

				/* 画像の幅(列数)と高さ(行数) */
  count_limit=0;
  sprintf(str_buf,"#");
  while((str_buf[0]=='#') || (str_buf[0]=='\n') || (str_buf[0]=='\t') || (str_buf[0]==' '))
  {
    fgets(str_buf,1024,fp);
    count_limit++;
    if(count_limit > 1000)
    {
      fprintf(stderr,"ERROR: Irregal file format.\n");
      exit(1);
    }
  }
  sscanf(str_buf,"%d %d",&width,&height);
  if((height != DIM1) || (width != DIM2)) /* 画像のサイズが仮定した値と異なる場合はエラーで終了 */
  {
    fprintf(stderr, "ERROR: Dimension dosenot match.\n");
    exit(1) ;
  }

				/* 最大値 */
  count_limit=0;
  sprintf(str_buf,"#");
  while((str_buf[0]=='#') || (str_buf[0]=='\n') || (str_buf[0]=='\t') || (str_buf[0]==' '))
  {
    fgets(str_buf,1024,fp);
    count_limit++;
    if(count_limit > 1000)
    {
      fprintf(stderr,"ERROR: Irregal file format.\n");
      exit(1);
    }
  }
  sscanf(str_buf,"%d",&max_val);
  if(max_val > 255)
  {
    fprintf(stderr, "ERROR: Irregal max value.\n");
    exit(1) ;
  }

				/* 画像データをrepeat(+1)個に分割して読み込み */
  {
    int i;
    int repeat=(DIM1*DIM2*DIM3)/(R_BYTE); /* 分割数 */
    int rest=DIM1*DIM2*DIM3-repeat*(R_BYTE); /* 余りデータ量 */
    UCHAR *pt=(UCHAR *)&data_buf[0][0][0]; /* 読み込みデータの位置を持つポインタ */

    for(i=0;i<repeat;i++)
    {
      fread(pt, sizeof(UCHAR), R_BYTE, fp);
      pt += R_BYTE;
    }
    if(rest > 0)
      fread(pt, sizeof(UCHAR), rest, fp);
  }

  /* ファイルを閉じる */    
  fclose(fp);
}

[h2so5]

平滑化の処理については係数がすべて同じなので、分配法則を用いてピクセルあたり9回の乗算を1回に減らせると思います。

[loss]

ありがとうございます。
できれば、どのような書き方をすると効率的かも具体的に教えてくださると助かります。参考にしたいです。

[usao]

そもそも，フィルタカーネルの内容に対して特化した話でもよいのか，
フィルタカーネル内容によらない汎用関数filtering()を高速化したい　という話なのか，どっちなんでしょう？

前者であれば，「平滑化なら～」という話ができますが，
後者だとすれば，まず，現状だと何が（どこが原因で）遅いの？ということを示してもらわないと答えられないんじゃないかなぁ？
（入力や出力の配列へのアクセスが遅いとかいう話なら→なるべく配列アクセス回数を減らす　とかいう方針が立つかもしれない）