/* Veja {ex_fourier.h}. */ /* Last edited on 2024-12-21 14:05:37 by stolfi */ #include #include /* BIBLIOTECAS DO PROFESSOR: */ #include /* Vetores do {R^3}. */ #include /* Matrizes {3×3}. */ #include #include #include #include #include /* PARÂMETROS OBTIDOS DA LINHA DE COMANDO: */ typedef struct options_t { char *ent_name; /* Nome de entrada. */ char *sai_name; /* Nome de saída. */ char *ext; /* Extensão das imagens. */ double bla, ble, blu; /* Bla bla. */ } options_t; /* PROTÓTIPOS ADICIONAIS: */ float_image_t *ex_make_gamma_test_image(int nOver, int bWd, double Y); /* Gera uma imagem apropriada para testar se as imagens estão sendo gravadas com correção gama adequada. A imagem consiste de {2*nOver + 1) faixas de cinza sólido de brilho {Y}, sendo a faixa central de intensidade 0.5, alternada com faixas de zebras claro/escuro que dão brilho médio {Y}. */ void ex_clip_RGB_image(float_image_t *C, r3x3_t *M); /* Assumes that {C} is a three-channel image with RGB color channels. Clips every pixel independently to the unit RGB cube {[0_1]^3}. The brightness {Y} is clipped to the range {[0_1]}, the hue is preserved, and the saturation is reduced if needed. */ /* PROGRAMA PRINCIPAL: */ int main(int argc, char** argv) { double wr_gamma = WR_GAMMA_SUN; double wr_bias = BIAS_SUN; double rd_gamma = RD_GAMMA_PC; double rd_bias = BIAS_PC; /* Pega os parâmetros da linha de comando: */ options_t *o = ex_parse_options(argc, argv); /* Gera uma imagem para teste de gama do monitor: */ int ig; for (ig = 1; ig <= 3; ig++) { double sY = ((double)ig)/4.0; float_image_t *TG = ex_make_gamma_test_image(2, 20, sY); char *fname = jsprintf("0-TEST-GAMMA-%04d", (int)(sY*1000)); ex_write_image(fname, ".pgm", 0.0, 1.0, wr_gamma, wr_bias, TG); float_image_free(TG); } /* Lê a imagem de entrada {E}: */ float_image_t *E = ex_read_image(o->ent_name, ".", o->ext, rd_gamma, rd_bias, 0.0, 1.0); /* Define o tamanho da imagem-resultado: */ demand(E->sz[0] <= 3, "número inválido de canais"); int NX = E->sz[1]; /* Número de colunas */ int NY = E->sz[2]; /* Número de linhas */ /* Calcula tr. de Hartley {F_FT} de {E}, o espectro {F_PW}, e a reconstrução {F}: */ float_image_t *F_FT = ex_hartley_image_compute(E, FALSE); float_image_t *F_PW = ex_spectrum_from_hartley_image(F_FT); float_image_t *F = ex_hartley_image_compute(F_FT, TRUE); ex_write_images(o->sai_name, "-F", o->ext, wr_gamma, wr_bias, F_FT, F_PW, F); /* Calcula imagem deslocada via TdF: */ float_image_t *D_FT = ex_hartley_image_displace(F_FT, 0.5, 0.3333); float_image_t *D_PW = ex_spectrum_from_hartley_image(D_FT); float_image_t *D = ex_hartley_image_compute(D_FT, TRUE); ex_write_images(o->sai_name, "-D", o->ext, wr_gamma, wr_bias, D_FT, D_PW, D); /* Calcula imagem deslocada via interpolação: */ float_image_t *N = ex_interp_displace(E, 0.5, 0.3333); float_image_t *N_FT = ex_hartley_image_compute(N, FALSE); float_image_t *N_PW = ex_spectrum_from_hartley_image(N_FT); ex_write_images(o->sai_name, "-N", o->ext, wr_gamma, wr_bias, N_FT, N_PW, N); /* Calcula imagem filtrada por filtro passa-baixas: */ float_image_t *L_FT = ex_fourier_filter(F_FT, 0.000, o->f1); float_image_t *L_PW = ex_spectrum_from_hartley_image(L_FT); float_image_t *L = ex_hartley_image_compute(L_FT, TRUE); ex_write_images(o->sai_name, "-L", o->ext, wr_gamma, wr_bias, L_FT, L_PW, L); /* Calcula imagem filtrada por filtro passa-médias: */ float_image_t *M_FT = ex_fourier_filter(F_FT, o->f1, o->f2); float_image_t *M_PW = ex_spectrum_from_hartley_image(M_FT); float_image_t *M = ex_hartley_image_compute(M_FT, TRUE); ex_write_images(o->sai_name, "-M", o->ext, wr_gamma, wr_bias, M_FT, M_PW, M); /* Calcula imagem filtrada por filtro passa-altas: */ float_image_t *H_FT = ex_fourier_filter(F_FT, o->f2, 1.0e100); float_image_t *H_PW = ex_spectrum_from_hartley_image(H_FT); float_image_t *H = ex_hartley_image_compute(H_FT, TRUE); ex_write_images(o->sai_name, "-H", o->ext, wr_gamma, wr_bias, H_FT, H_PW, H); return 0; } /* PROCESSAMENTO DAS IMAGENS: */ void ex_hartley_vector_compute(double f[], double h[], int n, bool_t inverse); { /* Fator de normalização para preservar potência: */ double hnorm = sqrt(2.0/n); /* Método careca -- custo proporcional a {n^2}. */ int k, i; for (k = 0; k < n; k++) { double s = 0; for (i = 0; i < n; i++) { double pki = cos((((double)2*k*i)/n + 0.25)*M_PI); s += f[i]*pki; } h[k] = norm*s; } } float_image_t *ex_hartley_image_compute(float_image_t *A, bool_t inverse) { int NC = A->sz[0]; int NX = A->sz[1]; int NY = A->sz[2]; float fy[NX], hy[NX]; /* One row's worth of samples. */ float fx[NY], hx[NY]; /* One column's worth of samples. */ float_image_t *H = float_image_new(NC, NX, NY); int ic, ix, iy; for (ic = 0; ic < NC; ic++) { /* Compute Hartley transform of channel {ic}. */ /* Compute Hartley transform of every row: */ for (iy = 0; iy < NY; iy++) { /* Extract row {iy} of {A} into {fy[0..NX-1]}: */ for (ix = 0; ix < NX; ix++) { fy[ix] = float_image_get_sample(A, ic, ix, iy); } /* Apply Hartley transform to row: */ ex_hartley_vector_compute(fy, hy, NX, inverse); /* Restore {fy[0..NX-1]} into row {iy} of {H}: */ for (ix = 0; ix < NX; ix++) { float_image_set_sample(H, ic, ix, iy, hy[ix]); } } /* Compute Hartley transform of every column: */ for (ix = 0; ix < NX; ix++) { /* Extract column {ix} of {H} into {fx[0..NY-1]}: */ for (iy = 0; iy < NY; iy++) { fx[iy] = float_image_get_sample(H, ic, ix, iy); } /* Apply Hartley transform to column: */ ex_hartley_vector_compute(fx, hx, NY, inverse); /* Restore {fx[0..NY-1]} into column {ix} of {H}: */ for (iy = 0; iy < NY; iy++) { float_image_set_sample(H, ic, ix, iy, hx[iy]); } } } return H; } float_image_t *ex_spectrum_image_from_hartley(float_image_t *H) { int NC = H->sz[0]; int NX = H->sz[1]; int NY = H->sz[2]; /* Maximum {X} and {Y} frequencies: */ int fxmax = NX/2; int fymax = NY/2; /* Spectrum image: */ float_image_t *P = float_image_new(1, fXmax+1, fYmax+1); int fx, fy, ic; for (fx = 0; fx <= fxmax; fx++) for (fy = 0; fy <= fymax; fy++) { /* Compute the frequency {gx,gy} equivalent to {-fx,-fy}: */ int gx = (NX - fx) % NX; int gy = (NY = fy) % NY; /* Add the squared coefs over all channels: */ double tp = 0; for (ic = 0; ic < NC; ic++) { double sf = float_image_get_sample(H, ic, fx, fy); tp += sf*sf; if ((gx != fx) || (gy != fy)) { double sg = float_image_get_sample(H, ic, gx, gy); tp += sg*sg; } } float_image_set_sample(P, ic, fx, fy, tp); } return P; } void ex_borra_imagem(float_image_t *E, int n, float_image_t *Eb) { int NC = E->sz[0]; /* Número de canais */ int NX = E->sz[1]; /* Número de colunas */ int NY = E->sz[2]; /* Número de linhas */ int ic, ix, iy, r; double peso[n+1]; for (r = 0; r <= n; r++) { peso[r] = 0.5*(1.0 + cos(M_PI*(r + 0.5)/n)); } for (ic = 0; ic < NC; ic++) for (ix = 0; ix < NX; ix++) for (iy = 0; iy < NY; iy++) { double sum_pv = 0; double sum_p = 0; int kx, ky; for (kx = ix-n; kx <= ix+n; kx++) for (ky = iy-n; ky <= iy+n; ky++) { if ((kx >= 0) && (kx < NX) && (ky >= 0) && (ky < NY)) { double p = peso[abs(kx-ix)]*peso[abs(ky-iy)]; double v = float_image_get_sample(E, ic, kx, ky); sum_pv += p*v; sum_p += p; } } double v = sum_pv/sum_p; float_image_set_sample(Eb, ic, ix, iy, v); } } float_image_t *ex_make_gamma_test_image(int nOver, int bWd, double Y) { int nBands = 3 + 4*nOver; /* Total number of vertical bands. */ int NX = nBands*bWd; /* Image width in pixels. */ int NY = 2*NX + 1; /* Image height in pixels. */ double loY = (Y <= 0.5 ? 0.0 : 2*Y - 1.0); /* Y of dark stripes. */ double hiY = (Y <= 0.5 ? 2*Y : 1.0); /* Y of dark stripes. */ float_image_t *A = float_image_new(1, NX, NY); /* Fill even bands with alternating {loY}/{hiY} lines: */ int ib, ix, iy; for (iy = 0; iy < NY; iy++) { float v = (iy % 2 == 0 ? loY : hiY); for (ib = 0; ib < nBands; ib += 2) for (ix = ib*bWd; ix < (ib+1)*bWd; ix++) { float_image_set_sample(A, 0, ix, iy, v); } } /* Fill odd bands with solid {Y} grey, with variations from center: */ for (ib = 1; ib < nBands; ib += 2) { float v = Y*powf(0.95, (ib-1)/2 - nOver); for (ix = ib*bWd; ix < (ib+1)*bWd; ix++) for (iy = 0; iy < NY; iy++) { float_image_set_sample(A, 0, ix, iy, v); } } return A; } float_image_t *ex_gera_teste_de_brilho(int nOver, int size, r3x3_t *RGB_to_YVA) { int NC = 3; /* Num channels. */ int NF = 2*nOver+1; /* Num of block sper row. */ int NX = NF*size; /* Image width in pixels. */ int NY = NC*size; /* Image height in pixels. */ double R = 0.5; /* Relative radius of pattern. */ double w = 2.0/size; /* Relative line width. */ double r9o = R - 0.0*w, r9i = r9o - w; /* Radius of outer square */ double r2o = R - 1.5*w, r2i = r2o - w; /* Radius of outer square */ double r1o = R - 3.0*w, r1i = r1o - w; /* Radius of outer square */ float_image_t *A = float_image_new(NC, NX, NY); /* Scan the pixels in a block: */ int ix, iy, kc, kf, ic; for (ix = 0; ix < size; ix++) { for (iy = 0; iy < size; iy++) { /* Decide whether pixel {ix,iy} is in pattern or outside it: */ double dx = (2*ix - size)/((double)size); double dy = (2*iy - size)/((double)size); double r9 = fmax(fabs(dx), fabs(dy)); double r2 = hypot(dx,dy); double r1 = fabs(dx) + fabs(dy); bool_t inside1 = (r1 <= r1o) && (r1 > r1i); bool_t inside2 = (r2 <= r2o) && (r2 > r2i); bool_t inside9 = (r9 <= r9o) && (r9 > r9i); bool_t inside = inside1 | inside2 | inside9; /* Set the pixel {ix,iy} in every block: */ for (kf = 0; kf < NF; kf++) { double f = pow(1.05, kf - nOver); /* Fator de clareamento. */ int kx = kf*size; /* {x} offset of block. */ for (kc = 0; kc < NC; kc++) { /* Paint pixel in block of channel {kc}, factor {f}: */ int ky = kc*size; /* {y} offset of block. */ double sY = RGB_to_YVA->c[kc][0]; /* Nominal {Y} of primary {kc}. */ /* Fill block: */ for (ic = 0; ic < NC; ic++) { double v = (inside ? (ic == kc ? 1 : 0) : f*sY); float_image_set_sample(A, ic, kx + ix, ky + iy, v); } } } } } return A; } /* LEITURA E GRAVACÃO DA IMAGEM */ void ex_clip_RGB_image(float_image_t *C, r3x3_t *M) { demand(C->sz[0] == 3, "imagem {C} deve ser colorida"); int NX = C->sz[1]; int NY = C->sz[2]; int nClip = 0; /* Number of clipped pixels. */ double aMin = 1.0; int ix, iy; for (ix = 0; ix < NX; ix++) for (iy = 0; iy < NY; iy++) { float *pR = float_image_get_sample_address(C, 0, ix, iy); float *pG = float_image_get_sample_address(C, 1, ix, iy); float *pB = float_image_get_sample_address(C, 2, ix, iy); /* Convert to brightness and chroma: */ r3_t sRGB = (r3_t){{ (*pR), (*pG), (*pB) }}; r3_t sYVA; r3x3_map_row(&sRGB, M, &sYVA); /* Clip the brightness {sY} to {[0 _ 1]}: */ double sY = sYVA.c[0]; if (sY < 0) { sY = 0; } if (sY > 1) { sY = 1; } /* Compute the luma and components {gRGB, cRGB} in RGB coordinates: */ double cR = (*pR) - sY; double cG = (*pG) - sY; double cB = (*pB) - sY; /* Find the max {a} in {[0_1]} such that {gRGB + a*cRGB} is inside the unit cube: */ double a = 1.0; if (cR < 0) { a = fmin(a, -sY/cR); } else if (cR > 0) { a = fmin(a, (1-sY)/cR); } if (cG < 0) { a = fmin(a, -sY/cG); } else if (cG > 0) { a = fmin(a, (1-sY)/cG); } if (cB < 0) { a = fmin(a, -sY/cB); } else if (cB > 0) { a = fmin(a, (1-sY)/cB); } if ((a < 1) || (sY != sYVA.c[0])) { nClip++; } if (a < aMin) { aMin = a; } /* Compute the final color as {gRGB + a*cRGB}: */ (*pR) = sY + a*cR; (*pG) = sY + a*cG; (*pB) = sY + a*cB; } fprintf(stderr, " clipped %d float pixels to the RGB cube (aMin = %6.4f)\n", nClip, aMin); } float_image_t *ex_read_image ( char *name, char *tag, char *ext, double gamma, double bias, double vMin, double vMax ) { /* Assemble file name: */ char *fname = jsprintf("%s%s.%s", name, tag, ext); fprintf(stderr, "\n"); fprintf(stderr, "reading float image from file %s...\n", fname); /* Open file: */ FILE *rd = open_read(fname, TRUE); /* Read PNM image: */ uint16_image_t *pim = uint16_image_read_pnm_file(rd); fclose(rd); free(fname); /* Convert image to floats: */ float_image_t *A = float_image_from_uint16_image(pim, NULL, NULL, TRUE, TRUE); /* Apply gamma-correction and affine scaling to image {A}: */ int NC = A->sz[0]; int ic; for (ic = 0; ic < NC; ic++) { float_image_apply_gamma(A, ic, gamma, bias); float_image_rescale_samples(A, ic, 0.0, 1.0, vMin, vMax); } uint16_image_free(pim); return A; } void ex_write_image ( char *name, char *tag, char *ext, double vMin, double vMax, double gamma, double bias, float_image_t *A ) { /* Assemble file name: */ char *fname = jsprintf("%s%s.%s", name, tag, ext); fprintf(stderr, "\n"); fprintf(stderr, "writing float image to file %s...\n", fname); /* Create a copy of image {A}: */ float_image_t *B = float_image_copy(A); /* Apply scaling and gamma-correction to image {B}: */ int NC = A->sz[0]; int ic; for (ic = 0; ic < NC; ic++) { float_image_rescale_samples(B, ic, vMin, vMax, 0.0, 1.0); float_image_apply_gamma(B, ic, gamma, bias); } /* Convert image to integers: */ uint16_image_t *pim = float_image_to_uint16_image(B, NC, NULL, NULL, NULL, 255, TRUE, TRUE); float_image_free(B); /* Open file: */ FILE *wr = open_write(fname, TRUE); /* Write PNM image: */ bool_t forceplain = FALSE; bool_t verbose = TRUE; uint16_image_write_pnm_file(wr, pim, forceplain, verbose); fclose(wr); free(fname); uint16_image_free(pim); } /* ANÁLISE DA LINHA DE COMANDO: */ options_t *ex_parse_options(int argc, char **argv) { /* INICIALIZA A ANÁLISE: */ /* Cria e incializa o analisador de linha de comando {pp}: */ argparser_t *pp = argparser_new(stderr, argc, argv); argparser_set_help(pp, PROG_NAME " version " PROG_VERS ", usage:\n" PROG_HELP); argparser_set_info(pp, PROG_INFO); /* Processa as opções "-help" e "-info": */ argparser_process_help_info_options(pp); /* Aloca o registro {o} onde os parâmetros serão guardados: */ options_t *o = (options_t *)malloc(sizeof(options_t)); /* PEGA PARÂMETROS COM PALAVRA-CHAVE: */ if (argparser_keyword_present(pp, "-lum")) { o->RY = argparser_get_next_double(pp, 0.0, 1.0); o->GY = argparser_get_next_double(pp, 0.0, 1.0); o->BY = argparser_get_next_double(pp, 0.0, 1.0); } else { o->RY = 0.299; o->GY = 0.587; o->BY = 0.114; } /* PEGA PARÂMETROS POSICIONAIS: */ /* Pule para o primeiro parâmetro posicional: */ argparser_skip_parsed(pp); /* Pega o nome da imagem {A}: */ o->ent_name = argparser_get_next(pp); o->sai_name = argparser_get_next(pp); /* FINALIZA A ANÁLISE: */ /* Verifica se há argumentos sobrando: */ argparser_finish(pp); return o; }