/*  2025/06/15 初版：機能強化版 C99画像処理用フレームワーク

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

画像処理の学習、アルゴリズム検証用の強化版フレームワーク、 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 のループ変数は、自動的に     //
//        プライベートになるので宣言は省略できるが、明示した方が誤解を回避できる     //
//                                                                                   //
#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>

#define HERE                            fprintf(stderr,"% 6d: %s\n",__LINE__, __func__);
#define IF                              if          // ソースコードの行間での位置合わせ用
#define EI                              else if     // 縦に揃えて可読性を向上させるため。
#define EE                              else        // 人間の可読性向上用エイリアス。
//////////////////////////////////////////////////////////////////////////////////////////
// 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;
}
//  #undef DrawText    This is necessary if you want to implement DrawText yourself.
//  #undef DrawTextW   This is necessary if you want to implement DrawTextW yourself.
#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);
void print_elapsed_time(double microsec) {
    IF(microsec < 1000.0) printf("Elapsed time: %.3f usec\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 { 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)

                                    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++)

void **alloc2Darray( int, int, int );
#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>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#undef  PI
#define PI  3.1415926535897932384626433832795L  // 円周率(π)

typedef struct COMPLEX { double r, i; } COMPLEX;

#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 AppendCOMPLEX(a,b) { (a).r += (b).r, (a).i += (b).i ; } 

enum { DFT, IDFT };
enum { FFT=DFT, IFFT=IDFT };

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;
    }
}
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;
}
void fourierTransform_horizontal( COMPLEX **d, int elements, int lines, COMPLEX *rotor, int flag ) {
    if( ft_size_check(elements) ){// Processing by FFT
        #pragma omp parallel for
        for(int e=0; e<lines; e++ ){ COMPLEX  W, temp;
            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 ){ temp=d[e][j], d[e][j]=d[e][i], d[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 0rad.
                    W=rotor[ir]; 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(d[e][a], d[e][b]);  // For butterfly operation
                        d[e][a] = AddCOMPLEX(d[e][a], d[e][b]);  // a ->- (a + b) ->- a
                        d[e][b] = MulCOMPLEX(temp, W);           // b ->- (a-b)*W ->- b
            }   }   }   double reciprocal=1./elements;
            if(flag==IFFT)for(int i=0; i<elements; d[e][i]=ImeMulCOMPLEX(d[e][i],reciprocal),i++);
        }
    } else {// Processing by DFT
        #pragma omp parallel for
        for(int e=0; e<lines; e++ ){ COMPLEX *tmp=malloc( elements*sizeof(COMPLEX));
            if(!tmp){ perror("Memory allocation error for temp.\n"); exit(1); }
            for(int u=0; u<elements; u++ ){ COMPLEX t = (COMPLEX){0,0};
                for(int v=0;v<elements;v++) AppendCOMPLEX( t, MulCOMPLEX(d[e][v],rotor[(u*v)%elements]));
                tmp[u] = t;
            }  double reciprocal = 1./elements;
            if(flag==IDFT) for(int i=0; i<elements; i++ ) d[e][i] = ImeMulCOMPLEX(tmp[i],reciprocal);
            else           for(int i=0; i<elements; i++ ) d[e][i] = tmp[i];
            free(tmp);
        }
    }
}
int fourierTransform( COMPLEX **d, int lx, int ly, int flag ) {
    //       rotorは垂直処理で再利用。max(lx, ly)に十分なサイズで確保済み
    COMPLEX *rotor = (COMPLEX *)calloc( max(lx,ly), sizeof(COMPLEX) );
    if(!rotor){ fprintf(stderr,"Memory allocation failure"); return 1; }

    int     elements=lx;
    double  unit = (flag==IDFT) ? 2*PI/elements : -2*PI/elements;
    for(int e=0; e<elements; rotor[e]=RotorCOMPLEX(unit*e), e++ );
    
    fourierTransform_horizontal( d, elements, ly, rotor, flag );
    ////////////////////////////////////////////////////////////////////////////////////////////////////
    // Below is the vertical processing of the two-dimensional Fourier transform..    以下、縦方向の処理
    COMPLEX **transpose =(COMPLEX**)alloc2Darray(sizeof(COMPLEX),ly,lx);
    if((elements=ly)!=lx){ unit = (flag==IDFT) ? 2*PI/elements : -2*PI/elements;
        for(int e=0; e<elements; rotor[e]=RotorCOMPLEX(unit*e), e++);
    }// データーが連続配置されるとアクセス速度が著しく向上、可能な時は優先する
    if( transpose ){// Processing by FFT, Consecutive arrangement of data
        PROCESS_ALL_PIXELS transpose[x][y] = d[y][x];
        fourierTransform_horizontal( transpose, elements, lx, rotor, flag );
        PROCESS_ALL_PIXELS d[y][x] = transpose[x][y];
    } else if( ft_size_check(elements) ){
        #pragma omp parallel for
        for(int x=0; x<lx; x++ ){ COMPLEX  W, temp;   //  W: Rotation factor
            for(int i=0, j=1; j<elements-1; j++ ) {// Processing by FFT
                for(int k=elements>>1; k>(i^=k); k>>=1 );
                if( j<i ){ temp=d[j][x], d[j][x]=d[i][x], d[i][x]=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 rad.
                    W=rotor[ir]; 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(d[a][x], d[b][x]);    // For butterfly operation
                        d[a][x] = AddCOMPLEX(d[a][x], d[b][x]);    // a ->- (a + b) ->- a
                        d[b][x] = MulCOMPLEX(temp, W);             // b ->- (a-b)*W ->- b
            }   }   }   double reciprocal=1./elements;
            if(flag==IFFT)for(int i=0; i<elements; d[i][x]=ImeMulCOMPLEX(d[i][x],reciprocal),i++);
        }
    } else {// Processing by DFT
        #pragma omp parallel for
        for(int x=0; x<lx; x++ ){ COMPLEX *tmp=malloc( elements*sizeof(COMPLEX));
            if(!tmp){ perror("Memory allocation error for temp.\n"); exit(1); }
            for(int u=0; u<elements; u++ ){ COMPLEX t = (COMPLEX){0,0};
                for(int v=0;v<elements;v++) AppendCOMPLEX( t, MulCOMPLEX(d[v][x],rotor[(u*v)%elements]));
                tmp[u] = t;
             }  double reciprocal = 1./elements;
             if(flag==IDFT) for(int e=0; e<elements; e++ ) d[e][x] = ImeMulCOMPLEX(tmp[e],reciprocal);
             else           for(int e=0; e<elements; e++ ) d[e][x] = tmp[e];
            free(tmp);
        }
    }   
    free(transpose);

    free(rotor);
    return 0;
}
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 );
}
////////////////////////////////////////////////////////////////////////////////
// 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 ){
  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;  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;
}
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
//
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 ネガ画像：輝度反転画像）
    //  - 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
    Ex.  //if(verbose) printf("Using %d threads.\n", omp_get_max_threads());
    #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|
     
     // ラベリング     確保するデータ配列のデータ部、”連続エリアとして確保”されます。必須条件
     double start_time = get_micro_time();
     int labels = labeling( b, w, lx, ly, connectivity );
     double end_time   = get_micro_time();
     print_elapsed_time(end_time - start_time); // 処理時間とラベル数の表示
     printf("Number of objects (labels): %d\n", labels);

     if( save_binary && binary_file && save_label_data(binary_file, w, lx, ly) )
         fprintf(stderr, "Failed to save binary label data.\n"); // バイナリ保存（オプション）
     PROCESS_ALL_PIXELS d[y][x] = int2PseudoColor(w[y][x]);      // 疑似カラー変換
     save_without_conversion = 1;                                // 疑似カラー画像を直接保存
     */
  
    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;
                          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] [-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