以下の機能は、テンプレート中には記述されていません。
１，１次元データ配列を２次以上の高次配列として扱う定義群
// 2次元配列用 （この定義のみ、既に存在している）
#define v2D(type, lx, name) ((type(*)[(lx)])name)
// 3次元配列用
#define v3D(type, lx, ly, name) ((type(*)[(lx)][(ly)])name)
// 4次元配列用
#define v4D(type, lx, ly, lz, name) ((type(*)[(lx)][(ly)][(lz)])name)
// 5次元配列用
#define v5D(type, lx, ly, lz, lt, name) ((type(*)[(lx)][(ly)][(lz)][(lt)])name)
// 6次元配列用
#define v6D(type, lx, ly, lz, lt, lq, name) ((type(*)[(lx)][(ly)][(lz)][(lt)][(lq)])name)
使用例
typedef struct COMPLEX { double r, i; } COMPLEX;
#define v4D(type, lx, ly, lz, name) ((type(*)[(lx)][(ly)][(lz)])name)
#define v4Dcpx(lx, ly, lz, name) v4D(COMPLEX, lx, ly, lz, name)
void func_4D(COMPLEX *s, int lx, int ly, int lz, int lt) {
#undef A_s // s は、COMPLEX型の 1次元配列
#define A_s (v4Dcpx(lx, ly, lz, s)) // A_s は、lx列ly行lz層lt時間を持つ4次元配列
// 例: A_s[t][z][y][x].r で実数成分にアクセス
A_s[t][z][y][x].i = 0.5; // 配列サイズが、lt, lz, ly, lx の複素データ 虚数成分にアクセス
}
補足用
void func_4Dp(COMPLEX ****s, int lx, int ly, int lz, int lt) {
:
// 例: s[t][z][y][x].r で実数成分にアクセス
s[t][z][y][x].i = 0.5; // 配列サイズが、lt, lz, ly, lx の複素データ 虚数成分にアクセス
}
二つの関数の違いが判らず不思議に思う方は、生成AI に尋ねると、 詳しく説明してくれます。
２，抽出系フィルター
////////////////////////////////////////////////////////////////////////////////
// 最小値、中央値、最大値 を抽出する 平均値はおまけ（利用手順の共通化）
//
enum { MIN_EXTRACTION, MED_EXTRACTION, MAX_EXTRACTION, AVE_EXTRACTION }EXTRACTION_FILTER_TYPE;
int compare_int_table(const void *a, const void *b){ return (*(int*)a) - (*(int*)b); }
void ExtractionFilter(int **s, int **d, int lx, int ly, int Kx, int Ky, int type) {
int here=0, size=Kx*Ky, ku = Kx / 2, kv = Ky / 2;// カーネル中心のオフセット
switch(type){
case AVE_EXTRACTION:
#pragma omp parallel for
for(int y=0;y<ly;y++)for(int x=0;x<lx;x++){ int sum = 0;
for(int j=0;j<Ky;j++)for(int i=0;i<Kx;i++){
int cx = x +i -kv;
int cy = y +j -ku;
sum += (cx>=0 && cx<lx && cy>=0 && cy<ly) ? s[cy][cx] : 0;
} d[y][x] = sum/size; // ゼロパディング：境界外は0とみなす
} return;
case MIN_EXTRACTION: here=0; break;
case MED_EXTRACTION: here=(size)>>1; break;
case MAX_EXTRACTION: here=size -1; break;
default: fprintf(stderr,"An unknown operation was specified!\n");
return;
}
#pragma omp parallel for
for(int y=0;y<ly;y++){ int *table=(int*)malloc(sizeof(int)*(size+1));
if(table){
for(int x=0;x<lx;x++){
for(int e=0, j=0;j<Ky;j++)for(int i=0;i<Kx;i++){
int cy = y +j -ku;
int cx = x +i -kv; // ゼロパディング：境界外は0とみなす
table[e++] = (cx>=0 && cx<lx && cy>=0 && cy<ly) ? s[cy][cx] : 0;
} qsort(table, size, sizeof(int), compare_int_table);
d[y][x] = table[here];
} free(table);
}else perror("Memory allocation failure!");
}// s != d は必須条件。ソースデータは、処理が完了するまで変更されてはいけない
}
#define Minimum_filter(s,d,lx,ly,Kx,Ky) ExtractionFilter(s,d,lx,ly,Kx,Ky,MIN_EXTRACTION)
#define Median_filter(s,d,lx,ly,Kx,Ky) ExtractionFilter(s,d,lx,ly,Kx,Ky,MED_EXTRACTION)
#define Maximum_filter(s,d,lx,ly,Kx,Ky) ExtractionFilter(s,d,lx,ly,Kx,Ky,MAX_EXTRACTION)
#define Smoothing_filter(s,d,lx,ly,Kx,Ky) ExtractionFilter(s,d,lx,ly,Kx,Ky,AVE_EXTRACTION)
USER
次の画像処理、データ処理用のテンプレートを評価してください。 少し長めですか、各関数の目的と使用方法を例を加え詳しく説明してください。/* Enhanced C99 Template 2025/06/30 version gcc -O3 -std=c99 -mavx -W -Wall -Wextra -o Template Template.c -fopenmp -lm 環境依存のコードが生成されますが、速度重視の場合は検討してください。 gcc -O3 -std=c99 -march=native -mavx -W -Wall -Wextra -o Template Template.c -fopenmp 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. // // // #define TEMPLATES_USE_AS_INCLUDE_FILE or gcc -DTEMPLATES_USE_AS_INCLUDE_FILE #include "Template.c" */ #ifndef IMAGE_PROCESSING_TEMPLATES_C99 #define IMAGE_PROCESSING_TEMPLATES_C99 //////// https://github.com/nothings/stb/blob/master/stb_image.h #define STB_IMAGE_IMPLEMENTATION #include "stb_image.h" #include <stdio.h> #include <stdlib.h> #include <stdint.h> #include <ctype.h> #include <stdint.h> #include <omp.h> // gcc -fopenmp -O3 ( -O3: Strongly recommendation ) #define HERE fprintf(stderr,"% 6d: %s\n",__LINE__, __func__); #define IF if // ソースコードの行間での位置合わせ用 #define EI else if // 縦に揃えて可読性を向上させるため。 #define EE else // 人間の可読性向上用エイリアス。 extern int verbose; ////////////////////////////////////////////////////////////////////////////////////////// // Timer section #include <time.h> #ifdef _WIN32 // If you'll implement DrawText or DrawTextW yourself. #include <windows.h> // You need, #undef DrawText or #undef DrawTextW 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; // print_elapsed_time(Elapsed_time); print_elapsed_time(E_N_D_time-START_time); #define print_elapsed_time(microsec) fprint_elapsed_time(stdout, (microsec)) void fprint_elapsed_time(FILE *fp, double microsec) { IF(microsec < 1000.0) fprintf(fp,"Elapsed time: %.3f usec\n", microsec); EI(microsec < 1e6 ) fprintf(fp,"Elapsed time: %.3f msec\n", microsec / 1000.0); EE fprintf(fp,"Elapsed time: %.3f sec\n", microsec / 1e6); } void fprintf_elapsed_time(FILE *fp, double microsec, int verbose) { unsigned long long total_us = (unsigned long long)(microsec + 0.5); // 丸め unsigned long long us = total_us % 1000; unsigned long long total_ms = total_us / 1000; unsigned long long ms = total_ms % 1000; unsigned long long total_sec = total_ms / 1000; unsigned long long sec = total_sec % 60; unsigned long long total_min = total_sec / 60; unsigned long long min = total_min % 60; unsigned long long hour = total_min / 60; static const char *time_format[]={ "Elapsed time: ", "%2lluh %2llum %2llus %3llums %3lluus\n", "%lluh %llum %llus %llums %lluus\n" }; fprintf(fp,time_format[0]); if( verbose ){ // 全項目を常に表示（0でも） // fprintf(fp,format[2], hour, min, sec, ms, us); fprintf(fp,time_format[1], hour, min, sec, ms, us); } else { // 必要な項目のみ表示（上位から下位へ） int printed = 0; if (hour > 0 ) { fprintf(fp,"%lluh ", hour); printed = 1; } if ( min > 0 || printed) { fprintf(fp,"%llum ", min); printed = 1; } if ( sec > 0 || printed) { fprintf(fp,"%llus ", sec); printed = 1; } if ( ms > 0 || printed) { fprintf(fp,"%llums ", ms); printed = 1; } printf("%lluus\n", us); } } //////////////////////////////////////////////////////////////////////////////// // Basic structure definition // #pragma pack(push,1) typedef struct PIXEL { uint8_t B, G, R, A; } PIXEL ; typedef struct STBPX { uint8_t 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 DDOUBLE long double #define INF_ZERO06 (1.0e-6) #define INF_ZERO10 (1.0e-10) #define INF_ZERO (1.0e-12) #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) typedef struct HashDB{ char *key; union { unsigned long long here; long long value; void *p_data; }; } HashDB; void **alloc2Darray(int, int, int); int saveDataAsBMP32( char *, void*, int, int ); #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++) #define mtxMul(m,p) (POINTf){\ ((m[0][0]*(p.x)+m[0][1]*(p.y)+m[0][2])/(m[2][0]*(p.x)+m[2][1]*(p.y)+m[2][2])),\ ((m[1][0]*(p.x)+m[1][1]*(p.y)+m[1][2])/(m[2][0]*(p.x)+m[2][1]*(p.y)+m[2][2])) } typedef struct POINTi { int x, y; } POINTi; typedef struct POINTf { float x, y; } POINTf; typedef struct POINTd { double x, y; } POINTd; #include <math.h> #undef PI #define PI 3.1415926535897932384626433832795L // 円周率(π) #define M_PI PI typedef struct COMPLEX { double r, i; } COMPLEX; #define ConjCOMPLEX(a) (COMPLEX){ (a).r, -((a).i) } #define RotorCOMPLEX(a) (COMPLEX){ cos(a), sin(a) } #define AddCOMPLEX(a,b) (COMPLEX){ (a).r + (b).r, (a).i + (b).i } #define SubCOMPLEX(a,b) (COMPLEX){ (a).r - (b).r, (a).i - (b).i } #define MulCOMPLEX(a,b) (COMPLEX){ ((a).r)*((b).r) - ((a).i)*((b).i), \ ((a).r)*((b).i) + ((a).i)*((b).r) } #define ImeDivCOMPLEX(a,b) (COMPLEX){ (a).r / (b), (a).i / (b) } #define ImeMulCOMPLEX(a,b) (COMPLEX){ (a).r * (b), (a).i * (b) } #define ImmediateMulCOMPLEX(b,a) (COMPLEX){ (a).r * (b), (a).i * (b) } #define AppendCOMPLEX(a,b) { (a).r += (b).r, (a).i += (b).i ; } enum { DFT, IDFT, FFT=DFT, IFFT=IDFT, CZT=FFT, ICZT=IFFT }; #define fft_Multiline_with_RF fft_Multiline_with_RotationFactors static inline void fft_Multiline_with_RotationFactors( const COMPLEX *, COMPLEX **, int, int, int); void interchange2D( COMPLEX **d, int lx, int ly ){ COMPLEX t; for(int hly=ly/2, hlx=lx/2, y=0; y<hly; y++ )for(int x=0; x<hlx; x++ ){ t = d[y][x] ; d[y][x] = d[y+hly][x+hlx]; d[y+hly][x+hlx] = t; t = d[y][x+hlx]; d[y][x+hlx] = d[y+hly][x]; d[y+hly][x] = t; } } #define INF_ZERO (1.0e-12) void PowerSpectrum( char *fname, COMPLEX **u, PIXEL **a, int lx, int ly ) { double temp, min, max; min=max=u[0][0].r*u[0][0].r + u[0][0].i*u[0][0].i; for(int y=0; y<ly; y++ )for(int x=0; x<lx; x++ ) { if( max<(temp=u[y][x].r*u[y][x].r + u[y][x].i*u[y][x].i) ) max = temp; else if( min>temp ) min = temp; } min = sqrt(min); max = sqrt(max); double mag = 255.0 / (max>min && max-min>INF_ZERO ? log10(max - min) : 1.0); for(int y=0; y<ly; y++ )for(int x=0; x<lx; x++ ) { int gray = mag*log10(sqrt(u[y][x].r*u[y][x].r + u[y][x].i*u[y][x].i)-min); if(255<gray) gray=255; else if( gray<0 ) gray=0; a[y][x] = (PIXEL){gray,gray,gray,0}; } if(fname)saveArray2bmp32( fname, a, lx, ly ); } //////////////////////////////////////////////////////////////////////////////// //// checks if the argument is a power of 2 ( Returns 1 if it is a power of 2 ) int ft_size_check( int n ){ if( n<4 ) return 0; for(int e=n;e>1;e>>=1 ) if( e&1 ) return 0; return 1; } //////////////////////////////////////////////////////////////////////////////// //// A function that returns the smallest power of 2 that is greater than or //// equal to twice its argument int get_czt_size( int n ){ if( ft_size_check(n) ) return n; int a=1; for( a=1; a<(n<<1); a<<=1 ); return a; } //////////////////////////////////////////////////////////////////////////////// //// A function that is slow but can be used to obtain accurate calculations. void dft_TwoDimensional(COMPLEX **d, int lx, int ly, int flag){int merr=0, elements=lx; if (verbose) printf("execution: %s\t%d\n", __func__, flag); COMPLEX *rotor = (COMPLEX *)calloc(lx > ly ? lx : ly, sizeof(COMPLEX)); if(!rotor){ fprintf(stderr, "Memory allocation failed\n"); return; } else{ elements = lx; // 回転因子の計算 double unit = (flag == IDFT) ? 2 * M_PI / elements : -2 * M_PI / elements; for(int e=0; e<elements; e++) rotor[e] = RotorCOMPLEX(unit * e); } // 行ごとのDFT #pragma omp parallel for for(int y=0; y<ly; y++){ COMPLEX *tmp = malloc(elements * sizeof(COMPLEX)); if( !tmp ){ merr=1; } else { // 各列に対してDFT計算 for(int u=0; u<elements; u++){ COMPLEX t = {0, 0}; for (int v = 0; v < elements; v++) AppendCOMPLEX(t, MulCOMPLEX(d[y][v], rotor[(u * v) % elements])); tmp[u] = t; } if(flag==IDFT){ double reciprocal = 1.0 / elements; //Tips to avoid division for(int e=0;e<elements;e++) d[y][e] = ImeMulCOMPLEX(tmp[e], reciprocal); }EE for(int e=0;e<elements;e++) d[y][e] = tmp[e]; } free(tmp); } // 列ごとのDFT if( lx != ly ){ elements = ly; // 回転因子の再計算 double unit = (flag == IDFT) ? 2 * M_PI / elements : -2 * M_PI / elements; for(int e=0; e<elements; e++) rotor[e] = RotorCOMPLEX(unit * e); } #pragma omp parallel for for(int x=0; x<lx; x++){ COMPLEX *tmp = malloc(elements * sizeof(COMPLEX)); if( !tmp ){ merr=1; } else { for(int u=0; u<elements; u++){ COMPLEX t = {0, 0}; for (int v=0; v<elements; v++) AppendCOMPLEX(t, MulCOMPLEX(d[v][x], rotor[(u * v) % elements])); tmp[u] = t; } if(flag==IDFT){ double reciprocal = 1.0 / elements; //Tips to avoid division for(int e=0;e<elements;e++) d[e][x] = ImeMulCOMPLEX(tmp[e], reciprocal); }EE for(int e=0;e<elements;e++) d[e][x] = tmp[e]; } free(tmp); } if (merr)fprintf(stderr,"Processing interrupted due to memory allocation failure.\n"); free(rotor); } //////////////////////////////////////////////////////////////////////////////// // Fourier Transform CONTROL BLOCK // // samples : number of actual data (Fixed input and output) // conforming_fft: Number of data used in CZT conversion (power of 2) // W[2] : For FFT/IFFT, Twiddle factor table // w : The number of weight data elements is (samples) // v : The number of impulse response data elements(conforming_fft) // // Specialized in image processing, output resizing function has been discarded // N + M -1, (N + M -1)< 2**n. n is .... //////////////////////////////////////////////////////////////////////////////// // ROTATION_FACTOR typedef struct ROTATION_FACTOR { int elements, flag; COMPLEX *table; } ROTATION_FACTOR, ROTATION_TABLE; typedef struct FT_CONTROL { int samples; // 標本点の数 int conforming_fft; // 2の整数乗値：flag_fft?samples:2*get_czt_size(samples) ROTATION_FACTOR W[2];// 回転因子用テーブル COMPLEX *w; // 重みデータ 要素数は(samples) COMPLEX *v; // インパルス応答データ 要素数は(conforming_fft) } FT_CONTROL; //////////////////////////////////////////////////////////////////////////////// //// A function that generates a structure for Chirp Z-transform for a given //// integer value and returns its pointer. //// You MUST call "free()" with this pointer when it is no longer needed. //// FT_CONTROL *GenerateFT_CONTROL( int n ){ int samples = n; int conforming_fft = ft_size_check(n) ? n : get_czt_size(n), nx = conforming_fft; FT_CONTROL *a = (FT_CONTROL*)malloc( sizeof(FT_CONTROL) + sizeof(COMPLEX)*( 2*conforming_fft ) + sizeof(COMPLEX)*( samples + (conforming_fft +1)) ); if(!a){ return NULL;} // +1 for dummy *a = (FT_CONTROL){ // 標本点の数, ２の冪：２の整数乗値 .samples = samples, .conforming_fft = conforming_fft, .W = { { nx, FFT, (COMPLEX*)&a[1] }, // 順変換FFT用回転因子値テーブル { nx, IFFT, &((COMPLEX*)&a[1])[nx] }, // 逆変換FFT用回転因子値テーブル }, .w = &((COMPLEX*)&a[1])[2*nx], // 重みデータ 要素数は(samples) .v = &((COMPLEX*)&a[1])[n+2*nx], // インパルス応答データ 要素数は(conforming_fft) }; COMPLEX *Forward = a->W[0].table, *Reverse = a->W[1].table; // 0:FFT, 1:IFFT double unit = -2*PI/nx; for(int e=0;e<nx;e++) Reverse[e]=ConjCOMPLEX( Forward[e]=RotorCOMPLEX(unit*e) ); COMPLEX *w=a->w, *v=a->v; double theta = PI/n; for(int i=n;i<conforming_fft;i++) v[i]=(COMPLEX){ 0, 0 }; for(int tail=conforming_fft, n2=n<<1, here=0, i=0;i<n;i++ ){ double s, c; sincos( theta*here, &s,&c); // 1.0*cos(), 1.0*sin() amplitude and phase here = ( here +1 + (i<<1) ) % n2 ; // 角度情報 データー数依存なので 1例 factor.mp4 w[i] = (COMPLEX){ c, -s }; // complex conjugate: 0<= { w[i] <-> v[i] } < n v[i] = (COMPLEX){ c, s }; v[tail-i] = (COMPLEX){ c, s }; // stretching ; v[tail] was keeded as dummy area } fft_Multiline_with_RotationFactors( a->W[0].table, &v, conforming_fft, 1, FFT ); return a; } static inline void fft_Multiline_with_RotationFactors( const COMPLEX *rotor, COMPLEX **c, int elements, int lines, int flag ){ //2の冪(べき)のみ、power of two, ONLY if( !ft_size_check(elements) ){ fprintf(stderr,"power of two, ONLY. %d\n",elements); return; } #pragma omp parallel for for(int e=0; e<lines; e++ ){ for(int i=0, j=1; j<elements-1; j++ ) { //// scrambler, for initial placement for(int k=elements>>1; k>(i^=k); k>>=1 ); if( j<i ){COMPLEX temp=c[e][j]; c[e][j]=c[e][i], c[e][i]=temp; } // exchange } for(int m, mh=1; (m=mh<<1)<=elements; mh=m ) { //// Operations for(int ir=0, i=0; i<elements; i+=m ) { // ir=0: Rotation angle is 0 radian. COMPLEX temp, W=rotor[ir]; // W: Rotation factor for(int k=elements>>2; k>(ir^=k); k>>=1 ); for(int j=i; j<mh+i; j++ ) { int a = j, b = j + mh; temp = SubCOMPLEX(c[e][a],c[e][b]);// For butterfly operation c[e][a] = AddCOMPLEX(c[e][a],c[e][b]);// a ->- (a + b) ->- a c[e][b] = MulCOMPLEX(temp, W); // b ->- (a-b)*W ->- b } } } double reciprocal = 1./elements; if(flag==IFFT)for(int i=0; i<elements; c[e][i]=ImeMulCOMPLEX(c[e][i],reciprocal),i++); } } void FourierTransform( COMPLEX **c, int lx, int ly, int flag ){ int st = 1; FT_CONTROL *CtlX = NULL, *CtlY = NULL; COMPLEX **tx = NULL, **ty = NULL; if(!(CtlX = GenerateFT_CONTROL(lx))) goto GOOD_BY; if(!(CtlY = GenerateFT_CONTROL(ly))) goto GOOD_BY; int fft_lx = CtlX->conforming_fft, fft_ly = CtlY->conforming_fft; if(!( tx = (COMPLEX**)alloc2Darray(sizeof(COMPLEX), fft_lx, ly))) goto GOOD_BY; if(!( ty = (COMPLEX**)alloc2Darray(sizeof(COMPLEX), fft_ly, lx))) goto GOOD_BY; int fftX = ft_size_check(lx), fftY = ft_size_check(ly); COMPLEX *RotorXf = CtlX->W[0].table, *RotorXr = CtlX->W[1].table; COMPLEX *RotorYf = CtlY->W[0].table, *RotorYr = CtlY->W[1].table; if( fftX ){ fft_Multiline_with_RF( CtlX->W[(flag==FFT?0:1)].table, c, lx, ly, flag ); for(int y=0;y<ly;y++)for(int x=0;x<lx;x++) ty[x][y] = c[y][x]; }else{ int samples = CtlX->samples, conforming_fft = CtlX->conforming_fft; COMPLEX *w=CtlX->w, *v=CtlX->v; if(flag==CZT){ #pragma omp parallel for for(int y=0;y<ly;y++)for(int e=0; e<samples; e++) tx[y][e] = MulCOMPLEX( c[y][e], w[e] ); fft_Multiline_with_RF( RotorXf, tx, conforming_fft, ly, FFT ); #pragma omp parallel for for(int y=0;y<ly;y++)for(int e=0; e<conforming_fft; e++) tx[y][e] = MulCOMPLEX(tx[y][e], v[e] ); fft_Multiline_with_RF( RotorXr, tx, conforming_fft, ly, IFFT ); #pragma omp parallel for for(int y=0;y<ly;y++)for(int e=0; e<samples; e++) ty[e][y] = MulCOMPLEX(tx[y][e], w[e] ); }else{ double m = 1.0/samples; // ICZT #pragma omp parallel for for(int y=0;y<ly;y++)for(int e=0; e<samples; e++) tx[y][e] = MulCOMPLEX( c[y][e], ConjCOMPLEX(w[e]) ); fft_Multiline_with_RF( RotorXf, tx, conforming_fft, ly, FFT ); #pragma omp parallel for for(int y=0;y<ly;y++)for(int e=0; e<conforming_fft; e++) tx[y][e] = MulCOMPLEX(tx[y][e], ConjCOMPLEX(v[e]) ); fft_Multiline_with_RF( RotorXr, tx, conforming_fft, ly, IFFT ); #pragma omp parallel for for(int y=0;y<ly;y++)for(int e=0; e<samples; e++) ty[e][y] = ImmediateMulCOMPLEX( m, MulCOMPLEX(tx[y][e], ConjCOMPLEX(w[e]) )); } } if( fftY ){ fft_Multiline_with_RF( CtlY->W[(flag==FFT?0:1)].table, ty, ly, lx, flag ); for(int y=0;y<ly;y++)for(int x=0;x<lx;x++) c[y][x] = ty[x][y]; }else{ int samples = CtlY->samples, conforming_fft = CtlY->conforming_fft; COMPLEX *w=CtlY->w, *v=CtlY->v; if(flag==CZT){ #pragma omp parallel for for(int x=0;x<lx;x++)for(int e=0; e<samples; e++) ty[x][e] = MulCOMPLEX(ty[x][e], w[e] ); fft_Multiline_with_RF( RotorYf, ty, conforming_fft, lx, FFT ); #pragma omp parallel for for(int x=0;x<lx;x++)for(int e=0; e<conforming_fft; e++) ty[x][e] = MulCOMPLEX(ty[x][e], v[e] ); fft_Multiline_with_RF( RotorYr, ty, conforming_fft, lx, IFFT ); #pragma omp parallel for for(int x=0;x<lx;x++)for(int e=0; e<samples; e++) c[e][x] = MulCOMPLEX(ty[x][e], w[e] ); }else{ double m = 1.0/samples; // ICZT #pragma omp parallel for for(int x=0;x<lx;x++)for(int e=0; e<samples; e++) ty[x][e] = MulCOMPLEX(ty[x][e], ConjCOMPLEX(w[e]) ); fft_Multiline_with_RF( RotorYf, ty, conforming_fft, lx, FFT ); #pragma omp parallel for for(int x=0;x<lx;x++)for(int e=0; e<conforming_fft; e++) ty[x][e] = MulCOMPLEX(ty[x][e], ConjCOMPLEX(v[e]) ); fft_Multiline_with_RF( RotorYr, ty, conforming_fft, lx, IFFT ); #pragma omp parallel for for(int x=0;x<lx;x++)for(int e=0; e<samples; e++) c[e][x] = ImmediateMulCOMPLEX( m, MulCOMPLEX(ty[x][e], ConjCOMPLEX(w[e]) ) ); } } st = 0; GOOD_BY:; if(st)fprintf(stderr,"FourierTransform() failed to allocate working memory area.\n"); free(ty); free(tx); free(CtlY); free(CtlX); } //////////////////////////////////////////////////////////////////////////////// // Global variables int verbose=0, nega=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 ){ // alloc2D: 既存のメモリブロックに2Dインデックスを作成 // data: 初期化済みのメモリブロック（未初期化の場合、未定義動作） // s: 既存ブロックの１要素分のサイズ, lx: 幅, 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){ if(s<1 || lx<1 || ly<1) return NULL; long long w = lx, h = ly, next = s*w, space = h*sizeof(void *) + w*h*s + 64; void **d = (void **)calloc(1, space); if(!d){fprintf(stderr,"alloc2Darray failed: %lld bytes\n",space); return NULL;} void *base = (void*)((uintptr_t)((char*)d + h * sizeof(void*) + 63) & ~63LL); for(long long i=0; i<h; i++) d[i] = (void*)((char*)base + next * i); return d; // C99: uintptr_t same as (void *) } 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; if(div==1) min-=127; PIXEL **pd = (PIXEL**)d; for(int value, y=0;y<ly;y++)for(int x=0;x<lx;x++){ value = 255*(s[y][x]-min)/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; } //////////////////////////////////////////////////////////////////////////////// // パブリックドメイン利用の32ビットのハッシュ関数 セキュリティ用の強度は無い事に注意 // // Public domain. FNV-1a based on http://www.isthe.com/chongo/tech/comp/fnv/ // Public domain. MurmurHash2 based on https://github.com/aappleby/smhasher, // modified with expansion for collision resistance // And, following code same too. #include <stdint.h> // FNV-1a uint32_t fnv1a_hash(const void *key, size_t len){ if( !key ) return 0; uint32_t hash = 0x811C9DC5; const uint8_t *data = (const uint8_t *)key; for(size_t i=0; i<len; i++){ hash ^= data[i]; hash *= 0x01000193; } return hash; } // MurmurHash2: 高速かつ高品質な非暗号学的ハッシュ関数 // key: ハッシュ対象のデータ // len: データのバイト長 // seed: ハッシュの初期シード値 uint32_t MurmurHash2(const void *key, int len, uint32_t seed) { const uint32_t m = 0x5bd1e995; // 'm' と 'r' はミキシング定数。 const int r = 24; // 特定の値は衝突を最小化するために選ばれている uint32_t h = seed ^ len; // ハッシュ値をシードで初期化 // データポインタをバイト単位で処理 const unsigned char *data = (const unsigned char *)key; while( len>=4 ){// 4バイト単位でデータを処理 // 4バイトをuint32_tとして取得（リトルエンディアン） uint32_t k = *(uint32_t *)data; k *= m; // ミキシングステップ k ^= k >> r; k *= m; h *= m; // ハッシュ値にミックス h ^= k; data += 4; // 次の4バイトへ len -= 4; } // 残りのバイトを処理（末尾処理） switch( len ){ case 3: h ^= data[2] << 16; // Fallthrough case 2: h ^= data[1] << 8; // Fallthrough case 1: h ^= data[0]; h *= m; } h ^= h >> 13; // 最終ミキシング（アバランシェ効果を強化） h *= m; h ^= h >> 15; return h; } // Create a hash data block for MurmurHash 32bit version (MAX HashDB: 64G Bytes) HashDB* make_HashDBmax(){ return (HashDB*)calloc(0x100000000,sizeof(HashDB)); } // MAX 40億件 #define make_HashDB512M() make_HashDB( 0x20000000 ) // MAX 5.3億件 8G #define make_HashDB256M() make_HashDB( 0x10000000 ) // MAX 2.5億件 4G HashDB* make_HashDB(uint32_t size){ return (HashDB*)calloc(size,sizeof(HashDB)); } // uint32_t HDB512 = (MurmurHash2( const void *key, size_t len, uint32_t seed ))>>3; 広く分布させる使い方 // uint32_t HDB256 = (MurmurHash2( const void *key, size_t len, uint32_t seed ))>>4; 広く分布させる使い方 // 上位ビットを捨てると狭い範囲に集中し易くなる、コンフリクトを起こしやすくなる。 void HomographyTrans( PIXEL **c, int lx, int ly, PIXEL **s, int slx, int sly, double inv[3][3] ){ #pragma omp parallel for for(int y=0;y<ly;y++)for(int x=0;x<lx;x++){ POINTf destination = (POINTf){x,y}; POINTf source = mtxMul( inv, destination ); int cx=source.x +.5, cy=source.y +.5; c[y][x] = (cx<0 || cy<0 || slx<=cx || sly<=cy) ? (PIXEL){255,255,255,0} : s[cy][cx]; } } void HomographyTransMask( PIXEL **c, int lx, int ly, PIXEL **s, int slx, int sly, int **mask, double inv[3][3] ){ #pragma omp parallel for for(int y=0;y<ly;y++)for(int x=0;x<lx;x++){ POINTf destination = (POINTf){x,y}; POINTf source = mtxMul( inv, destination ); int cx=source.x +.5, cy=source.y +.5; if( (cx<0 || cy<0 || slx<=cx || sly<=cy) || mask[cy][cx] )continue; c[y][x] = s[cy][cx]; } } ///////////////////////////////////////////////////////////////////////////////////////////////// // DrawXXXX functions // void DrawLine(PIXEL **canvas, int lx, int ly, int x0, int y0, int x1, int y1, PIXEL pen) { int dx = abs(x1 - x0), sx = (x0 < x1) ? 1 : -1; int dy = abs(y1 - y0), sy = (y0 < y1) ? 1 : -1; int err = dx - dy; while(1){ if(x0>=0 && x0<lx && y0>=0 && y0<ly) canvas[y0][x0] = pen; if (x0 == x1 && y0 == y1) break; int e2 = 2 * err; if (e2 > -dy) { err -= dy; x0 += sx; } if (e2 < dx) { err += dx; y0 += sy; } } } void DrawLineExor(PIXEL **canvas, int lx, int ly, int x0, int y0, int x1, int y1, PIXEL pen) { int dx = abs(x1 - x0), sx = (x0 < x1) ? 1 : -1; int dy = abs(y1 - y0), sy = (y0 < y1) ? 1 : -1; int **c=(int**)canvas, p=*(int*)&pen, err = dx - dy; while(1){ if(x0>=0 && x0<lx && y0>=0 && y0<ly) c[y0][x0] ^= p; if (x0 == x1 && y0 == y1) break; int e2 = 2 * err; if (e2 > -dy) { err -= dy; x0 += sx; } if (e2 < dx) { err += dx; y0 += sy; } } } void DrawPolyline(PIXEL **canvas, int lx, int ly, int num_points, POINTi *points, int close, PIXEL pen) { if( num_points < 2 ) return; // エッジケースのチェック int t = num_points -1; for(int i=0; i<t; i++ ) DrawLine(canvas,lx,ly, points[i].x, points[i].y, points[i+1].x, points[i+1].y,pen); if(close&&num_points>2) DrawLine(canvas,lx,ly, points[t].x, points[t].y, points[0].x, points[0].y, pen); } void DrawCircle(PIXEL **canvas, int lx, int ly, int cx, int cy, int r, PIXEL pen) { for(int x=0, y=r, d=1-r; x<=y; x++) { if(cx + x >= 0 && cx + x < lx && cy + y >= 0 && cy + y < ly) canvas[cy + y][cx + x] = pen; if(cx - x >= 0 && cx - x < lx && cy + y >= 0 && cy + y < ly) canvas[cy + y][cx - x] = pen; if(cx + x >= 0 && cx + x < lx && cy - y >= 0 && cy - y < ly) canvas[cy - y][cx + x] = pen; if(cx - x >= 0 && cx - x < lx && cy - y >= 0 && cy - y < ly) canvas[cy - y][cx - x] = pen; if(cx + y >= 0 && cx + y < lx && cy + x >= 0 && cy + x < ly) canvas[cy + x][cx + y] = pen; if(cx - y >= 0 && cx - y < lx && cy + x >= 0 && cy + x < ly) canvas[cy + x][cx - y] = pen; if(cx + y >= 0 && cx + y < lx && cy - x >= 0 && cy - x < ly) canvas[cy - x][cx + y] = pen; if(cx - y >= 0 && cx - y < lx && cy - x >= 0 && cy - x < ly) canvas[cy - x][cx - y] = pen; if(d<0) d = d + 2 * x + 3; else { d = d + 2 * (x - y) + 5; y--; } } } void DrawArc(PIXEL **canvas, int lx, int ly, int cx, int cy, int r, int start_angle, int end_angle, PIXEL pen) { start_angle = ((start_angle % 360) + 360) % 360; end_angle = (( end_angle % 360) + 360) % 360; int wrap_around = (start_angle > end_angle); for(int x=0, y=r, d=1-r; x<=y; x++) { // 8つの対称点をチェック int points[8][2] = { {cx + x, cy + y}, {cx - x, cy + y}, {cx + x, cy - y}, {cx - x, cy - y}, {cx + y, cy + x}, {cx - y, cy + x}, {cx + y, cy - x}, {cx - y, cy - x} }; for(int i = 0; i < 8; i++) { int px = points[i][0]; int py = points[i][1]; if (px >= 0 && px < lx && py >= 0 && py < ly) { double angle = atan2(py - cy, px - cx) * 180.0 / M_PI; if (angle < 0) angle += 360; if ((!wrap_around && angle >= start_angle && angle <= end_angle) || (wrap_around && (angle >= start_angle || angle <= end_angle))) { canvas[py][px] = pen; } } } if(d<0) d = d + 2 * x + 3; else { d = d + 2 * (x - y) + 5; y--; } } } void DrawEllipse(PIXEL **canvas, int lx, int ly, int x0, int y0, int x1, int y1, PIXEL pen) { int left = (x0 < x1) ? x0 : x1; int right = (x0 > x1) ? x0 : x1; int top = (y0 < y1) ? y0 : y1; int bottom = (y0 > y1) ? y0 : y1; int cx = left + ( right - left) / 2; int cy = top + (bottom - top ) / 2; int a = (right - left) / 2; int b = (bottom - top ) / 2; if (a <= 0 || b <= 0) return; const double step_deg = 0.1; for(double angle = 0; angle < 360; angle += step_deg) { double rad = angle * M_PI / 180.0; int px = cx + round(a * cos(rad)); int py = cy + round(b * sin(rad)); if(px >= 0 && px < lx && py >= 0 && py < ly) canvas[py][px] = pen; } } ////////////////////////////////////////////////////////////////////////////////////////// // Function set for algorithm verification // #define dabs(a) ((a)<0? -(a):(a)) void lu_decomposition_DoublePrecision(double** A, long double** L, long double** U, long long* P, int n) { for(int i=0; i<n; i++){ P[i] = i; for(int j=0; j<n; j++) { L[i][j] = 0.0; U[i][j] = A[i][j]; } } for(int k=0; k<n; k++) { double max = 0.0; int k_prime = k; for(int i=k; i<n; i++) { if(dabs(U[i][k]) > max) { max = dabs(U[i][k]); k_prime = i; } } if(k_prime != k){ long double* temp_row; // 入れ替え long long tmp = P[k]; P[k] = P[k_prime]; P[k_prime] = tmp; temp_row = U[k]; U[k] = U[k_prime]; U[k_prime] = temp_row; temp_row = L[k]; L[k] = L[k_prime]; L[k_prime] = temp_row; } for(int i=k +1; i<n; i++){ L[i][k] = U[i][k] / U[k][k]; for(int j=k; j<n; j++) U[i][j] -= L[i][k] * U[k][j]; } for(int j=0; j<=k; j++) L[k][j] = (k==j) ? 1.0 : L[k][j]; } } void solve_LU_DoublePrecision(long double** L, long double** U, long long* P, long double* b, long double* x, int n) { long double* y = (long double*)malloc(n * sizeof(long double)); long double* pb = (long double*)malloc(n * sizeof(long double)); if( y && pb ){ for(int i=0; i<n; i++) pb[i] = b[P[i]]; // 前進代入（Ly = Pb） for(int i=0; i<n; i++) { y[i] = pb[i]; for(int j=0; j<i; j++) y[i] -= L[i][j] * y[j]; } // 後退代入（Ux = y） for(int i= n -1; i>=0; i--) { x[i] = y[i]; for(int j=i +1; j<n; j++) x[i] -= U[i][j] * x[j]; x[i] /= U[i][i]; } } else { perror("Memory allocation error in solve_LU_DoublePrecision\n"); exit(1); } free(y); free(pb); } void inverse_matrix_DoublePrecision(double** A, double** A_inv, int n) { long double** L = (long double**)malloc(n * sizeof(long double*)); long double** U = (long double**)malloc(n * sizeof(long double*)); long long* P = (long long*)malloc(n * sizeof(long long)); if( L && U && P ){ for(int i=0; i<n; i++) { L[i] = (long double*)calloc(n, sizeof(long double)); U[i] = (long double*)calloc(n, sizeof(long double)); if( !L[i] && !U[i] ){ perror("Memory allocation error L[i],U[i] in inverse_matrix\n"); exit(1); } } lu_decomposition_DoublePrecision(A, L, U, P, n); // 並列化：列ごとに並列に解く #pragma omp parallel for schedule(dynamic) for(int i=0; i<n; i++) { long double* e = (long double*)calloc(n, sizeof(long double)); long double* x = (long double*)calloc(n, sizeof(long double)); if( e && x ){ e[i] = 1.0; solve_LU_DoublePrecision(L, U, P, e, x, n); for(int j=0; j<n; j++) A_inv[j][i] = x[j]; } else { perror("Memory allocation error for solve_LU_DoublePrecision\n"); exit(1); } free(e); free(x); }// if( n>=1000 )printf("\n" ); } else { perror("Memory allocation error in inverse_matrix\n"); exit(1); } for(int i=0; i<n; i++){ free(L[i]); free(U[i]); } free(L); free(U); free(P); } void lu_decomposition(double** A, double** L, double** U, int* P, int n) { for(int i=0; i<n; i++){ P[i] = i; for(int j=0; j<n; j++) { L[i][j] = 0.0; U[i][j] = A[i][j]; } } for(int k=0; k<n; k++) { double max = 0.0; int k_prime = k; for(int i=k; i<n; i++) { if(fabs(U[i][k]) > max) { max = fabs(U[i][k]); k_prime = i; } } if(k_prime != k){ double* temp_row; // 入れ替え int tmp = P[k]; P[k] = P[k_prime]; P[k_prime] = tmp; temp_row = U[k]; U[k] = U[k_prime]; U[k_prime] = temp_row; temp_row = L[k]; L[k] = L[k_prime]; L[k_prime] = temp_row; } for(int i=k +1; i<n; i++){ L[i][k] = U[i][k] / U[k][k]; for(int j=k; j<n; j++) U[i][j] -= L[i][k] * U[k][j]; } for(int j=0; j<=k; j++) L[k][j] = (k==j) ? 1.0 : L[k][j]; } } void solve_LU(double** L, double** U, int* P, double* b, double* x, int n) { double* y = (double*)malloc(n * sizeof(double)); double* pb = (double*)malloc(n * sizeof(double)); if( y && pb ){ for(int i=0; i<n; i++) pb[i] = b[P[i]]; // 前進代入（Ly = Pb） for(int i=0; i<n; i++) { y[i] = pb[i]; for(int j=0; j<i; j++) y[i] -= L[i][j] * y[j]; } // 後退代入（Ux = y） for(int i= n -1; i>=0; i--) { x[i] = y[i]; for(int j=i +1; j<n; j++) x[i] -= U[i][j] * x[j]; x[i] /= U[i][i]; } } else { perror("Memory allocation error in solve_LU\n"); exit(1); } free(y); free(pb); } void inverse_matrix(double** A, double** A_inv, int n) { double** L = (double**)malloc(n * sizeof(double*)); double** U = (double**)malloc(n * sizeof(double*)); int* P = (int*)malloc(n * sizeof(int)); if( L && U && P ){ for(int i=0; i<n; i++) { L[i] = (double*)calloc(n, sizeof(double)); U[i] = (double*)calloc(n, sizeof(double)); if( !L[i] && !U[i] ){ perror("Memory allocation error L[i],U[i] in inverse_matrix\n"); exit(1); } } lu_decomposition(A, L, U, P, n); // 並列化：列ごとに並列に解く #pragma omp parallel for for(int i=0; i<n; i++) { double* e = (double*)calloc(n, sizeof(double)); double* x = (double*)calloc(n, sizeof(double)); if( e && x ){ e[i] = 1.0; solve_LU(L, U, P, e, x, n); for(int j=0; j<n; j++) A_inv[j][i] = x[j]; } else { perror("Memory allocation error for solve_LU\n"); exit(1); } free(e); free(x); } } else { perror("Memory allocation error in inverse_matrix\n"); exit(1); } for(int i=0; i<n; i++){ free(L[i]); free(U[i]); } free(L); free(U); free(P); } // 行列積を計算 void multiplication_matrix(double **A, double **B, double **C, int n) { #pragma omp parallel for for(int i=0; i<n; i++){ for(int j=0; j<n; j++){ double sum = 0.0; for(int k=0; k<n; k++) sum += A[i][k] * B[k][j]; C[i][j] = sum; } } } void dump_array(FILE *fp, char *fmt, char *msg, double ** A, int N) { if(msg) fprintf(fp, "\n%s\n",msg); else fprintf(fp, "\n"); for(int y=0; y<N; y++) { fprintf(fp, "% 6d: \t", y); for(int x=0; x<N; x++) fprintf(fp, fmt, A[y][x] + 0.0000001); fprintf(fp, "\n"); } fprintf(fp, "\n"); } void HomographyCore( POINTf sp[4], POINTf dp[4], double hkl[3][3] ){ POINTf SA = sp[0], SB = sp[1], SC = sp[2], SD = sp[3]; POINTf DA = dp[0], DB = dp[1], DC = dp[2], DD = dp[3]; double tranp[8] = { DA.x, DA.y, DB.x, DB.y, DC.x, DC.y, DD.x, DD.y }; double abcdefgh1[3*3]={0,0,0,0,0,0,0,0,1}; //Coefficient of 3x3 matrix double hg[8][8] = { { SA.x, SA.y, 1, 0, 0, 0, -SA.x*DA.x, -SA.y*DA.x }, { 0, 0, 0, SA.x, SA.y, 1, -SA.x*DA.y, -SA.y*DA.y }, { SB.x, SB.y, 1, 0, 0, 0, -SB.x*DB.x, -SB.y*DB.x }, { 0, 0, 0, SB.x, SB.y, 1, -SB.x*DB.y, -SB.y*DB.y }, { SC.x, SC.y, 1, 0, 0, 0, -SC.x*DC.x, -SC.y*DC.x }, { 0, 0, 0, SC.x, SC.y, 1, -SC.x*DC.y, -SC.y*DC.y }, { SD.x, SD.y, 1, 0, 0, 0, -SD.x*DD.x, -SD.y*DD.x }, { 0, 0, 0, SD.x, SD.y, 1, -SD.x*DD.y, -SD.y*DD.y }, },invhg[8][8], *phg[8], *pinvhg[8]; // inverseHomography ////////////////////////////////////////////////////////////////////// for(int e=0; e<8; e++){ phg[e]=hg[e]; pinvhg[e]=invhg[e]; } inverse_matrix( phg, pinvhg, 8 ); // Calculate the inverse matrix for(int y=0;y<8;y++)for(int x=0;x<8;x++) abcdefgh1[y] += invhg[y][x]*tranp[x]; for(int z=0,y=0;y<3;y++)for(int x=0;x<3;x++) hkl[y][x] = abcdefgh1[z++]; } void Homography( PIXEL **s, int slx, int sly, int **mask, POINTf sp[4], PIXEL **d, int lx, int ly, POINTf dp[4] ){ ////////////////////////////////////////////////////////////////////// double homographyMatrix[3][3], transMatrix[3][3]; double *p_homographyMatrix[3], *p_transMatrix[3]; for(int e=0; e<3; e++){ p_homographyMatrix[e]=homographyMatrix[e]; p_transMatrix[e] =transMatrix[e]; } HomographyCore( sp, dp, homographyMatrix ); ////////////////////////////////////////////////////////////////////// // Calculate the inverse matrix for image transportation inverse_matrix( p_homographyMatrix, p_transMatrix, 3 ); if(!mask) HomographyTrans( d, lx, ly, s, slx, sly, transMatrix ); else HomographyTransMask( d, lx, ly, s, slx, sly, mask, transMatrix ); } //////////////////////////////////////////////////////////////////////////////// // ユニオン・ファインド構造体 typedef struct { int *parent; // 親ノード int *rank; // ランク（木の高さの近似） int size; // 配列サイズ } UnionFind; // ユニオン・ファインドの初期化 UnionFind *uf_init(int size) { UnionFind *uf = (UnionFind *)malloc(sizeof(UnionFind)); if(!uf) return NULL; uf->parent = (int *)malloc(size * sizeof(int)); uf->rank = (int *)malloc(size * sizeof(int)); if (!uf->parent || !uf->rank) { free(uf->parent); free(uf->rank); free(uf); return NULL; } uf->size = size; for(int i=0; i<size; i++) { uf->parent[i] = i; // 各ノードは自分自身を親に uf->rank[i] = 0; } return uf; } // ルート検索（経路圧縮） int uf_find(UnionFind *uf, int x) { if(uf->parent[x] != x) uf->parent[x] = uf_find(uf,uf->parent[x]);// 経路圧縮 return uf->parent[x]; } // 結合（ランク付き） void uf_union(UnionFind *uf, int x, int y) { int root_x = uf_find(uf, x); int root_y = uf_find(uf, y); if(root_x == root_y) return; if(uf->rank[root_x] < uf->rank[root_y]) uf->parent[root_x] = root_y; EI(uf->rank[root_x] > uf->rank[root_y]) uf->parent[root_y] = root_x; EE{ uf->parent[root_y] = root_x; uf->rank[root_x]++; } } // ユニオン・ファインドの解放 void uf_free(UnionFind *uf) { if(uf){ free(uf->parent); free(uf->rank); free(uf); } } // ラベル数を取得 int uf_count_labels(UnionFind *uf, int *labels, int size) { int *counted = (int *)calloc(uf->size, sizeof(int)); int label_count = 0; for(int i = 0; i < size; i++) { if (labels[i] > 0) { int root = uf_find(uf, labels[i]); if (!counted[root]) { counted[root] = 1; label_count++; } } } free(counted); return label_count; } // バイナリファイルにラベルデータを保存 int save_label_data(char *filename, int **labels, int lx, int ly) { FILE *fp = fopen(filename, "wb"); if (!fp) { fprintf(stderr, "Cannot create binary file '%s'.\n", filename); return -1; } for(int y = 0; y < ly; y++) { if (fwrite(labels[y], sizeof(int), lx, fp) != (size_t)lx) { fprintf(stderr, "Write error to '%s'.\n", filename); fclose(fp); return -1; } } fclose(fp); return 0; } int labeling(int **b, int **labels, int lx, int ly, int connectivity) { if( !b || !labels || lx<1 || ly<1 ) return -1; if( connectivity!=4 && connectivity!=8 ) return -1; PROCESS_ALL_PIXELS labels[y][x]=0; // ラベル配列（初期値0） UnionFind *uf = uf_init(lx * ly); // ユニオン・ファインドの初期化 int *label_map = (int *)calloc(lx*ly+1,sizeof(int)); // ラベルを連続番号に再割り当て用 if( !uf || !label_map){ if(!uf) fprintf(stderr, "Union-Find initialization failed.\n"); EE fprintf(stderr, "Memory allocation failure for re-labels.\n");; free(label_map); uf_free(uf); return -1; } // 1パス目：Union-Find操作（uf_union, uf_find）の競合を避けるため逐次処理。 label:ラベル番号のカウンタ if( connectivity==4 ){ for(int label=1, y=0; y<ly; y++)for(int x=0; x<lx; x++) { if( !b[y][x] ) continue; // 背景ならスキップ、前景（白色ピクセル）でないとき int min_label = 0; if( y>0 && b[y-1][x] && labels[y-1][x]>0 ) min_label = labels[y-1][x]; // 上 if( x>0 && b[y][x-1] && labels[y][x-1]>0 ) // 左 if( min_label==0 || labels[y][x-1]<min_label ) min_label = labels[y][x-1]; if(min_label==0) labels[y][x] = label++; // 新しいラベルを割り当て else { labels[y][x] = min_label; // 近傍のラベルを結合 if( y>0 && b[y-1][x] && labels[y-1][x]>0 && labels[y-1][x] != min_label ) uf_union(uf, labels[y-1][x], min_label); if( x>0 && b[y][x-1] && labels[y][x-1]>0 && labels[y][x-1] != min_label ) uf_union(uf, labels[y][x-1], min_label); } } } else { for(int label=1, y=0; y<ly; y++)for(int x=0; x<lx; x++) { if( !b[y][x] ) continue; // 背景ならスキップ、前景（白色ピクセル）でないとき int min_label = 0; if( y>0 && b[y-1][x] && labels[y-1][x]>0 ) min_label = labels[y-1][x]; // 上 if( x>0 && b[y][x-1] && labels[y][x-1]>0 ) // 左 if( min_label==0 || labels[y][x-1]<min_label ) min_label = labels[y][x-1]; if( y>0 && x>0 && b[y-1][x-1] && labels[y-1][x-1]>0 ) //8: 左上 if( min_label==0 || labels[y-1][x-1]<min_label ) min_label = labels[y-1][x-1]; //8: if( y>0 && x < lx-1 && b[y-1][x+1] && labels[y-1][x+1]>0 ) //8: 右上 if( min_label==0 || labels[y-1][x+1]<min_label ) min_label = labels[y-1][x+1]; //8: if(min_label==0) labels[y][x] = label++; // 新しいラベルを割り当て else { labels[y][x] = min_label; // 近傍のラベルを結合 if( y>0 && b[y-1][x] && labels[y-1][x]>0 && labels[y-1][x] != min_label ) uf_union(uf, labels[y-1][x], min_label); if( x>0 && b[y][x-1] && labels[y][x-1]>0 && labels[y][x-1] != min_label ) uf_union(uf, labels[y][x-1], min_label); if( y>0 && x>0 && b[y-1][x-1] && labels[y-1][x-1]>0 && labels[y-1][x-1] != min_label ) uf_union(uf, labels[y-1][x-1], min_label);//8: if( y>0 && x<lx-1 && b[y-1][x+1] && labels[y-1][x+1]>0 && labels[y-1][x+1] != min_label ) uf_union(uf, labels[y-1][x+1], min_label);//8: } } } // 2パス目：ラベルをルートに統一 #pragma omp parallel for collapse(2) for(int y=0; y<ly; y++) { for(int x=0; x<lx; x++) if(labels[y][x]>0) labels[y][x] = uf_find(uf, labels[y][x]); } // ラベルを連続番号に再割り当て int new_label = 1; PROCESS_ALL_PIXELS { if(labels[y][x] > 0 && label_map[labels[y][x]]==0)label_map[labels[y][x]] = new_label++; labels[y][x] = label_map[labels[y][x]]; } // メモリ解放 free(label_map); uf_free(uf); return new_label - 1; } //////////////////////////////////////////////////////////////////////////////// // 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; #pragma omp parallel for 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; // カーネル中心のオフセット #pragma omp parallel for 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 }; } } #define int2PC int2PseudoColor // ~ : ビット反転(ビット毎に反転)、 ! : 否定(0 なら1、以外 0) void image_set(int v, int **dst, int lx, int ly){ PROCESS_ALL_PIXELS dst[y][x] = v; } void image_pc( int **a, int **dst, int lx, int ly){ PROCESS_ALL_PIXELS dst[y][x] = int2PC(a[y][x]); } void image_copy(int **a, int **dst, int lx, int ly){ PROCESS_ALL_PIXELS dst[y][x] = a[y][x]; } void image_not(int **a, int **dst, int lx, int ly){ PROCESS_ALL_PIXELS dst[y][x] = ~a[y][x]; } void image_and(int **a, int **b, int **dst, int lx, int ly){ PROCESS_ALL_PIXELS dst[y][x] = a[y][x] & b[y][x];} void image_or( int **a, int **b, int **dst, int lx, int ly){ PROCESS_ALL_PIXELS dst[y][x] = a[y][x] | b[y][x];} void image_xor(int **a, int **b, int **dst, int lx, int ly){ PROCESS_ALL_PIXELS dst[y][x] = a[y][x] ^ b[y][x];} void image_add(int **a, int **b, int **dst, int lx, int ly){ PROCESS_ALL_PIXELS dst[y][x] = a[y][x] + b[y][x];} void image_sub(int **a, int **b, int **dst, int lx, int ly){ PROCESS_ALL_PIXELS dst[y][x] = a[y][x] - b[y][x];} void image_mul(int **a, int **b, int **dst, int lx, int ly){ PROCESS_ALL_PIXELS dst[y][x] = a[y][x] * b[y][x];} void image_div(int **a, int **b, int **dst, int lx, int ly){ PROCESS_ALL_PIXELS 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) { PROCESS_ALL_PIXELS dst[y][x] = abs(a[y][x]) + abs(b[y][x]) + abs(c[y][x]); } void clamp_and_convert(int **src, PIXEL **dst, int lx, int ly) { PROCESS_ALL_PIXELS{int v = src[y][x]<0 ? 0:(src[y][x]>255 ? 255:src[y][x]); dst[y][x] = (PIXEL){v, v, v, 0}; } } // templates use as include file. SET => TEMPLATES_USE_AS_INCLUDE_FILE #ifndef TEMPLATES_USE_AS_INCLUDE_FILE ////////////////////////////////////////////////////////////////////////////////////////// // Enhanced version // Algorithm implementation area, start here // // 関数内で画像データを作成した時 整数=>PIXEL（BMP形式）変換を行わない int save_without_conversion = 0; // コンボリューション、ガウス関数等 int kernel = 5; float sigma = 1.0; // コマンドラインオプション ラベリング処理用 // デフォルトは8連結 バイナリ保存フラグ バイナリファイル名 //static int connectivity = 8, save_binary = 0; static char *binary_file = NULL; int image_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 に二値化したデータ return 0; } // Algorithm implementation area. ends here ////////////////////////////////////////////////////////////////////////////////////////// #include <string.h> void argumentAnalysis(int argc, char *argv[], char **iname, char **oname){ int status=1; 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 'n': if(argv[i][2]) goto GOOD_BY; nega=1; break; //case 'c': if(argv[i][3]) goto GOOD_BY; // ラベリング処理用 // if(argv[i][2]=='4') connectivity = 4; // EI(argv[i][2]=='8') connectivity = 8; // EE goto GOOD_BY; // 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 'b': if( argv[i][2] ) binary_file = &argv[i][2]; // EI( ++i>=argc ) goto GOOD_BY; // EE binary_file = argv[i]; //ラベリング処理用 // save_binary=1; break; case 'h': status=0; // fallthrough default : goto GOOD_BY; // In gcc, you can use the magic comment "// fallthrough" to } // suppress the warning. without "goto GOOD_BY;" or "break;" } // warning: this statement may fall through [-Wimplicit-fallthrough=] 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] [-n] [-t<num>] [-h]"; // "[-k|-kernel <num>] [-s|-sigma <num>] [-c{4|8}] [-b filename]" 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-n: Invert the brightness of the input image.\n" ); //fprintf( stderr, "\t\t-c: -c{4|8}. Specify the connection method (4 or 8).\n" ); //fprintf( stderr, "\t\t-k: -k|-kernel <num>. For specifying kernels such as convolution.\n" ); //fprintf( stderr, "\t\t-s: -s|-sigma <num>. For specifying standard deviation etc.\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( status ); } void PIXEL_Fitting( PIXEL *s, int lx, int ly ) { #undef A_stb #define A_stb (v2Dp(lx,s)) for(int ey=(ly +1)/2, ty=ly-1, y=0; y<ey; y++ )for(int 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 )及び 左上原点系 から 左下原点系への変換 } 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'.%s\n",ifname,stbi_failure_reason());return 1; }else PIXEL_Fitting( stb, lx, ly ); int32_t **s=NULL, **d=NULL, **w=NULL, **t=NULL, **g=NULL, **b=NULL, st=1; const char *err_msg="Memory allocation failure"; if( !(s=(int32_t**)alloc2Darray(sizeof(int32_t),lx,ly)) ){ perror(err_msg); goto GOOD_BY; } if( !(d=(int32_t**)alloc2Darray(sizeof(int32_t),lx,ly)) ){ perror(err_msg); goto GOOD_BY; } if( !(w=(int32_t**)alloc2Darray(sizeof(int32_t),lx,ly)) ){ perror(err_msg); goto GOOD_BY; } if( !(t=(int32_t**)alloc2Darray(sizeof(int32_t),lx,ly)) ){ perror(err_msg); goto GOOD_BY; } if( !(g=(int32_t**)alloc2Darray(sizeof(int32_t),lx,ly)) ){ perror(err_msg); goto GOOD_BY; } if( !(b=(int32_t**)alloc2Darray(sizeof(int32_t),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] = nega ? (PIXEL){255-stb[e].B,255-stb[e].G,255-stb[e].R,0}:stb[e]; g[y][x] = gray = (306*ps[y][x].R + 601*ps[y][x].G + 117*ps[y][x].B)>>10; b[y][x] = gray>127 ? 255 : 0; pt[y][x] = (PIXEL){ gray, gray, gray, 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; } #endif //TEMPLATES_USE_AS_INCLUDE_FILE #endif //IMAGE_PROCESSING_TEMPLATES_C99
AI
double start = get_micro_time();
// Perform some image processing
for (int y = 0; y < ly; y++) for (int x = 0; x < lx; x++) d[y][x] = s[y][x];
double end = get_micro_time();
print_elapsed_time(end - start); // Output: Elapsed time: 0.123 msec
fprintf_elapsed_time(stdout, end - start, 1); // Output: Elapsed time: 0h 0m 0s 0ms 123usPIXEL **img = (PIXEL **)alloc2Darray(sizeof(PIXEL), 100, 100);
img[0][0] = (PIXEL){255, 128, 0, 255}; // Orange pixel with full opacity
POINTf p = (POINTf){10.5, 20.5}; // Floating-point coordinate
COMPLEX c = (COMPLEX){1.0, 2.0}; // Complex number: 1 + 2iint **array = (int **)alloc2Darray(sizeof(int), 200, 200);
PROCESS_ALL_PIXELS array[y][x] = x + y; // Initialize with gradient
free(array); // Must free the base pointerPIXEL **img = (PIXEL **)alloc2Darray(sizeof(PIXEL), 100, 100);
PROCESS_ALL_PIXELS img[y][x] = (PIXEL){255, 0, 0, 255}; // Red image
saveArray2bmp32("output.bmp", img, 100, 100);
free(img);COMPLEX **data = (COMPLEX **)alloc2Darray(sizeof(COMPLEX), 256, 256);
PROCESS_ALL_PIXELS data[y][x] = (COMPLEX){rand() % 256, 0}; // Random real data
FourierTransform(data, 256, 256, FFT); // Perform 2D FFT
PIXEL **out = (PIXEL **)alloc2Darray(sizeof(PIXEL), 256, 256);
PowerSpectrum("spectrum.bmp", data, out, 256, 256); // Save power spectrum
free(data); free(out);double **A = (double **)alloc2Darray(sizeof(double), 3, 3);
double **A_inv = (double **)alloc2Darray(sizeof(double), 3, 3);
A[0][0] = 1; A[0][1] = 2; A[0][2] = 3;
A[1][0] = 4; A[1][1] = 5; A[1][2] = 6;
A[2][0] = 7; A[2][1] = 8; A[2][2] = 10;
inverse_matrix(A, A_inv, 3);
dump_array(stdout, "%.2f ", "Inverse Matrix", A_inv, 3);
free(A); free(A_inv);PIXEL **src = (PIXEL **)alloc2Darray(sizeof(PIXEL), 100, 100);
PIXEL **dst = (PIXEL **)alloc2Darray(sizeof(PIXEL), 120, 120);
POINTf src_pts[4] = {{0, 0}, {100, 0}, {100, 100}, {0, 100}};
POINTf dst_pts[4] = {{10, 10}, {110, 10}, {100, 100}, {20, 110}};
Homography(src, 100, 100, NULL, src_pts, dst, 120, 120, dst_pts);
saveArray2bmp32("transformed.bmp", dst, 120, 120);
free(src); free(dst);int **binary = (int **)alloc2Darray(sizeof(int), 100, 100);
int **labels = (int **)alloc2Darray(sizeof(int), 100, 100);
PROCESS_ALL_PIXELS binary[y][x] = (y < 50 && x < 50) ? 255 : 0; // Top-left quadrant
int num_labels = labeling(binary, labels, 100, 100, 8);
printf("Number of components: %d\n", num_labels);
save_label_data("labels.bin", labels, 100, 100);
free(binary); free(labels);float **src = (float **)alloc2Darray(sizeof(float), 100, 100);
float **dst = (float **)alloc2Darray(sizeof(float), 100, 100);
float kernel[3][3] = {{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}}; // Blur
PROCESS_ALL_PIXELS src[y][x] = rand() % 256;
convolution_float(src, dst, 100, 100, 3, kernel);
free(src); free(dst);char *key = "example";
uint32_t hash = MurmurHash2(key, strlen(key), 0);
HashDB *db = make_HashDB256M();
db[hash >> 4].key = strdup(key);
db[hash >> 4].value = 42;
printf("Hash: %u, Value: %lld\n", hash, db[hash >> 4].value);
free(db[hash >> 4].key); free(db);PIXEL **canvas = (PIXEL **)alloc2Darray(sizeof(PIXEL), 100, 100);
PIXEL red = (PIXEL){255, 0, 0, 255};
DrawCircle(canvas, 100, 100, 50, 50, 20, red); // Red circle at (50,50)
saveArray2bmp32("circle.bmp", canvas, 100, 100);
free(canvas);int **src = (int **)alloc2Darray(sizeof(int), 100, 100);
int **dst = (int **)alloc2Darray(sizeof(int), 100, 100);
PROCESS_ALL_PIXELS src[y][x] = rand() % 256;
image_not(src, dst, 100, 100); // Invert pixel values
PIXEL **out = (PIXEL **)alloc2Darray(sizeof(PIXEL), 100, 100);
clamp_and_convert(dst, out, 100, 100);
saveArray2bmp32("inverted.bmp", out, 100, 100);
free(src); free(dst); free(out);./Template -v -n input.bmp output.bmp// In image_processing_program:
int image_processing_program(int **s, int **d, int **w, int **t, int **g, int **b, int lx, int ly) {
image_not(s, d, lx, ly); // Invert input image
return 0;
}USER
日本語で回答してください、
AI
double start = get_micro_time();
// 画像処理の実行
for (int y = 0; y < ly; y++) for (int x = 0; x < lx; x++) d[y][x] = s[y][x];
double end = get_micro_time();
print_elapsed_time(end - start); // 出力例: Elapsed time: 0.123 msec
fprintf_elapsed_time(stdout, end - start, 1); // 出力例: Elapsed time: 0h 0m 0s 0ms 123usPIXEL **img = (PIXEL **)alloc2Darray(sizeof(PIXEL), 100, 100);
img[0][0] = (PIXEL){255, 128, 0, 255}; // オレンジ色のピクセル
POINTf p = (POINTf){10.5, 20.5}; // 浮動小数点の座標
COMPLEX c = (COMPLEX){1.0, 2.0}; // 複素数: 1 + 2iint **array = (int **)alloc2Darray(sizeof(int), 200, 200);
PROCESS_ALL_PIXELS array[y][x] = x + y; // グラデーションで初期化
free(array); // ベースポインタを解放PIXEL **img = (PIXEL **)alloc2Darray(sizeof(PIXEL), 100, 100);
PROCESS_ALL_PIXELS img[y][x] = (PIXEL){255, 0, 0, 255}; // 赤い画像
saveArray2bmp32("output.bmp", img, 100, 100);
free(img);COMPLEX **data = (COMPLEX **)alloc2Darray(sizeof(COMPLEX), 256, 256);
PROCESS_ALL_PIXELS data[y][x] = (COMPLEX){rand() % 256, 0}; // ランダムな実数データ
FourierTransform(data, 256, 256, FFT); // 2D FFT
PIXEL **out = (PIXEL **)alloc2Darray(sizeof(PIXEL), 256, 256);
PowerSpectrum("spectrum.bmp", data, out, 256, 256); // パワースペクトル保存
free(data); free(out);double **A = (double **)alloc2Darray(sizeof(double), 3, 3);
double **A_inv = (double **)alloc2Darray(sizeof(double), 3, 3);
A[0][0] = 1; A[0][1] = 2; A[0][2] = 3;
A[1][0] = 4; A[1][1] = 5; A[1][2] = 6;
A[2][0] = 7; A[2][1] = 8; A[2][2] = 10;
inverse_matrix(A, A_inv, 3);
dump_array(stdout, "%.2f ", "逆行列", A_inv, 3);
free(A); free(A_inv);PIXEL **src = (PIXEL **)alloc2Darray(sizeof(PIXEL), 100, 100);
PIXEL **dst = (PIXEL **)alloc2Darray(sizeof(PIXEL), 120, 120);
POINTf src_pts[4] = {{0, 0}, {100, 0}, {100, 100}, {0, 100}};
POINTf dst_pts[4] = {{10, 10}, {110, 10}, {100, 100}, {20, 110}};
Homography(src, 100, 100, NULL, src_pts, dst, 120, 120, dst_pts);
saveArray2bmp32("transformed.bmp", dst, 120, 120);
free(src); free(dst);int **binary = (int **)alloc2Darray(sizeof(int), 100, 100);
int **labels = (int **)alloc2Darray(sizeof(int), 100, 100);
PROCESS_ALL_PIXELS binary[y][x] = (y < 50 && x < 50) ? 255 : 0; // 左上領域
int num_labels = labeling(binary, labels, 100, 100, 8);
printf("連結成分数: %d\n", num_labels);
save_label_data("labels.bin", labels, 100, 100);
free(binary); free(labels);float **src = (float **)alloc2Darray(sizeof(float), 100, 100);
float **dst = (float **)alloc2Darray(sizeof(float), 100, 100);
float kernel[3][3] = {{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}}; // ぼかし
PROCESS_ALL_PIXELS src[y][x] = rand() % 256;
convolution_float(src, dst, 100, 100, 3, kernel);
free(src); free(dst);char *key = "example";
uint32_t hash = MurmurHash2(key, strlen(key), 0);
HashDB *db = make_HashDB256M();
db[hash >> 4].key = strdup(key);
db[hash >> 4].value = 42;
printf("ハッシュ: %u, 値: %lld\n", hash, db[hash >> 4].value);
free(db[hash >> 4].key); free(db);PIXEL **canvas = (PIXEL **)alloc2Darray(sizeof(PIXEL), 100, 100);
PIXEL red = (PIXEL){255, 0, 0, 255};
DrawCircle(canvas, 100, 100, 50, 50, 20, red); // (50,50)に赤い円
saveArray2bmp32("circle.bmp", canvas, 100, 100);
free(canvas);int **src = (int **)alloc2Darray(sizeof(int), 100, 100);
int **dst = (int **)alloc2Darray(sizeof(int), 100, 100);
PROCESS_ALL_PIXELS src[y][x] = rand() % 256;
image_not(src, dst, 100, 100); // ピクセル値を反転
PIXEL **out = (PIXEL **)alloc2Darray(sizeof(PIXEL), 100, 100);
clamp_and_convert(dst, out, 100, 100);
saveArray2bmp32("inverted.bmp", out, 100, 100);
free(src); free(dst); free(out);./Template -v -n input.bmp output.bmp// image_processing_program内:
int image_processing_program(int **s, int **d, int **w, int **t, int **g, int **b, int lx, int ly) {
image_not(s, d, lx, ly); // 入力画像を反転
return 0;
} 使用例
typedef struct COMPLEX { double r, i; } COMPLEX;
#define v4D(type, lx, ly, lz, name) ((type(*)[(lx)][(ly)][(lz)])name)
#define v4Dcpx(lx, ly, lz, name) v4D(COMPLEX, lx, ly, lz, name)
void func_4D(COMPLEX *s, int lx, int ly, int lz, int lt) {
#undef A_s // s は、COMPLEX型の 1次元配列
#define A_s (v4Dcpx(lx, ly, lz, s)) // A_s は、lx列ly行lz層lt時間を持つ4次元配列
// 例: A_s[t][z][y][x].r で実数成分にアクセス
A_s[t][z][y][x].i = 0.5; // 配列サイズが、lt, lz, ly, lx の複素データ 虚数成分にアクセス
}
補足用
void func_4Dp(COMPLEX ****s, int lx, int ly, int lz, int lt) {
:
// 例: s[t][z][y][x].r で実数成分にアクセス
s[t][z][y][x].i = 0.5; // 配列サイズが、lt, lz, ly, lx の複素データ 虚数成分にアクセス
}
二つの関数の違いが判らず不思議に思う方は、生成AI に尋ねると、詳しく説明してくれます。
確かに、func_4D と func_4Dp の違いは少し微妙で、最初は理解しにくいかもしれません。
これらの関数の違いをもう少し詳しく掘り下げて説明します。
func_4D と func_4Dp の基本的な違いfunc_4Dfunc_4D では、s は 1次元配列のポインタ です。しかし、v4D マクロを使って、s を仮想的に
4次元配列 としてアクセスできるようにします。
v4Dcpx(lx, ly, lz, s) によって、s を COMPLEX 型の4次元配列のように扱うことができますが、
実際にはメモリは1次元配列として確保されています。
A_s[t][z][y][x] のようにアクセスし、r や i の成分にアクセスします。
これは 仮想的な4次元配列 にアクセスしているということです。
func_4Dpfunc_4Dp では、s は 4次元配列を指すポインタのポインタのポインタのポインタ です。 つまり、s はすでに実際に 4次元配列 のポインタを指しており、s[t][z][y][x] のようにインデックスで直接アクセスできます。
ここでの s は、COMPLEX[lt][lz][ly][lx] のように4次元の配列が実際にメモリ上に格納されているという点で、func_4D
とは異なります。
配列の渡し方:
func_4D: 配列 s は 1次元配列 として渡され、v4D マクロで4次元配列のようにアクセスします。
実際にはメモリは1次元配列として扱われます。
func_4Dp: 配列 s は 4次元配列を指すポインタのポインタのポインタのポインタ として渡され、
s[t][z][y][x] のように直接4次元配列にアクセスします。
メモリの構造:
func_4D: メモリ上では1次元配列として格納され、アクセスの際にマクロを通じて仮想的に4次元アクセスします。
func_4Dp: 実際に4次元配列がメモリ上に格納されており、ポインタを通じて4次元配列にアクセスします。
可読性:
func_4D: マクロを使って仮想的に4次元配列を扱っているため、少し分かりにくいかもしれませんが、
配列の管理は比較的簡単です。
func_4Dp: 4次元配列を直接扱っているため、アクセスは直感的で理解しやすいです。
パフォーマンス:
func_4D: マクロを使うことでインデックス計算を動的に行っているため、インデックス計算が少し複雑になりますが、
メモリ管理は比較的単純です。
func_4Dp: ポインタの参照を直接使うため、アクセスのパフォーマンスはシンプルで効率的です。
ただし、ポインタの多重参照（****s）が少しオーバーヘッドになる可能性があります。
func_4D は 1次元配列 を仮想的に 4次元配列 として扱う方法で、マクロによってアクセスします。
この方法はメモリ管理がシンプルであり、メモリを1次元で確保しながら、配列アクセスを4次元のように振る舞わせます。
func_4Dp は、4次元配列 を直接操作する方法で、ポインタの多重参照を使って直接アクセスします。
この方法は、より直感的で可読性が高いですが、ポインタ操作において少し複雑になることがあります。
結論としては、どちらの方法を選ぶかは、メモリ管理の方法やコードの可読性、パフォーマンス要件によって決まります。