Merge pull request #13 from ramanakumars/dev

Fixing docs and adding projection and mosaicing

README.md

@@ -1,6 +1,6 @@
 # Processing and projecting JunoCam images
 
-[<img src="https://readthedocs.org/projects/junocamprojection/badge/?version=latest&style=flat-default">](<LINK>)
+[<img src="https://readthedocs.org/projects/junocamprojection/badge/?version=latest&style=flat-default">](https://junocamprojection.readthedocs.io/en/latest/)
 
 This is a tool to process and project JunoCam images onto a lat/lon grid. 
 

docs/requirements.txt

@@ -13,3 +13,4 @@ sphinx-rtd-theme
 ipykernel
 lxml_html_clean
 myst-parser
+netCDF4

projection/mosaic.py

@@ -1,53 +1,72 @@
 import numpy as np
-import healpy as hp
-from skimage import color
+from scipy.ndimage import gaussian_filter
 import tqdm
 
 
-class Mosaic:
-    def __init__(self, fnames, n_side=512):
-        self.fnames = fnames
-        self.nframes = len(fnames)
-        self.n_side = n_side
+def blend_maps(maps: np.ndarray, blur_sigma: float = 50, trim_threshold: float = 0.25) -> np.ndarray:
+    '''
+    Blend a series of projected maps (all on the same coordinate system) together to reduce seams
 
-    def load_data(self):
-        self.maps = np.zeros((len(self.fnames), hp.nside2npix(self.n_side), 3))
+    :param maps: A numpy array containing each individual JunoCam image projected on the same coordinate system
+    :param blur_sigma: the Gaussian blur radius in pixels for removing large-scale luminance gradient (should be larger than the largest scale feature that needs to be resolved)
+    :param trim_threshold: the standard deviation threshold above which to trim pixels. This is useful for reducing color noise. A large threshold will account for pixels where one color dominates over the others, e.g., when one filter is saturated)
 
-        for i, fname in enumerate(tqdm.tqdm(self.fnames)):
-            map = np.load(fname)
-            if hp.nside2npix(self.n_side) != map.shape[0]:
-                map = hp.ud_grade(map.T, self.n_side).T
+    :return: the blended image in the same size as the original set of images
+    '''
 
-            map_lab = color.rgb2lab(map)
+    # open the file and load the variables
+    nfiles, ny, nx, _ = maps.shape
+    filtered = np.zeros_like(maps, dtype=np.float32)
 
-            # clean up the fringes. these are negative a* and b* in the CIELAB space
-            map_lab[(map_lab[:, 1] < -0.001) | (map_lab[:, 2] < -0.001), 0] = 0
+    for i in tqdm.tqdm(range(nfiles), desc='Applying filter'):
+        mapi = maps[i]
+        Ls = mapi.mean(-1)
 
-            self.maps[i, :] = color.lab2rgb(map_lab)
-        self.maps[~np.isfinite(self.maps)] = 0.
+        # find the image footprint for each map
+        footprint = np.asarray(Ls > np.percentile(Ls, 1), dtype=np.float32)
 
-    def stack(self, radius=2):
-        m_normed = self.maps.copy()
-        # first, norma
-        for j in range(self.nframes):
-            m_normed[j, :] = m_normed[j, :] / np.percentile(m_normed[j, :], 99)
+        # we will blur this so we can trim the edges after we apply the high-pass filter
+        footprint = np.repeat(gaussian_filter(footprint, sigma=blur_sigma, mode=['nearest', 'wrap'])[:, :, np.newaxis], 3, axis=-1)
+        footprint[footprint < 0.8] = 0
+        footprint[footprint >= 0.8] = 1
 
-        count = np.sum(np.min(m_normed, axis=-1) > 1.e-6, axis=0)
-        m_hsv = color.rgb2hsv(m_normed)
-        v_ave = np.mean(m_hsv[:, count > 0, 2])
+        # get the low frequency luminance gradient by applying a Gaussian blur on the luminance
+        low_freq = np.repeat(gaussian_filter(Ls, sigma=blur_sigma, mode=['nearest', 'wrap'])[:, :, np.newaxis], 3, axis=-1)
 
-        pix_inds = np.where(count > 0)[0]
-        v_loc = np.zeros_like(m_normed[:, :, 0])
-        vecs = np.asarray(hp.pix2vec(self.n_side, ipix=pix_inds)).T
+        # the filtered image is divided by the low frequency data (i.e., retains high frequency info)
+        filtered[i] = mapi / (low_freq + 1e-6) * footprint
 
-        for i, vec in tqdm.tqdm(zip(pix_inds, vecs), total=len(pix_inds)):
-            neighbours = hp.query_disc(self.n_side, vec=vec, radius=np.radians(radius))
-            m_hsv_neigh = m_hsv[:, neighbours, 2]
-            v_loc[:, i] = np.mean(m_hsv_neigh[:, count[neighbours] > 0], axis=1)
+        # flatten the image histogram
+        filtered[i] = filtered[i] / np.percentile(filtered[i][filtered[i] > 0], 99)
 
-        m_new = np.zeros_like(m_normed[0, :])
-        for i in range(3):
-            m_new[:, i] = np.sum(m_normed[:, :, i] * (v_ave / v_loc), axis=0) / count
-        m_new[~np.isfinite(m_new)] = 0
+    filtered_Ls = filtered.mean(-1)
 
-        return m_new
+    IMG = np.zeros((ny, nx, 3))
+    alpha = np.zeros_like(filtered[:, :, :, 0])
+    cut_threshold = np.nanpercentile(filtered_Ls[filtered_Ls > 0.0].flatten(), 1)
+
+    for kk in range(nfiles):
+        # find the standard deviation in the image axes, and cut off pixels which are above the
+        # trim threshold for color variance
+        std = np.nanstd(filtered[kk], axis=(-1))
+        mask_k = (filtered_Ls[kk] > cut_threshold) & (std / filtered_Ls[kk] < trim_threshold)
+        # weight each image by its relative brightness wrt to the global mean
+        if len(filtered_Ls[kk, :, :][mask_k]) > 0:
+            alpha[kk][mask_k] = np.nanmean(filtered_Ls[kk, :, :][mask_k])
+            alpha[kk][alpha[kk, :] != 0] = 1. / alpha[kk, :][alpha[kk, :] != 0.]
+
+    alpha_sum = np.sum(alpha, axis=0)
+    alpha_sum[alpha_sum == 0] = np.nan
+
+    IMGs_sub = filtered.copy()
+    for c in range(3):
+        IMGs_sub[:, :, :, c] = IMGs_sub[:, :, :, c] * alpha
+
+    # remove really dim pixels which can occur on the edges and will pull the median down
+    IMGs_sub[IMGs_sub <= np.percentile(filtered_Ls, 40)] = np.nan
+
+    IMG = np.nanmedian(IMGs_sub, axis=0)
+
+    IMG[~np.isfinite(IMG)] = 0.
+
+    return IMG

projection/project.c

@@ -17,6 +17,7 @@ double aperture = 1.;
  */
 double pixel_size = 7.4e-6;
 double focal_length = 0.001095637;
+double c = 3e5; // light speed in km/s
 
 struct Camera
 {
@@ -128,10 +129,10 @@ void furnish(char *file) { furnsh_c(file); }
 void process(double eti, int cam, double cam2jup[3][3], double *lon, double *lat,
              double *inclin, double *emis, double *fluxcal)
 {
-  double pix[2], spoint[3], srfvec[3], pixvec[3], pos_jup[3], scloc[3];
+  double pix[2], spoint[3], srfvec[3], pixvec[3], pos_jup[3], scloc[3], radii[3];
   double disti, loni, lati, phase, inc, emission, trgepoch, lt, plate_scale,
-      ang_size, footprint;
-  int found;
+      ang_size, footprint, f;
+  int found, n;
 
   plate_scale = (1. / focal_length) / (pixel_size);
 
@@ -141,6 +142,7 @@ void process(double eti, int cam, double cam2jup[3][3], double *lon, double *lat
 
   spkpos_c("JUNO", eti, "IAU_JUPITER", "CN+S", "JUPITER", scloc, &lt);
 
+
   for (int jj = 0; jj < FRAME_HEIGHT; jj++)
   {
     for (int ii = 23; ii < FRAME_WIDTH - 17; ii++)
@@ -206,10 +208,10 @@ void project_midplane(double eti, int cam, double tmid, double *lon,
                       double *lat, double *incid, double *emis, double *coords,
                       double *fluxcal)
 {
-  double pix[2], pix_transformed[2], spoint[3], srfvec[3], pixvec[3], scloc[3],
-      vec_transformed[3], vec_iau[3], pxfrm_mid[3][3], pxfrm_iau[3][3], disti,
-      loni, lati, phase, inc, emission, trgepoch, lt, plate_scale, cosalpha;
-  int found;
+  double pix[2], pix_transformed[2], spoint[3], srfvec[3], pixvec[3], scloc[3], radii[3],
+      vec_transformed[3], vec_iau[3], pxfrm_mid[3][3], pxfrm_iau[3][3], alti,
+      loni, lati, phase, inc, emission, trgepoch, lt, plate_scale, cosalpha, f;
+  int found, n;
 
   struct Camera camera, cam0;
   initialize_camera(&camera, cam);
@@ -220,6 +222,10 @@ void project_midplane(double eti, int cam, double tmid, double *lon,
   pxfrm2_c("JUNO_JUNOCAM", "JUNO_JUNOCAM", eti, tmid, pxfrm_mid);
 
   pxform_c("JUNO_JUNOCAM", "IAU_JUPITER", eti, pxfrm_iau);
+
+  bodvrd_c("JUPITER", "RADII", 3, &n, radii);
+
+  f = (radii[0] - radii[2]) / radii[0];
 
   for (int jj = 0; jj < FRAME_HEIGHT; jj++)
   {
@@ -238,7 +244,7 @@ void project_midplane(double eti, int cam, double tmid, double *lon,
       if (found)
       {
         // find the coordinate of the surface point
-        reclat_c(spoint, &disti, &loni, &lati);
+        recpgr_c("JUPITER", spoint, radii[0], f, &loni, &lati, &alti);
 
         // calculate the illumination angle to do
         // solar flux correction
@@ -271,15 +277,20 @@ void project_midplane(double eti, int cam, double tmid, double *lon,
 void get_pixel_from_coords(double *lon, double *lat, int npoints, double et,
                            double *extents, double *pix)
 {
-  double scloc[3], spoint[3], dvec[3], dvec_jcam[3], pixi[2], lt,
-      pxfrm[3][3], x0, x1, y0, y1, percentage;
+  double scloc[3], spoint[3], dvec[3], spoint_et[3], dvec_jcam[3], pixi[2], lt,
+      pxfrm[3][3], pxfrm_et[3][3], x0, x1, y0, y1, percentage, radii[3], f, dist, et2;
+  int n;
   struct Camera cam0;
   initialize_camera(&cam0, 1);
 
   x0 = extents[0];
   x1 = extents[1];
   y0 = extents[2];
   y1 = extents[3];
+
+  bodvrd_c("JUPITER", "RADII", 3, &n, radii);
+
+  f = (radii[0] - radii[2]) / radii[0];
 
   spkpos_c("JUNO", et, "IAU_JUPITER", "CN+S", "JUPITER", scloc, &lt);
   pxform_c("IAU_JUPITER", "JUNO_JUNOCAM", et, pxfrm);
@@ -289,9 +300,8 @@ void get_pixel_from_coords(double *lon, double *lat, int npoints, double et,
     double lati, loni;
     lati = lat[i];
     loni = lon[i];
-    srfrec_c(JUPITER, loni, lati, spoint);
+    pgrrec_c("JUPITER", loni, lati, 0, radii[0], f, spoint);
     subtract3D(spoint, scloc, dvec);
-
     matmul3D(pxfrm, dvec, dvec_jcam);
     vec2pix(&cam0, dvec_jcam, pixi);