#include #include #include #include #include #define filength 7 void wt1_mz (double *a, double *h, double *c, int N, int D, int scale); void filter (double *input, double *filtx, double *filty, double *output, int W, int H) { int i, j, kanaal; double *Row_In, *Row_Out, *Column_In, *Column_Out; int scale = 0; double Filt_Length = filength; Row_In = malloc (W * sizeof (double)); Column_In = malloc (H * sizeof (double)); Row_Out = malloc (W * sizeof (double)); Column_Out = malloc (H * sizeof (double)); for (j = 0; j < H; j++) { for (i = 0; i < W; i++) { int index = i + W * j; Row_In[i] = input[index]; } wt1_mz (Row_In, filtx, Row_Out, W, Filt_Length, scale); for (i = 0; i < W; i++) { int index = i + W * j; output[index] = Row_Out[i]; } } for (i = 0; i < W; i++) { for (j = 0; j < H; j++) { int index = i + W * j; Column_In[j] = output[index]; } wt1_mz (Column_In, filty, Column_Out, H, Filt_Length, scale); for (j = 0; j < H; j++) { int index = i + W * j; output[index] = Column_Out[j]; } } free (Row_In); free (Row_Out); free (Column_In); free (Column_Out); } void wt1_mz (double *a, double *h, double *c, int N, int D, int scale) { int i, j, Dext, I1, I2, len, t; double *ap, *he; t = pow (2, scale); ap = malloc (3 * N * sizeof (double)); Dext = (D - 1) * t + 1; he = malloc (Dext * sizeof (double)); /* Solving the border problems: Input vector should be made periodical and mirrored*/ for (i = 0; i < N; i++) { ap[i + N] = a[i]; ap[3 * N - 1 - i] = a[i]; ap[N - 1 - i] = a[i]; } /* Extending the impulse response h: Putting 2^p-1 zeros between the coefficients of h*/ for (i = 0; i < Dext; i++) he[i] = 0; for (i = 0; i < D; i++) he[i * t] = h[i]; /* Convolution: c=conv(ap,he) */ I2 = ((D - 1) / 2) * t; /*index of the element h(0) in the vector he */ for (j = 0; j < N; j++) { c[j] = 0; I1 = N + j; /*index of the element a(0+j) in the vector ap */ if ((I1 + I2 + 1) >= Dext) len = Dext; else len = I1 + I2 + 1; for (i = 0; i < len; i++) { c[j] += he[i] * ap[I1 + I2 - i]; } } free (ap); free (he); } double argum(double x, double y) { double out, temp1, temp2; //first quadrant if ((x > 0) && (y > 0)) {out = atan(y/x);} //second quadrant if ((x < 0) && (y > 0)) {out = M_PI + atan(y/x);} //third qudrant if ((x < 0) && (y < 0)) {out = M_PI + atan(y/x);} //fourth quadrant if ((x > 0) && (y < 0)) {out = 2 * M_PI + atan(y/x);} //limit cases if ((x == 0) && (y > 0)) {out = M_PI / 2;} if ((x == 0) && (y < 0)) {out = 3 * M_PI /2;} if ((x > 0) && (y == 0)) {out = 0.0;} //printf("out = %lf\n",out); return out; } output (double *imin, int W, int H, ImlibData *idout, ImlibImage *imout, char *nameout, int range) { int i,j; double minin, maxin; minin = maxin = 0.0; for (j=0; j maxin) {maxin = imin[index];} if (imin[index] < minin) {minin = imin[index];} } } if (range == 0) { for (j=0; jrgb_data[index] = imout->rgb_data[index + 1] = imout->rgb_data[index + 2] = 127 + (127 / (maxin - minin)) * imin[index2]; } } } else { for (j=0; jrgb_data[index] = imout->rgb_data[index + 1] = imout->rgb_data[index + 2] = (255 / (maxin - minin)) * (imin[index2] - minin); } } } Imlib_save_image (idout, imout, nameout, NULL); } int main (int argc, char **argv) { double g2f1[filength] = { 0.001933, 0.118119, 0.338927, -0.921300, 0.338927, 0.118119, 0.001933}; double g2f2[filength] = { 0.000123, 0.018316, 0.367879, 1.000000, 0.367879, 0.018316, 0.000123}; double g2f3[filength] = {-0.000503, -0.049745, -0.499580, 0.000000, 0.499580, 0.049745, 0.000503}; double h2f1[filength] = {-0.002443, -0.062551, 0.451172, 0.000000, -0.451172, 0.062551, 0.002443}; double h2f2[filength] = { 0.000123, 0.018316, 0.367879, 1.000000, 0.367879, 0.018316, 0.000123}; double h2f3[filength] = {-0.000370, -0.036631, -0.367879, 0.000000, 0.367879, 0.036631, 0.000370}; double h2f4[filength] = { 0.001018, 0.059498, 0.091418, -0.751500, 0.091418, 0.059498, 0.001018}; double g2f3m[9] = { 0.0028, 0.0480, 0.3020, 0.5806, 0.0000, -0.5806, -0.3020, -0.0480, -0.0028}; double lp2[9] = { 0.0000, 0.0000, 0.0038, 0.9923, 0.0038, 0.0000, 0.0000, 0.0000, 0.0000}; double theta; char scl[2], namerec[15], namelp; char nameH2a[15], nameH2b[15], nameH2c[15], nameH2d[15], name045[15], name135[15]; char nameG2a[15], nameG2b[15], nameG2c[15], namemag[15], nameori[15]; ImlibData *id, *idH2a, *idH2b, *idH2c, *idH2d, *id045, *id135, *idrec, *idori, *idmag, *idG2a, *idG2b, *idG2c; ImlibImage *imori, *immag; ImlibImage *im, *imH2a, *imH2b, *imH2c, *imH2d, *im045, *im135, *imG2a, *imG2b, *imG2c; FILE *fp, *fp2; double *input, *orig, *rec, *C1, *C2, *C3, *ori, *mag; double *H2a, *H2b, *H2c, *H2d, *t000, *t045, *t090, *t135, mmax, mmin, *G2a, *G2b, *G2c; double *G2aa, *G2bb, *G2cc, *lp, *lpaa; int W, H, i, j, k; double drempel; Display *disp; if (argc < 2) { printf("Usage: %s 'image_in' 'drempelwaarde (optioneel)'\n",argv[0],argv[1]); return 0; } if (argc > 2) {drempel = atof(argv[2]);} else {drempel = 5000.0;} disp = XOpenDisplay (NULL); id = idH2a = idH2b = idH2c = idH2d = idG2a = idG2b = idG2c = idori = idmag = id045 = idrec = Imlib_init (disp); //idori = idmag = Imlib_init (disp); im = Imlib_load_image (id, argv[1]); imH2a = Imlib_load_image (idH2a, argv[1]); imH2b = Imlib_load_image (idH2b, argv[1]); imH2c = Imlib_load_image (idH2c, argv[1]); imH2d = Imlib_load_image (idH2d, argv[1]); imG2a = Imlib_load_image (idH2b, argv[1]); imG2b = Imlib_load_image (idH2c, argv[1]); imG2c = Imlib_load_image (idH2d, argv[1]); imori = Imlib_load_image (idori, argv[1]); immag = Imlib_load_image (idmag, argv[1]); im045 = Imlib_load_image (id045, argv[1]); W = im->rgb_width; H = im->rgb_height; input = malloc (W * H * sizeof (double)); rec = malloc (W * H * sizeof (double)); G2a = malloc (W * H * sizeof (double)); G2b = malloc (W * H * sizeof (double)); G2c = malloc (W * H * sizeof (double)); G2aa = malloc (W * H * sizeof (double)); G2bb = malloc (W * H * sizeof (double)); G2cc = malloc (W * H * sizeof (double)); lp = malloc (W * H * sizeof (double)); lpaa = malloc (W * H * sizeof (double)); t045 = malloc (W * H * sizeof (double)); H2a = malloc (W * H * sizeof (double)); H2b = malloc (W * H * sizeof (double)); H2c = malloc (W * H * sizeof (double)); H2d = malloc (W * H * sizeof (double)); C1 = malloc (W * H * sizeof (double)); C2 = malloc (W * H * sizeof (double)); C3 = malloc (W * H * sizeof (double)); ori = malloc (W * H * sizeof (double)); mag = malloc (W * H * sizeof (double)); for (j = 0; j < H; j++) {for (i = 0; i < W; i++) {input[i + j * W] = im->rgb_data[3 * (i + j * W)];}} filter (input, h2f1, h2f2, H2a, W, H); filter (input, h2f4, h2f3, H2b, W, H); filter (input, h2f3, h2f4, H2c, W, H); filter (input, h2f2, h2f1, H2d, W, H); filter (input, g2f1, g2f2, G2a, W, H); filter (input, g2f3, g2f3, G2b, W, H); filter (input, g2f2, g2f1, G2c, W, H); //filter (G2a, g2f1, g2f2, G2aa, W, H); //filter (G2b, g2f3m,g2f3m, G2bb, W, H); //filter (G2c, g2f2, g2f1, G2cc, W, H); //filter (input, lp2, lp2, lp, W, H); //filter (lp, lp2, lp2, lpaa, W, H); theta = 45.0 * M_PI / 360.0; for (j = 0; j < H; j++) { for (i = 0; i < W; i++) { int index = i + j * W; C1[index] = 0.5 * pow(G2b[index],2) + 0.25 * G2a[index] * G2c[index] + 0.375 * (pow(G2a[index],2) + pow(G2c[index],2)) + 0.3125 * (pow(H2a[index],2) + pow(H2d[index],2)) + 0.5625 * (pow(H2b[index],2) + pow(H2c[index],2)) + 0.375 * ((H2a[index] * H2c[index]) + (H2b[index] * H2d[index])); C2[index] = 0.5 * (pow(G2a[index],2) - pow(G2c[index],2)) + 0.46875 * (pow(H2a[index],2) - pow(H2d[index],2)) + 0.28125 * (pow(H2b[index],2) - pow(H2c[index],2)) + 0.1875 * (H2a[index] * H2c[index] - H2b[index] * H2d[index]); C3[index] = - G2a[index] * G2b[index] - G2b[index] * G2c[index] - 0.9375 * (H2c[index] * H2d[index] + H2a[index] * H2b[index]) - 1.6875 * H2b[index] * H2c[index] - 0.1875 * H2a[index] * H2d[index]; /*t045[index] = cos(theta) * cos(theta) * G2a[index] - 2 * cos(theta) * sin(theta) * G2b[index] + sin(theta) * sin(theta) * G2c[index]; rec[index] = G2aa[index] + G2bb[index] + G2cc[index]; //+ 10 * lpaa[index];*/ } } mmax = mmin = 0; for (j = 0; j < H; j++) { for (i = 0; i < W; i++) { int index = i + j * W; mag[index] = sqrt(pow(C2[index], 2) + pow(C3[index], 2)); ori[index] = argum(C2[index], C3[index]) * 0.5; } } strcpy(nameH2a, "imH2a.tif"); strcpy(nameH2b, "imH2b.tif"); strcpy(nameH2c, "imH2c.tif"); strcpy(nameH2d, "imH2d.tif"); strcpy(nameG2a, "imG2a.tif"); strcpy(nameG2b, "imG2b.tif"); strcpy(nameG2c, "imG2c.tif"); strcpy(namemag, "immag.tif"); strcpy(name045, "im045.tif"); output (H2a, W, H, idH2a, imH2a, nameH2a, 0); output (H2b, W, H, idH2b, imH2b, nameH2b, 0); output (H2c, W, H, idH2c, imH2c, nameH2c, 0); output (H2d, W, H, idH2d, imH2d, nameH2d, 0); output (G2a, W, H, idG2a, imG2a, nameG2a, 0); output (G2b, W, H, idG2b, imG2b, nameG2b, 0); output (G2c, W, H, idG2c, imG2c, nameG2c, 0); output (t045, W, H, id045, im045, name045, 0); for (j = 0; j < H; j++) { for (i = 0; i < W; i++) { int index1 = i + j * W; int index2 = 3 * (i + j * W); if (fabs(mag[index1]) > drempel) { if ((ori[index1] >= 0) && (ori[index1] <= M_PI / 3)) { imori->rgb_data[index2] = 255 - (ori[index1] / M_PI) * 3 * 255; imori->rgb_data[index2 + 1] = 0.0; imori->rgb_data[index2 + 2] = (ori[index1] / M_PI) * 3 * 255; } if ((ori[index1] >= M_PI / 3) && (ori[index1] <= 2 * M_PI / 3)) { imori->rgb_data[index2] = 0.0; imori->rgb_data[index2 + 1] = ((ori[index1] - M_PI / 3) / M_PI) * 3 * 255; imori->rgb_data[index2 + 2] = 255 - ((ori[index1] - M_PI / 3) / M_PI) * 3 * 255; } if ((ori[index1] >= 2 * M_PI / 3) && (ori[index1] <= M_PI)) { imori->rgb_data[index2] = ((ori[index1] - 2 * M_PI / 3) / M_PI) * 3 * 255; imori->rgb_data[index2 + 1] = 255 - ((ori[index1] - 2 * M_PI / 3) / M_PI) * 3 * 255; imori->rgb_data[index2 + 2] = 0.0; } if ((ori[index1] >= M_PI) && (ori[index1] <= 4 * M_PI / 3)) { imori->rgb_data[index2] = 255 - ((ori[index1] - M_PI) / M_PI) * 3 * 255; imori->rgb_data[index2 + 1] = 0.0; imori->rgb_data[index2 + 2] = ((ori[index1] - M_PI) / M_PI) * 3 * 255; } if ((ori[index1] >= 4 * M_PI / 3) && (ori[index1] <= 5 * M_PI / 3)) { imori->rgb_data[index2] = 0.0; imori->rgb_data[index2 + 1] = ((ori[index1] - 4 * M_PI / 3) / M_PI) * 3 * 255; imori->rgb_data[index2 + 2] = 255 - ((ori[index1] - 4 * M_PI / 3) / M_PI) * 3 * 255; } if ((ori[index1] >= 5 * M_PI / 3) && (ori[index1] <= 2 * M_PI)) { imori->rgb_data[index2] = ((ori[index1] - 5 * M_PI / 3) / M_PI) * 3 * 255; imori->rgb_data[index2 + 1] = 255 - ((ori[index1] - 5 * M_PI / 3) / M_PI) * 3 * 255; imori->rgb_data[index2 + 2] = 0.0; } } else { imori->rgb_data[index2] = imori->rgb_data[index2 + 1] = imori->rgb_data[index2 + 2] = 0; } } } Imlib_save_image (idori, imori, "orient.tif", NULL); output (mag, W, H, idmag, immag, "immag.tif", 1); //output (rec, W, H, idrec, imrec, "imrec.tif", 1); }