include # CR_FIND -- Find cosmic ray candidates. # This procedure is an interface to special procedures specific to a given # window size. procedure cr_find (cr, threshold, data, nc, nl, col, line, sf1, sf2, x, y, z, w) pointer cr # Cosmic ray list real threshold # Detection threshold pointer data[ARB] # Data lines int nc # Number of columns int nl # Number of lines int col # First column int line # Center line pointer sf1, sf2 # Surface fitting real x[ARB], y[ARB], z[ARB], w[ARB] # Surface arrays pointer a, b, c, d, e, f, g begin switch (nl) { case 5: a = data[1] b = data[2] c = data[3] d = data[4] e = data[5] call cr_find5 (cr, threshold, col, line, Memr[a], Memr[b], Memr[c], Memr[d], Memr[e], nc, sf1, sf2, x, y, z, w) case 7: a = data[1] b = data[2] c = data[3] d = data[4] e = data[5] f = data[6] g = data[7] call cr_find7 (cr, threshold, col, line, Memr[a], Memr[b], Memr[c], Memr[d], Memr[e], Memr[f], Memr[g], nc, sf1, sf2, x, y, z, w) } end # CR_FIND7 -- Find cosmic rays candidates in 7x7 window. # This routine finds cosmic rays candidates with the following algorithm. # 1. If the pixel is not a local maximum relative to it's 48 neighbors # go on to the next pixel. # 2. Identify the next strongest pixel in the 7x7 region. # This suspect pixel is excluded in the following. # 2. Compute the flux of the 7x7 region excluding the cosmic ray # candidate and the suspect pixel. # 3. The candidate must exceed the average flux per pixel by a specified # threshold. If not go on to the next pixel. # 4. Fit a plane to the border pixels (excluding the suspect pixel). # 5. Subtract the background defined by the plane. # 6. Determine a replacement value as the average of the four adjacent # pixels (excluding the suspect pixels). # 7. Add the pixel to the cosmic ray candidate list. procedure cr_find7 (cr, threshold, col, line, a, b, c, d, e, f, g, n, sf1, sf2, x, y, z, w) pointer cr # Cosmic ray list real threshold # Detection threshold int col # First column int line # Line real a[ARB], b[ARB], c[ARB], d[ARB] # Image lines real e[ARB], f[ARB], g[ARB] # Image lines int n # Number of columns pointer sf1, sf2 # Surface fitting real x[49], y[49], z[49], w[49] # Surface arrays real bkgd[49] int i1, i2, i3, i4, i5, i6, i7, j, j1, j2 real p, flux, replace, asumr() pointer sf begin for (i4=4; i4<=n-3; i4=i4+1) { # Must be local maxima. p = d[i4] if (p z[j2]) j2 = j } # Compute the flux excluding the extreme points. flux = (asumr (z, 49) - z[j1] - z[j2]) / 47 # Pixel must be exceed specified threshold. if (p < flux + threshold) next # Fit and subtract the background. if (j2 < 25) { w[j2] = 0 sf = sf2 call gsfit (sf, x, y, z, w, 24, WTS_USER, j) w[j2] = 1 } else { sf = sf1 call gsrefit (sf, x, y, z, w, j) } call gsvector (sf, x, y, bkgd, 49) call asubr (z, bkgd, z, 49) p = z[j1] # Compute the flux excluding the extreme points. flux = (asumr (z, 49) - z[j1] - z[j2]) / 47 # Determine replacement value from four nearest neighbors again # excluding the most deviant pixels. replace = 0 j = 0 if (j2 != 32) { replace = replace + c[i4] j = j + 1 } if (j2 != 36) { replace = replace + d[i3] j = j + 1 } if (j2 != 38) { replace = replace + d[i5] j = j + 1 } if (j2 != 42) { replace = replace + e[i4] j = j + 1 } replace = replace / j # Add pixel to cosmic ray list. flux = 100. * flux call cr_add (cr, col+i4-1, line, flux, flux/p, 0., replace, 0) i4 = i7 } end # CR_FIND5 -- Find cosmic rays candidates in 5x5 window. # This routine finds cosmic rays candidates with the following algorithm. # 1. If the pixel is not a local maximum relative to it's 24 neighbors # go on to the next pixel. # 2. Identify the next strongest pixel in the 5x5 region. # This suspect pixel is excluded in the following. # 2. Compute the flux of the 5x5 region excluding the cosmic ray # candidate and the suspect pixel. # 3. The candidate must exceed the average flux per pixel by a specified # threshold. If not go on to the next pixel. # 4. Fit a plane to the border pixels (excluding the suspect pixel). # 5. Subtract the background defined by the plane. # 6. Determine a replacement value as the average of the four adjacent # pixels (excluding the suspect pixels). # 7. Add the pixel to the cosmic ray candidate list. procedure cr_find5 (cr, threshold, col, line, a, b, c, d, e, n, sf1, sf2, x, y, z, w) pointer cr # Cosmic ray list real threshold # Detection threshold int col # First column int line # Line real a[ARB], b[ARB], c[ARB], d[ARB], e[ARB] # Image lines int n # Number of columns pointer sf1, sf2 # Surface fitting real x[25], y[25], z[25], w[25] # Surface arrays real bkgd[25] int i1, i2, i3, i4, i5, j, j1, j2 real p, flux, replace, asumr() pointer sf begin for (i3=3; i3<=n-2; i3=i3+1) { # Must be local maxima. p = c[i3] if (p z[j2]) j2 = j } # Compute the flux excluding the extreme points. flux = (asumr (z, 25) - z[j1] - z[j2]) / 23 # Pixel must be exceed specified threshold. if (p < flux + threshold) next # Fit and subtract the background. if (j2 < 17) { w[j2] = 0 sf = sf2 call gsfit (sf, x, y, z, w, 16, WTS_USER, j) w[j2] = 1 } else { sf = sf1 call gsrefit (sf, x, y, z, w, j) } call gsvector (sf, x, y, bkgd, 25) call asubr (z, bkgd, z, 25) p = z[j1] # Compute the flux excluding the extreme points. flux = (asumr (z, 25) - z[j1] - z[j2]) / 23 # Determine replacement value from four nearest neighbors again # excluding the most deviant pixels. replace = 0 j = 0 if (j2 != 18) { replace = replace + b[i3] j = j + 1 } if (j2 != 20) { replace = replace + c[i2] j = j + 1 } if (j2 != 22) { replace = replace + c[i4] j = j + 1 } if (j2 != 24) { replace = replace + d[i3] j = j + 1 } replace = replace / j # Add pixel to cosmic ray list. flux = 100. * flux call cr_add (cr, col+i3-1, line, flux, flux/p, 0., replace, 0) i3 = i5 } end