streaks_remove(inImage,                 // Input image
               inMask,                  // Input mask (optional)
               inWeight,                // Input weight map (optional)
               inAstrom,                // Input astrometry (optional)
               maskValue,               // Mask value to set
               recoveryRoot,            // Root name for recovery files
               streaks,                 // Streaks table
               warps                    // Warps table
               )
{
    // Lock Nebulous files to prevent access while changing in-place
    neb_lock(inImage);
    if (inMask) neb_lock(inMask);
    if (inWeight) neb_lock(inWeight);

    // Optionally load astrometry
    if (inAstrom) astrom = astrometry_read(inAstrom);

    // Iterate through components of image
    foreach extension(inImage) {
        // Read header and extract astrometry
        header = header_read(inImage);
        if (!inAstrom) astrom = astrometry_from_header(header);

        // Identify pixels to mask because of streaks
        streak_pixels = pixel_list_initialise(NONE); // List of pixels to be excised because of streak
        foreach streak (streaks) {
            // Calculate streak parameters on image
            image_streak = streak_on_component(streak, astrom, num_cols, num_rows);
            // Only bother with streaks that touch the image
            next if image_streak not defined;

            foreach pixel specified by local_streak {
                add_pixel(streak_pixels, pixel);
            }
        }

        // Identify pixels to mask because not in a warp
        warp_pixels = pixel_list_initialise(ALL); // List of pixels to be excised because not in a warp
        foreach warp (warps) {
            // Calculate warp region on image
            image_warp = warp_on_component(warp, astrom, num_cols, num_rows);
            foreach pixel specified by image_warp {
                remove_pixel(warp_pixels, pixel);
            }
        }

        // Mark the header as "destreaked"
        mark_header(inImage);
        if (inMask) mark_header(inImage);
        if (inWeight) mark_header(inImage);

        if (size(streak_pixels) > 0 || size(warp_pixels) > 0) {
            // Read the pixels
            if (num_streaks == 0) {
                image = image_read(inImage);
                if (inMask) mask = image_read(inMask);
                if (inWeight) weight = image_read(inWeight);

                recovery_image = initialise_image(0);
                if (inMask) recovery_mask = initialise_image(0);
                if (inWeight) recovery_weight = initialise_image(0);
            }

            // Mask the pixels
            foreach masked_pixel (streak_pixels, warp_pixels) {
                recovery_image[x,y] = image[x,y];
                image[x,y] = NAN;
                if (inMask) {
                    recovery_mask[x,y] = mask[x,y];
                    mask[x,y] = maskValue;
                }
                if (inWeight) {
                    recovery_weight[x,y] = weight[x,y];
                    weight[x,y] = NAN;
                }
            }

            // Save the outputs
            image_save(image);
            image_write_compressed(recoveryRoot + ".im.fits", recovery_image);

            if (inMask) {
                image_save(mask);
                image_write_compressed(recoveryRoot + ".mk.fits", recovery_mask);
            }
            if (inWeight) {
                image_save(mask);
                image_write_compressed(recoveryRoot + ".wt.fits", recovery_weight);
            }
        }
    }

    // Unlock Nebulous
    neb_unlock(inImage);
    if (inMask) neb_unlock(inMask);
    if (inWeight) neb_unlock(inWeight);

    return SUCCESS;
}
