/* @(#)sphough.c 17.1.1.1 (ES0-DMD) 01/25/02 17:56:15 */ /*=========================================================================== Copyright (C) 1995 European Southern Observatory (ESO) This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, MA 02139, USA. Corresponding concerning ESO-MIDAS should be addressed as follows: Internet e-mail: midas@eso.org Postal address: European Southern Observatory Data Management Division Karl-Schwarzschild-Strasse 2 D 85748 Garching bei Muenchen GERMANY ===========================================================================*/ /* Program : hough.c */ /* Author : P. Ballester - ESO Garching */ /* Date : 20.11.90 Version 1.0 */ /* 22.03.91 Version 2.0 */ /* 27.03.91 Version 3.0 */ /* */ /* Purpose : */ /* */ /* a) Apply Hough Transform on the central region of an image */ /* */ /* Input: */ /* - name of input image : IN_A */ /* - name of output image : IN_B */ /* - Inter-trace width : INPUTI(1) */ /* - Number of perpendicular traces : INPUTI(2) */ /* - Number of columns of Hough Tr. : INPUTI(3) */ /* - Number of rows of Hough Tr. : INPUTI(4) */ /* - Lower scan limit : INPUTI(5) */ /* - Upper scan limit : INPUTI(6) */ /* - Start of Hough Transform : INPUTD(1),(2) */ /* - Step of Hough Transform : INPUTD(3),(4) */ /* */ /* Limits: */ /* */ /* a) This algorithm must process only positive values. */ /* b) The row 0 of the transformed space is not computed */ /* */ /* Algorithm: * * a) Only a certain number of columns of the image are processed. * Given the distance between two columns and the number of columns * to process (INPUTI(1) and (2)), the positions of the different * columns are given by the function icol(ncol,trace,ntrace,inter). * * b) The image is first filtered and a constant value (background) * is subtracted. The filtering consists of a median estimate over a * kernel 3*5 pixels. The constant value corresponds to the minimum * of all median estimates. The result is stored in a buffer of size * ntrace * nrow (number of columns processed, number of rows of the * original image). The first and last row (1 and n) of the filtered * image are copies of the adjacent rows (2 and n-1). * * c) The Hough Transform is performed, in which the incremented value * corresponds to the value of the pixel in the processed image. For * more details on the Hough transform, see P. Ballester, 1991. Proceedings * of the 3rd Data Analysis Workshop and "Hough transform for robust * regression and automated detection" in Astronomy and Astrophysics * (March 1994). */ #include #include #include #include #include int inull; /* Representation of NULL values */ float rnull; double dnull; int next_prgs = 10, step_prgs = 10; main () { char catal[81], transf[81], line[81], text[81]; char colx[11], colw[11], mode[11], idmode[11]; int npix_hg[3], ndim, dim, kunit; double start_hg[3], step_hg[3]; int null, status, npar, iav, actvals; int inpi[2], inter, ntrace, scan[2]; int linx, catwav; int nbcol, nbrow, nsort, allcol, allrow; int nbcolcat, nbrowcat, nsortcat, allcolcat, allrowcat; int tidlin, tidcat, rowlin, rowcat; double x, wave, peak_coord[3], fct(); float unity = 1., rpar[5], avdisp, wcenter, delta; float range, xnorm, xc, findmax(), max; int col1, col2, col3, linid, lindif, lininf, linwav; double *colpos, *colcat, *colwav, *colid, *coldif, *colinf, absdif, mindif; int *select, ncat, nbsel, row; SCSPRO ("hough"); SCKRDI ("INPUTI", 4, 1, &actvals, &ndim, &kunit, &null); SCKRDI ("INPUTI", 1, 3, &actvals, npix_hg, &kunit, &null); SCKRDD ("INPUTD", 1, 3, &actvals, start_hg, &kunit, &null); SCKRDD ("INPUTD", 4, 3, &actvals, step_hg, &kunit, &null); SCKRDR ("INPUTR", 1, 4, &actvals, rpar, &kunit, &null); wcenter = rpar[0], avdisp = rpar[1], xc = rpar[2], range = rpar[3]; SCKGETC ("IN_A", 1, 60, &actvals, line); SCKGETC ("IN_B", 1, 60, &actvals, catal); SCKGETC ("OUT_A", 1, 60, &actvals, transf); SCKGETC ("OUT_B", 1, 10, &actvals, mode); SCKGETC ("OUT_B",10, 10, &actvals, idmode); SCKGETC ("INPUTC", 1, 10, &actvals, colx); SCKGETC ("INPUTC", 10, 10, &actvals, colw); TCTOPN(line, F_IO_MODE, &tidlin); TCIGET(tidlin, &nbcol, &nbrow, &nsort, &allcol, &allrow); TCCSER(tidlin, colx, &linx); if (catal[0] != '@') { TCTOPN(catal, F_I_MODE, &tidcat); TCIGET(tidcat, &nbcolcat, &nbrowcat, &nsortcat, &allcolcat, &allrowcat); TCCSER(tidcat, colw, &catwav); } else { nbrowcat = nbrow; tidcat = tidlin; TCCSER(tidlin, colw, &catwav); } TCCSER(tidlin, ":WAVEC", &linwav); if (linwav == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom", "WAVEC ", &linwav); colwav = (double *) osmmget((nbrow+1)*sizeof(double)); select = (int *) osmmget((nbrow+1)*sizeof(int)); colcat = (double *) osmmget((nbrowcat+1)*sizeof(double)); colpos = (double *) osmmget((nbrow+1)*sizeof(double)); ncat = read_col(tidcat, nbrowcat, catwav, colcat); sprintf(text,"Number of lines in catalog (total, selected) : %d , %d", nbrowcat, ncat), SCTPUT(text); nbsel = read_select(tidlin, nbrow, select); sprintf(text,"Number of lines in table line (total, selected) : %d , %d", nbrow, nbsel), SCTPUT(text); read_col(tidlin, nbrow, linx, colpos); nbrow = nbsel, nbrowcat = ncat; /* Generate the Hough transform accumulator */ create_hough(transf, npix_hg, start_hg, step_hg, ndim); for (rowlin=1; rowlin<=nbrow; rowlin++) { display_progress(rowlin,nbrow); x = colpos[rowlin]; if (catal[0] != '@') { for (rowcat=1; rowcat<=nbrowcat; rowcat++) { wave = colcat[rowcat]; increment_hough(x,wave,unity,mode,avdisp,range); } } else { wave = colcat[rowlin]; increment_hough(x,wave,unity,mode,avdisp,range); } } /* Explore the accumulator and determine the location of the maximum */ max = findmax(&col1,&col2,&col3); sprintf (text,"Found maximum %f at location: x=%d y=%d z=%d\n", max,col1,col2,col3); SCTPUT(text); /* for (dim = 3; dim > ndim; dim--) {start_hg[dim] = 0.; step_hg[dim] = 0.;} */ peak_coord[0] = start_hg[0] + (col1-1)*step_hg[0]; peak_coord[1] = start_hg[1] + (col2-1)*step_hg[1]; peak_coord[2] = start_hg[2] + (col3-1)*step_hg[2]; sprintf(text,"Coordinates: %f %f %f\n",peak_coord[0], peak_coord[1],peak_coord[2]); SCTPUT(text); SCKWRI("OUTPUTI",&col1,3,1,&kunit); SCKWRI("OUTPUTI",&col2,4,1,&kunit); SCKWRI("OUTPUTI",&col3,5,1,&kunit); SCKWRR("OUTPUTR",&max,1,1,&kunit); SCKWRD("OUTPUTD",peak_coord,1,3,&kunit); close_hough(); /* Write Computed wavelength in the :WAVEC column */ for (rowlin=1; rowlin<=nbrow; rowlin++) { x = colpos[rowlin]; switch(mode[0]) { case('1'): colwav[rowlin] = peak_coord[0] + avdisp*x; delta = 0.; break; case('L'): colwav[rowlin] = peak_coord[1] + peak_coord[0]*x; delta = 2.*peak_coord[0]*range; break; case('N'): colwav[rowlin] = peak_coord[1] + avdisp*x*(1.+peak_coord[0]*x); delta = 2.*avdisp*range*(1.+2*peak_coord[0]*x); break; case('3'): colwav[rowlin] = peak_coord[1] + x*peak_coord[0]*(1.+peak_coord[2]*x); delta = 2.*peak_coord[0]*range*(1.+2*peak_coord[2]*x); break; } /* printf("rowlin, select, wavec: %d %d %f\n",rowlin,select[rowlin], colwav[rowlin]); */ } write_dcol(tidlin, nbrow, select, linwav, colwav); /* Performs identification if idmode="IDENT" */ if (toupper(idmode[0]) == 'I') { /* Allocate arrays and read position and catalog values */ TCCSER(tidlin, ":IDENT", &linid ); if (linid == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom", "IDENT ", &linid ); TCCSER(tidlin, ":RESIDUAL", &lindif); if (lindif == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom", "RESIDUAL", &lindif); TCCSER(tidlin, ":INFLUENCE", &lininf); if (lininf == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom", "INFLUENCE", &lininf); TCMNUL(&inull, &rnull, &dnull); colid = (double *) osmmget((nbrow+1)*sizeof(double)); coldif = (double *) osmmget((nbrow+1)*sizeof(double)); colinf = (double *) osmmget((nbrow+1)*sizeof(double)); /* Write null values in all output columns */ for (row=0; row<=nbrow; row++) { colwav[row] = dnull; colid[row] = dnull; coldif[row] = dnull; colinf[row] = dnull; } write_dcol(tidlin, nbrow, select, linid, colid); write_dcol(tidlin, nbrow, select, lindif, coldif); write_dcol(tidlin, nbrow, select, lininf, colinf); /* Identifies the lines by association */ if (range < 0.) delta = (double) range; if (delta < 0.) delta *= -1.; delta /= step_hg[1]; if (catal[0] == '@') nbrowcat = nbrow; for (rowlin=1; rowlin<=nbrow; rowlin++) { for (rowcat=1; rowcat<=nbrowcat; rowcat++) { wave = colcat[rowcat]; coldif[rowlin] = colwav[rowlin] - wave; if (coldif[rowlin] < 0.) absdif = (-1.)*coldif[rowlin]; else absdif = coldif[rowlin]; if (rowcat == 1) mindif = absdif; if (absdif <= mindif) { colid[rowlin] = wave; mindif = absdif; } /* printf( "rowlin, rowcat, absdif, wave, ident, computed: %d, %d %f %f %f %f\n", rowlin, rowcat , absdif, wave, colid[rowlin], colwav[rowlin]); */ coldif[rowlin] = colwav[rowlin] - colid[rowlin]; colinf[rowlin] = fct(coldif[rowlin],delta); /* printf("wavec, ident, diff, inf: %f %f %f %f\n", colwav[rowlin],colid[rowlin],coldif[rowlin],colinf[rowlin]); */ }} /* Update and close all tables */ write_dcol(tidlin, nbrow, select, linid,colid); write_dcol(tidlin, nbrow, select, lindif, coldif); write_dcol(tidlin, nbrow, select, lininf,colinf); } /* Matches if (toupper(idmode[0] == 'I') .. */ if (catal[0] != '@') TCTCLO(tidcat); TCTCLO(tidlin), SCSEPI(); } read_col(tid, nb, colnb, col) int tid, colnb, nb; double col[]; { int i, k=0, null, sel; for (i=1; i<=nb; i++) { TCSGET(tid, i, &sel); if (sel) { TCERDD(tid, i, colnb, &col[++k], &null); if (null) col[k] = dnull; } } return(k); } read_select(tid, nb, col) int tid, nb, col[]; { int i, k=0, null, sel; for (i=1; i<=nb; i++) { TCSGET(tid, i, &sel); if (sel) col[++k] = i; } return(k); } write_dcol(tid, nb, select, colnb, col) int tid, nb, colnb, select[]; double col[]; { int i; for (i=1; i<=nb; i++) TCEWRD(tid, select[i], colnb, &col[i]); } display_progress(row, nrow) int row, nrow; { float prgs; extern int next_prgs, step_prgs; char date[28], text[80]; struct tm date_struct; prgs = 100. * (float) row / (float) nrow; if (prgs > next_prgs) { if (oshdate(date,&date_struct)) date[0] = '\0'; sprintf(text,"%s %d %% performed...", date, next_prgs); next_prgs += step_prgs; SCTPUT(text); } }