USER
// 次のテンプレートを使って、 // 画像処理の学習、アルゴリズムの実証・検証作業をアシストしてください。 // まずは、生成AIに、このテンプレートの意味と使い方を説明をお願いします。 ///* gcc -O3 -std=c99 -mavx -W -Wall -Wextra -o Template Template.c -fopenmp 環境依存のコードが生成されますが、速度重視の場合は検討してください。 gcc -O3 -std=c99 -march=native -mavx -W -Wall -Wextra -o Template Template.c -fopenmp 画像処理の学習、アルゴリズム検証用のフレームワーク、 C99 テンプレート 目的は、アルゴリズムの検証を簡単に行えるようにする事と、標準では外部ライブラリ （例外、stb_image.h）を使用しない環境を提供し、環境準備の構築工数を低減させる。 A professionally created C99 template for learning and verifying image processing The purpose is to make it easy to verify algorithms and to provide an environment that does not use external libraries (exception: stb_image.h) by default, reducing the labor required to set up the environment. // // // You have the right to: You are free to use this source code and modify it // // for your own use without any restrictions. In addition, we do not provide // // any warranty or support under any circumstances. // // // // Note: OpenMP 3.0 以降では、#pragma omp parallel for のループ変数は、自動的に // // プライベートになるので宣言は省略できるが、明示した方が誤解を回避できる // // // */ #include <stdio.h> #include <stdlib.h> #include <stdint.h> //////// https://github.com/nothings/stb/blob/master/stb_image.h #define STB_IMAGE_IMPLEMENTATION #include "stb_image.h" #include <omp.h> #define HERE fprintf(stderr,"% 6d: %s\n",__LINE__, __func__); #define IF if #define EI else if #define EE else ////////////////////////////////////////////////////////////////////////////////////////// // Timer section #ifdef _WIN32 #include <windows.h> double get_micro_time() { static LARGE_INTEGER freq; LARGE_INTEGER now; if (freq.QuadPart == 0) QueryPerformanceFrequency(&freq); QueryPerformanceCounter(&now); return (double)now.QuadPart * 1e6 / (double)freq.QuadPart; } #else #include <sys/time.h> double get_micro_time() { struct timeval tv; gettimeofday(&tv, NULL); return (double)tv.tv_sec * 1e6 + (double)tv.tv_usec; } #endif // double START_time = get_micro_time(); // double E_N_D_time = get_micro_time(); // double Elapsed_time = E_N_D_time - START_time; // printf_elapsed_time(Elapsed_time); printf_elapsed_time(E_N_D_time-START_time); void print_elapsed_time(double microsec) { IF(microsec < 1000.0) printf("Elapsed time: %.3f μsec\n", microsec); EI(microsec < 1e6 ) printf("Elapsed time: %.3f msec\n", microsec / 1000.0); EE printf("Elapsed time: %.3f sec\n", microsec / 1e6); } //////////////////////////////////////////////////////////////////////////////// // Basic structure definition // #pragma pack(push,1) typedef struct PIXEL { unsigned char B, G, R, A; } PIXEL ; typedef struct STBPX { unsigned char R, G, B, A; } STBPX ; #pragma pack(pop) // BMPヘッダ構造体 #pragma pack(push, 2) typedef struct HEADER { unsigned short bfType; // 'BM' (0x4D42) unsigned int bfSize; // ファイルサイズ unsigned short bfReserved1; // 0 unsigned short bfReserved2; // 0 unsigned int bfOffBits; // データオフセット unsigned int biSize; // 情報ヘッダサイズ (40) signed int biWidth; // 幅 signed int biHeight; // 高さ（負で上下反転） unsigned short biPlanes; // 1 unsigned short biBitCount; // 32 unsigned int biCompression; // 0 (BI_RGB) unsigned int biSizeImage; // 画像データサイズ signed int biXPelsPerMeter; // 0 signed int biYPelsPerMeter; // 0 unsigned int biClrUsed; // 0 unsigned int biClrImportant; // 0 unsigned char lut[]; // パレット（8ビット用） } HEADER; #pragma pack(pop) // 画像パラメータ #define BMPID 0x4D42 // 'BM' #define BMPIFHSIZE 40 // 情報ヘッダサイズ #define v2D(type,lx,name) ((type(*)[(lx)])name) #define v2Dp(lx,name) v2D(PIXEL,lx,name) #define v2Dc(lx,name) v2D(char,lx,name) #define v2Di(lx,name) v2D(int,lx,name) #define v2Df(lx,name) v2D(float,lx,name) #define v2Dd(lx,name) v2D(double,lx,name) #define saveArray2bmp32(name,a,w,h) saveDataAsBMP32((name),(&(a)[0][0]),(w),(h)) #define PROCESS_ALL_PIXELS for(int y=0;y<ly;y++)for(int x=0;x<lx;x++) //////////////////////////////////////////////////////////////////////////////// // Global variables int verbose=0, auto_adjust=0, max_num_threads=0 ; //////////////////////////////////////////////////////////////////////////////// // Primitive functions unsigned int int2PseudoColor(unsigned int v) { register unsigned int r=0, g=0, b=0 ; if( (v&0x00ffffff) == 0 ) return 0 ; if(v & 0x000001) b|=0x80;; if(v & 0x000002) g|=0x80;; if(v & 0x000004) r|=0x80; if(v & 0x000008) b|=0x40;; if(v & 0x000010) g|=0x40;; if(v & 0x000020) r|=0x40; if(v & 0x000040) b|=0x20;; if(v & 0x000080) g|=0x20;; if(v & 0x000100) r|=0x20; if(v & 0x000200) b|=0x10;; if(v & 0x000400) g|=0x10;; if(v & 0x000800) r|=0x10; if(v & 0x001000) b|=0x08;; if(v & 0x002000) g|=0x08;; if(v & 0x004000) r|=0x08; if(v & 0x008000) b|=0x04;; if(v & 0x010000) g|=0x04;; if(v & 0x020000) r|=0x04; if(v & 0x040000) b|=0x02;; if(v & 0x080000) g|=0x02;; if(v & 0x100000) r|=0x02; if(v & 0x200000) b|=0x01;; if(v & 0x400000) g|=0x01;; if(v & 0x800000) r|=0x01; return (r <<16) | (g << 8) | b ; } void **alloc2D( void *data, int s, int lx, int ly ){ long long i, w=lx, h=ly; if( !data || s<1 || lx<1 || ly<1 ) return NULL; void **d, *base=data ; if( (d=(void **)malloc(h*sizeof(void *))) == NULL ) return NULL ; for( i=0; i<ly; i++ ) d[i] = (void *)((unsigned char *)base + w*s*i) ; return d ; } void **alloc2Darray( int s, int lx, int ly ){ long long i, w=lx, h=ly; if(s<1 || lx<1 || ly<1) return NULL ; void **d, *base ; if( (d=(void **)calloc(1,h*sizeof(void *)+(w*h*s))) == NULL ) return NULL ; base = (unsigned char *)d + h*sizeof(void *) ; for( i=0; i<ly; i++ ) d[i] = (void *)((unsigned char *)base + w*s*i) ; return d ; } void int2PIXEL( int **s, int **d, int lx, int ly ){ if( !s || !d || lx<1 || ly<1) return ; PIXEL **pd = (PIXEL**)d; for(int value, y=0;y<ly;y++)for(int x=0;x<lx;x++){ value = s[y][x]<0 ? 0 : ( s[y][x]>255 ? 255 : s[y][x] ); pd[y][x] = (PIXEL){ value, value, value, 0 }; } } void int2PIXEL_adjust( int **s, int **d, int lx, int ly ){ if( !s || !d || lx<1 || ly<1) return ; int min=s[0][0], max=s[0][0]; PROCESS_ALL_PIXELS IF(min>s[y][x])min=s[y][x]; EI(max<s[y][x])max=s[y][x]; int div = min!=max ? (max-min) : 1 ; // divide number PIXEL **pd = (PIXEL**)d; for(int value, y=0;y<ly;y++)for(int x=0;x<lx;x++){ value = 255*s[y][x]/div; pd[y][x] = (PIXEL){ value, value, value, 0 }; } } int saveDataAsBMP32( char *fname, void *b, int lx, int ly ){ HEADER h = {BMPID, lx*abs(ly)*sizeof(PIXEL)+sizeof(HEADER), 0, 0, sizeof(HEADER), BMPIFHSIZE, lx, ly, 1, 32, 0, lx*abs(ly)*sizeof(PIXEL), 0, 0, 0, 0 }; FILE *fp=NULL ; int st=-1; if( (fp=fopen(fname, "wb")) == NULL ) { fprintf( stderr, "File '%s': cannot create.\n", fname ); } else if ( fwrite( &h, sizeof(HEADER), 1, fp ) != 1 ) { fprintf( stderr, "File '%s': Header write error.\n", fname ); } else if ( fwrite( b, h.biSizeImage, 1, fp ) != 1 ) { fprintf( stderr, "File '%s': Data write error.\n", fname ); } else st=0; if(fp)fclose(fp) ; return st; } ////////////////////////////////////////////////////////////////////////////////////////// // Function set for algorithm verification // // convolution // float版 畳み込み関数（ゼロパディング） void convolution_float(float **s, float **d, int lx, int ly, int n, float kernel[n][n]) { int u=n, v=n, ku = u / 2, kv = v / 2; #ifdef _OPENMP #pragma omp parallel for #endif for (int y = 0; y < ly; y++) { for (int x = 0; x < lx; x++) { float sum = 0.0f; for (int j = 0; j < u; j++) { for (int i = 0; i < v; i++) { int yy = y + j - ku; int xx = x + i - kv; if (xx >= 0 && xx < lx && yy >= 0 && yy < ly) { sum += s[yy][xx] * kernel[j][i]; } } } d[y][x] = sum; } } } void convolution(int **s, int **d, int lx, int ly, int n, int kernel[n][n]) { int u=n, v=n, ku = u / 2, kv = v / 2; // カーネル中心のオフセット #ifdef _OPENMP #pragma omp parallel for #endif for (int y = 0; y < ly; y++) { for (int x = 0; x < lx; x++) { int sum = 0; for (int j = 0; j < u; j++) { for (int i = 0; i < v; i++) { int yy = y + j - ku; int xx = x + i - kv; // ゼロパディング：境界外は0とみなす if (xx >= 0 && xx < lx && yy >= 0 && yy < ly) { sum += s[yy][xx] * kernel[j][i]; } } } d[y][x] = sum; } } } void image_PIXEL2rgb(int **a, int **r, int **g, int **b, int lx, int ly) { PIXEL **p=(PIXEL**)a; PROCESS_ALL_PIXELS{r[y][x]=p[y][x].R; g[y][x]=p[y][x].G; b[y][x]=p[y][x].B;} } void image_rgb2PIXEL(int **r, int **g, int **b, int **a, int lx, int ly) { PIXEL **p=(PIXEL**)a; PROCESS_ALL_PIXELS{ int vr = r[y][x]<0 ? 0 : ( r[y][x]>255 ? 255 : r[y][x] ); int vg = g[y][x]<0 ? 0 : ( g[y][x]>255 ? 255 : g[y][x] ); int vb = b[y][x]<0 ? 0 : ( b[y][x]>255 ? 255 : b[y][x] ); p[y][x] = (PIXEL){ vb, vg, vr, 0 }; } } void rgb2PIXEL_adjust( int **r, int **g, int **b, int **a, int lx, int ly ){ int rmin=r[0][0], rmax=r[0][0], gmin=g[0][0], gmax=g[0][0], bmin=b[0][0], bmax=b[0][0]; PROCESS_ALL_PIXELS{ IF(rmin>r[y][x])rmin=r[y][x]; EI(rmax<r[y][x])rmax=r[y][x]; IF(gmin>g[y][x])gmin=g[y][x]; EI(gmax<g[y][x])gmax=g[y][x]; IF(bmin>b[y][x])bmin=b[y][x]; EI(bmax<b[y][x])bmax=b[y][x]; } int min, max, div; min = (rmin<gmin) ? (rmin<bmin?rmin:bmin) : (gmin<bmin?gmin:bmin); max = (rmax>gmax) ? (rmax>bmax?rmax:bmax) : (gmax>bmax?gmax:bmax); div = min!=max ? (max-min) : 1 ; // divide number PIXEL **p = (PIXEL**)a; PROCESS_ALL_PIXELS{// (PIXEL){ B, G, R, 0 }; p[y][x] = (PIXEL){ 255*b[y][x]/div, 255*g[y][x]/div, 255*r[y][x]/div, 0 }; } } void image_add(int **a, int **b, int **dst, int lx, int ly) { for(int y=0;y<ly;y++)for(int x=0;x<lx;x++) dst[y][x] = a[y][x] + b[y][x]; } void image_sub(int **a, int **b, int **dst, int lx, int ly) { for(int y=0;y<ly;y++)for(int x=0;x<lx;x++) dst[y][x] = a[y][x] - b[y][x]; } void image_mul(int **a, int **b, int **dst, int lx, int ly) { for(int y=0;y<ly;y++)for(int x=0;x<lx;x++) dst[y][x] = a[y][x] * b[y][x]; } void image_div(int **a, int **b, int **dst, int lx, int ly) { for(int y=0;y<ly;y++)for(int x=0;x<lx;x++) dst[y][x] = (b[y][x] != 0) ? a[y][x] / b[y][x] : 0; } void image_add_abs3(int **a, int **b, int **c, int **dst, int lx, int ly) { for(int y=0;y<ly;y++)for(int x=0;x<lx;x++) dst[y][x] = abs(a[y][x]) + abs(b[y][x]) + abs(c[y][x]); } void image_and(int **a, int **b, int **dst, int lx, int ly) { for(int y=0;y<ly;y++)for(int x=0;x<lx;x++) dst[y][x] = a[y][x] & b[y][x]; } void image_or(int **a, int **b, int **dst, int lx, int ly) { for(int y=0;y<ly;y++)for(int x=0;x<lx;x++) dst[y][x] = a[y][x] | b[y][x]; } void image_xor(int **a, int **b, int **dst, int lx, int ly) { for(int y=0;y<ly;y++)for(int x=0;x<lx;x++) dst[y][x] = a[y][x] ^ b[y][x]; } void image_not(int **a, int **dst, int lx, int ly) { for(int y=0;y<ly;y++)for(int x=0;x<lx;x++) dst[y][x] = ~a[y][x]; } void clamp_and_convert(int **src, PIXEL **dst, int lx, int ly) { for (int y = 0; y < ly; y++) for (int x = 0; x < lx; x++) { int v = src[y][x]; if (v < 0) v = 0; else if (v > 255) v = 255; dst[y][x] = (PIXEL){v, v, v, 0}; } } ////////////////////////////////////////////////////////////////////////////////////////// // Algorithm implementation area, start here // int save_without_conversion = 0; int kernel = 5; float sigma = 1.0;int mage_processing_program( int **s, int **d, int **w, int **t, int **g, int **b, int lx, int ly ){ (void)s; (void)d; (void)w; (void)t; (void)g; (void)b; (void)lx; (void)ly; ////////////////////////////////////////////////////////////// // image_processing_program // - s: 入力画像（カラー画像 or ネガ画像：輝度反転画像）BMP互換のPIXEL形式 // - d: 出力画像（処理結果） // - w, t, g, b: 作業用ワーク領域（必要に応じて自由に使用可） // t: カラー画像から、グレー画像データに変換したBMPファイル形式データ // g: カラー画像から、グレー画像データに変換した値、演算処理用 // b: グレー画像データを更に、128以上、以下で 255,0 に二値化したデータ /* 汎用的な作業エリアしか用意していませんので、 必要な作業領域は自己で確保、使用後は破棄を行ってください。 Ex. float **f = (float**)alloc2Darray(sizeof(float),lx,ly)); : free(f); */ //// Dummy functions for template purposes //// /* If you want to create a multi-threaded program #ifdef _OPENMP if(verbose) printf("Using %d threads.\n", omp_get_max_threads()); #endif Ex. #ifdef _OPENMP #pragma omp parallel for collapse(2) #endif for(int y=0;y<ly;y++) for(int x=0;x<lx;x++) d[y][x] = (s[y][x] + g[y][x]) / 2; */ /* 実際の画像処理（Sobelフィルタ）の実装例 int sobelX[3][3] = { {-1, 0, 1}, {-2, 0, 2}, {-1, 0, 1} }; int sobelY[3][3] = { {-1, -2, -1}, { 0, 0, 0}, { 1, 2, 1} }; convolution(g, w, lx, ly, 3, sobelX); // Gx convolution(g, b, lx, ly, 3, sobelY); // Gy image_add_abs3(w, b, b, d, lx, ly); // G = |Gx| + |Gy| */ return 0; } // Algorithm implementation area. ends here ////////////////////////////////////////////////////////////////////////////////////////// #include <string.h> void argumentAnalysis(int argc, char *argv[], char **iname, char **oname){ for(int i=1; i<argc; i++ ){ // プログラム起動時の引数解析処理 if( argv[i][0]=='-' ) { switch(argv[i][1]){ case 'v': if( argv[i][2] ) goto GOOD_BY; verbose=1; break; case 'V': if( argv[i][2] ) goto GOOD_BY; verbose=2; break; case 'a': if( argv[i][2] ) goto GOOD_BY; auto_adjust=1; break; // If you need for the multi-threaded program. omp_set_num_threads(n); case 't': if( !isdigit(argv[i][2]) || argv[i][2]==0 ) goto GOOD_BY; max_num_threads = atoi(&argv[i][2]); if(max_num_threads > 0) omp_set_num_threads(max_num_threads); break; //case 'k': if( strcmp(argv[i],"-kernel") && strcmp(argv[i],"-k") ) goto GOOD_BY; // if( ++i>=argc ) goto GOOD_BY; // kernel = atoi(argv[i]); break; //case 's': if( strcmp(argv[i],"-sigma") && strcmp(argv[i],"-s") ) goto GOOD_BY; // if( ++i>=argc ) goto GOOD_BY; // sigma = atof(argv[i]); break; case 'h': goto GOOD_BY; default : goto GOOD_BY; } } else if( *iname==NULL ) *iname = argv[i]; else if( *oname==NULL ) *oname = argv[i]; else goto GOOD_BY; } return; GOOD_BY:; int i; char *p=argv[0], *op="[-v|V] [-a] [-t<num>] [-h]"; for(i=strlen(p); i>0 && p[i]!='\\' && p[i]!='/'; i-- );; if(i) i++; fprintf( stderr, "\nusage: %s %s [input_bmp [output_bmp]]\n\n", &p[i], op ); fprintf( stderr, "\t\t-v: Verbose mode, V shows more details.\n" ); fprintf( stderr, "\t\t-a: Adjust the integer data range to be 0 to 255.\n" ); fprintf( stderr, "\t\t-t: -t<num>. a positive integer\n" ); fprintf( stderr, "\t\t Specify the maximum number of threads to use.\n" ); fprintf( stderr, "\t\t-h: Display this message and exit.\n" ); exit(0); } void PIXEL_Fitting( PIXEL *s, int lx, int ly ) { int x, y, w, ey=(ly +1)/2, ty=ly-1; #undef A_stb #define A_stb (v2Dp(lx,s)) for( y=0; y<ey; y++ ) for( w=ty-y, x=0; x<lx; x++ ) { PIXEL p = (PIXEL){ A_stb[y][x].R, A_stb[y][x].G, A_stb[y][x].B, 0 }; A_stb[y][x] = (PIXEL){ A_stb[w][x].R, A_stb[w][x].G, A_stb[w][x].B, 0 }; A_stb[w][x] = p; }// 色要素順の変更( B ⇔ R )及び 左上原点系 から 左下原点系への変換 } char *Caution[]={ "This template may not work correctly because it is based on \n", "the assumption that sizeof(PIXEL)==sizeof(int). To port to \n", "an environment where int does not use 32 bits, use int32_t \n", "by introducing #include <stdint.h>.\n\n", NULL}; void PX_Warning(){for(int e=0;Caution[e];fprintf(stderr,"%s",Caution[e++]));} int PIXEL_warning=(sizeof(PIXEL)!=sizeof(int));//if(PXCEL_warning)PX_Warning(); int main(int argc, char *argv[]){ char *ifname=NULL, *ofname=NULL; argumentAnalysis( argc, argv, &ifname, &ofname ); if( ifname==NULL ) ifname = "sample.bmp"; if( ofname==NULL ) ofname = "result.bmp"; // 入力画像の読み込み int lx, ly, bpp; PIXEL *stb = (PIXEL *)stbi_load(ifname, &lx, &ly, &bpp, 4); if(!stb){fprintf(stderr, "Couldn't load the image '%s'.\n", ifname); return 1;} else PIXEL_Fitting( stb, lx, ly ); int **s=NULL, **d=NULL, **w=NULL, **t=NULL, **g=NULL, **b=NULL, st=1; const char *err_msg="Memory allocation failure"; if( !(s=(int**)alloc2Darray(sizeof(int),lx,ly)) ){ perror(err_msg); goto GOOD_BY; } if( !(d=(int**)alloc2Darray(sizeof(int),lx,ly)) ){ perror(err_msg); goto GOOD_BY; } if( !(w=(int**)alloc2Darray(sizeof(int),lx,ly)) ){ perror(err_msg); goto GOOD_BY; } if( !(t=(int**)alloc2Darray(sizeof(int),lx,ly)) ){ perror(err_msg); goto GOOD_BY; } if( !(g=(int**)alloc2Darray(sizeof(int),lx,ly)) ){ perror(err_msg); goto GOOD_BY; } if( !(b=(int**)alloc2Darray(sizeof(int),lx,ly)) ){ perror(err_msg); goto GOOD_BY; } /////////////////////////////////////////////////////////////////////////////////// // Calling the image processing program PIXEL **ps=(PIXEL**)s, **pt=(PIXEL**)t; for(int gray, e=0, y=0;y<ly;y++)for(int x=0; x<lx; x++, e++){ ps[y][x] = stb[e]; gray = (306*ps[y][x].R + 601*ps[y][x].G + 117*ps[y][x].B)>>10; pt[y][x] = (PIXEL){ gray, gray, gray, 0 }; g[y][x] = gray; b[y][x] = gray>127 ? 255 : 0; } if( !(st = image_processing_program( s, d, w, t, g, b, lx, ly )) ){ if( !save_without_conversion ){// int → PIXEL 変換 if(auto_adjust) int2PIXEL_adjust( d, d, lx, ly ); else clamp_and_convert(d,(PIXEL**)d,lx,ly); }// BMP保存 PIXEL画像 saveArray2bmp32( ofname, d, lx, ly ); // BMP保存 } //////////////////////////////////////////////////////////////////////////// // メモリ解放 GOOD_BY:; free(b); free(g); free(t); free(w); free(d); free(s); free(stb); return st; }まず、このテンプレートを解読し詳細を解説してください。続けて、image_processing_program()に アルゴリズムを実装してください。 ・二値画像のラベリング処理（4連結、8連結：デフォルト）を実装してもらえますか？ 関数引数の２次元配列 b に格納されている二値化画像が使用できると思います。 ・ラベリング結果は、保存オプションを指定した時にバイナリーデータとして保存してください。 ・出力ファイルは、ラベルデータを疑似カラー変換した画像とします。 関数内で疑似カラー画像を作成した時は、main関数で無変換で保存するように、 save_without_conversionを有効（０以外）にしてください。 ・コンソールに見つかった対象物数（ラベル数）と、処理時間も表示する様にしてください。 ・ユニオン・ファインドを使うなら、回帰呼び出しは使用しない、経路圧縮は必須で実装してください。