streaks_remove(inImage,                 // Input image  ?: char*?
               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    ?: table of RA/Dec endpoints?
               warps                    // Warps table ?: what is this format?
               )
{
    // Lock Nebulous files to prevent access while changing in-place
    neb_lock(inImage);
    if (inMask) neb_lock(inMask);   // ?: Is the call format, neb_lock (const char* path)?
    if (inWeight) neb_lock(inWeight);

    // Optionally load astrometry
    if (inAstrom) astrom = astrometry_read(inAstrom);   
    // ?: Is this a UH provided routine?  EAM : yes
    // ?: What is the structure for astrom? EAM : TBD
    // ?: Definition of higher order WCS keywords? : EAM : added to ICD
    
    // ?: For detrend and raw files, are the higher order terms stored in
    //    a separate file?  EAM : yes

    // Is the unique higher order WCS per extension or is it based on an image. EAM : see ICD

    // Iterate through components of image
    foreach extension(inImage) {
        // Read header and extract astrometry
        header = header_read(inImage);
        if (!inAstrom) astrom = astrometry_from_header(header); 
	// ?: Is this an affine transform only? 

	// EAM : no.  the optical distortion is third order.  However, the non-linear terms are
	// perturbations on the linear terms, and assuming there are no caustics in the optics
	// (should not be!), there is a one-to-one relationship between input and output
	// coordinates.  This gives us two simplifications: 1) we can use a linear approximation
	// for localized patches (we use this as an optimization in the warping, and I recommed
	// we use it here as well).  2) the inverse transformation can be closely represented
	// by a polynomial of the same order.  Or, to be more accurate, an approximate inverse
	// transformation of the same order only makes small errors.  in our analysis, we
	// define one direction (chip->focal plane->tangent plane->sky) to be the primary
	// transformation, and the inverse is understood to be lower accuracy and used only for
	// finding relationships. (ie, match reference star to detection).

        // 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(inMask);   
        if (inWeight) mark_header(inWeight); 

	// ?: mark_header(inMask);  // Record masking in main image or mask?
	// ?: mark_header(inWeight);// Record masking in weight image or mask?

	// EAM : these operations are mean to describe the addition of metadata to the headers
	// to record that magic has been applied (and any interesting stats, like number of
	// pixels removed -- NOTE: that means this should come after the following block).  we
	// would want to record the metadata in each of image, weight, and mask so that they
	// are all self-documenting.

        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;   
		// ?: Really record it as Not A Number (NAN)?
		// EAM : that is our preference: we use NAN to mark pixels with invalid data elsewhere

                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;
                }
            }

            // ?: No temporary files?

	    // EAM : temporary files are probably required, just ignored in this description.
	    // We are assuming that the process loops over the cells (extensions) in the image
	    // to minimize memory footprint, so we can't write the result to the same output
	    // file. 

            // Save the outputs (write for current component/image extension)
            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);
            }

            // ?: Add writing of streak_pixels and warp_pixels list for image?

	    // EAM : Paul Price has made an alternate proposal for the storage of the 'recovery
	    // image', which we used in defining this: make that image consist of only the
	    // masked pixels, with everything else set to 0.0.  His point is that such an
	    // image, if written RICE compressed, would compress extremely well.  It would also
	    // make recovery extremely simple: just read it in the two image and replace into
	    // the 'magicked image' any non-zero pixels from the recover image.

        }           // if (size(streak_pixels) > 0 || size(warp_pixels) > 0)
    }               // foreach extension(inImage)

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

    return SUCCESS;
}
