Veil Nebula, about 4 hours

Veil Nebula, about 4 hours

This post is outdated, this is the new article.

A dense starfield can be the amateur astronomer’s nightmare for it can effectively hide the faint nebula behind it. The more errors add up, the less is reveiled of the target. Why? Because the stars that are supposed to be pointlike get big, lens errors blow up the stars, bad seeing blows up the stars, any amount of dew on the lens or fog spoil the stars while the nebula is still there. So I thought of something: what if I could tune the volume of the stars. Another train of thought arose from an image artifact I keep getting – a panda eye like feature, ie a dark ring around the pointlike luminous stars. These two problems amazingly have almost the same solution: throw out the stars and use the remaining context. The point being aesthetics and not science, I started to brew a software.

The program

My algorithm does not produce a final picture, it just removes the stars and much work is still needed after it does its job. In layman’s terms it works like this:

0) since my program is 8 bits for now, it receives an input image that has almost been finished from other points of view
1) the program receives an input.jpg file saved with 100% quality
2) the program identifies pixels above an L luminosity threshold
3) it replaces the identified pixels and their R surroundings with black (that’s 0x000000 in RGB)
4) the program identifies pixels that are grayish (very low saturation) and more luminous than part of their surroundings by at least D difference
5) it replaces the identified pixels and their R surroundings with black
6) the program iterates taking the black pixels and if they have non-black  surroundings, replaces them with the average of the  surroundings, thus progressing from the margin of the black areas until the wounds are healed
7) there are two outputs: a) the healed picture with no stars and b) the “lighter than” (a max function) combination of the original and the one from point a).
8) the needed output goes back into the photo editing software as a layer

4) and 5) are meant to remove fainter stars. I don’t always use them.

Examples – before-after

While the difference is not always obvious on the downscaled images, I provide some zoom ins to:

input

input

carved

carved

healed

healed

before

before

after

after

input

input

carved

carved

healed

healed

before

before

after

after

before

before

after

after

before

before

after

after

input

input

carved

carved

healed

healed

before

before

after

after

before

before

after

after

Examples: final pictures

Caldwell 49 Rosette Nebula (20×150 sec, modified Canon 1100D @ Canon EF 200mm f/2.8L II USM at f/4, EQ3)

Caldwell 49 Rosette Nebula (20×150 sec, modified Canon 1100D @ Canon EF 200mm f/2.8L II USM at f/4, EQ3)

NGC 7000 North America Nebula

NGC 7000 North America Nebula

The Milky Way in Chepheus and Cygnus – 50×3 min at 50mm

The Milky Way in Chepheus and Cygnus – 50×3 min at 50mm

 

The source code

<!-- 
script to fill the blanks of input.jpg 
by extrapolating the context of the blank, 
saving the output as output*.jpg
-->
<div align="center">
<?php

report("=============starting the script=============");


$i = imagecreatefromstring(file_get_contents('input.jpg'));

report("created the input image");

//settings
$GLOBALS['minlumi'] = 5;
$GLOBALS['maxlumi'] = 140;
$GLOBALS['max_damage'] = 250;
$GLOBALS['csillaglumi'] = 150;
$GLOBALS["kill_highlight_radius"] = 2;  //2 or 3 

$stages_to_save = 25;

set_time_limit(60000);//majdnem 24 ora :P
ini_set('memory_limit', '2048M'); //2 giga, szoval legyszives
ignore_user_abort(true);

//========================================================================
function report($s){
   $s = $s.'<br />';
   $s = str_replace("<br /><br />", "<br />", $s);
   $s = date('Y.m.d. H:i:s ').$s;
   echo $s;
   $s = str_replace("<br />", "\r\n", $s)."\r\n";
   $s = str_replace("\r\n\r\n", "\r\n", $s);
   file_put_contents('log.txt', $s, FILE_APPEND | LOCK_EX);
}
//========================================================================
function numberaskilo($i){
   while (strlen($i) < 3){
      $i = "0".$i;
   }
   return $i;
}
//========================================================================
function int2rgb($i){
   if (is_array($i)){
      return $i;
   }
   $r = $i >> 16;
   $g = ($i >> 8) %256;
   $b = $i % 256;
   return  array($r, $g, $b);
}
//========================================================================
function rgb2int($r, $g=false, $b=false){
   if (is_array($r)){
      $b = $r[2];
      $g = $r[1];
      $r = $r[0];
   }
   return $r * 256*256 + $g*256 + $b;
}
//========================================================================
function luminanceof($i_or_rgb){
   if (!is_array($i_or_rgb)){
      $i_or_rgb = int2rgb($i_or_rgb);
   }
   return floor(($i_or_rgb[0] + $i_or_rgb[1]*2 + $i_or_rgb[2]) / 4); 
}
//========================================================================
function imagesetpixel_($i, $tx, $ty, $p){
   if (($tx >= 0) && ($tx < imagesx($i)) && ($ty >= 0) && ($ty < imagesy($i))){
      imagesetpixel($i, $tx, $ty, $p);
   };
}
//========================================================================
function imagecolorat_($i, $tx, $ty){
   if (($tx >= 0) && ($tx < imagesx($i)) && ($ty >= 0) && ($ty < imagesy($i))){
      return imagecolorat($i, $tx, $ty);
   };
   return 0;
}
//========================================================================
function maxcolororluminance($pix){
   return luminanceof($pix);
   ////
   $pix  = int2rgb($pix);
   $pix[] = floor(($pix[0] + $pix[1] + $pix[1] + $pix[2]) / 4); 
   return max($pix);
}
//========================================================================
function heal_pixel_if_possible($i, $x, $y){
       $xy = $x.'_'.$y;
       $pix = imagecolorat($i, $x, $y);
       if (maxcolororluminance($pix) >= $GLOBALS['minlumi']){
           return false;
       };
       $u = array();
       for ($tx = -1; $tx < 2; $tx++){
          for ($ty = -1; $ty < 2; $ty++){
             if (($tx==0)&&($ty==0)){
                // _self
             }else{
                $p = imagecolorat_($i, $tx+$x, $ty+$y);
                $l = maxcolororluminance($p);
                if (($l < $GLOBALS['maxlumi']) && ($l >= $GLOBALS['minlumi'])){
                   $u[] = int2rgb($p);
                }
             };   
          }
       }
       
       if (count($u) >= $GLOBALS['mincontext']){
          $uj = array(0,0,0);
          foreach ($u as $j){
             for ($q=0; $q<3; $q++){
                $uj[$q] = $uj[$q] + $j[$q];
             };                
          }
          for ($q=0; $q<3; $q++){
             $uj[$q] = floor($uj[$q] / count($u));
          };         
          imagesetpixel($i, $x, $y, rgb2int($uj));          
          return true;
       }
}
//========================================================================
function handle_this_pixel(&$i, &$x, &$y, &$healed, &$save_stage_at_x_heals, &$st){
   if (!file_exists("run.txt")){
      imagejpeg($i, "output_8_premature.jpg", 95);
      die();
   }
   if (heal_pixel_if_possible($i, $x, $y)){
      $healed++;
      if ($healed % $save_stage_at_x_heals == 0){
          imagejpeg($i, 'output_3_stage_'.numberaskilo($st).'.jpg', 95);
          $st++;
      }
   };
};
//========================================================================
function kill_starmap(&$i, &$star_map, $radius = -1){
   if ($radius == -1){   
      $radius = $GLOBALS["kill_highlight_radius"];
   };      
   foreach ($star_map as $star){
      $min = -$radius;
      $max = $radius;
      for ($rx = $min; $rx <= $max; $rx++){
         for ($ry = $min; $ry <= $max; $ry++){                   
            if (
                (($rx == $min) OR ($rx == $max)) 
                AND 
                (($ry == $min) OR ($ry == $max)) 
               ){
               //this is a corner, but lets make the hole a bit round
            }else{
               imagesetpixel_($i, $star[0]+$rx, $star[1]+$ry, 0);
            }
         }
      }
   }
}
//========================================================================
function kill_off_the_highlights(&$i){
   report("killing off the highlights");
   $star_map = array();
   for ($x = 0; $x < imagesx($i); $x++){
      for ($y = 0; $y < imagesy($i); $y++){
          $pix = imagecolorat($i, $x, $y);
          if (maxcolororluminance($pix) > $GLOBALS['csillaglumi']){
             $star_map[] = array($x, $y);
          };
      }
   }
   kill_starmap($i, $star_map);
   
   imagejpeg($i, "input_2_nohighlights.jpg", 100);
   report('done killing off the highlights');
};
//========================================================================
function kill_off_remaining_starlike_areas(&$i){
   ////ha minden iranyba ugyanaz vagy kisebb, akkor jelolt
   $lumidelta = 25;
   $max_coldelta = 45;
   $radius = 3;
   $kill = array();
   $oszto = 10;
   $progress = round(imagesx($i) / $oszto);
   $ko = 0;
   for ($x =0; $x < imagesx($i); $x++){
      if (($x+1) % $progress == 0){
           $ko++;
           file_put_contents("input_4_no_smallerstars_".numberaskilo($ko)."_of_".numberaskilo($oszto).".txt", 100);              
      }
      if (!file_exists('run.txt')){
         die("premature exit");
      }
      for ($y = 0; $y < imagesy($i); $y++){
         $p = imagecolorat($i, $x, $y);
         if ($p > 0){
             $rgb = int2rgb($p);
             $atlag = floor(($rgb[0] + $rgb[1]*2 + $rgb[2])/4);
             $coldelta = 0;
             for ($j =0; $j<=2; $j++){
                $coldelta = max($coldelta, abs($rgb[$j] - $atlag));
             }
             //based on hue, this could be a star, check its context
             if ($coldelta < $max_coldelta){
             
                $context = array();
                $context_is_lighter = false;
                for ($tx = -$radius; $tx <= $radius; $tx++){
                   for ($ty = -$radius; $ty <= $radius; $ty++){
                      if ((abs($tx)==$radius)AND(abs($ty)==$radius)){
                         //round the corners
                      }else{
                         $k = imagecolorat_($i, $x+$tx, $y+$ty);
                         $k_rgb = int2rgb($k);
                         $k_atlag = floor(($k_rgb[0] + $k_rgb[1]*2 + $k_rgb[2])/4);
                         if ($k_atlag > $atlag){
                            $context_is_lighter = true;
                         }
                         $context[] = $k;
                      }                  
                   }
                }
                           
                $fainter = 0;
                $black = 0;
                if (!$context_is_lighter){
                    $context_looks_good = 0;                 
                    foreach ($context as $c) if (($c!=0)OR(true)) {
                        if ($c ==0 ){
                           $black++;
                        } 
                        $c_rgb = int2rgb($c);
                        $c_atlag = floor(($c_rgb[0] + $c_rgb[1]*2 + $c_rgb[2])/4);
                        $c_coldelta = 0;
                        for ($j =0; $j<=2; $j++){
                           $c_coldelta = max($c_coldelta, abs($c_rgb[$j] - $c_atlag));
                        }
                        
                        $is_fainter = ($c_atlag < $atlag-$lumidelta);
                        if ($is_fainter){
                           $fainter++;
                        }
                        if (((abs($c_coldelta - $coldelta) < 15) AND (abs($c_atlag - $atlag) < 5)) OR ($is_fainter)){
                           $context_looks_good++;
                        }
                    }
                     
                    if ((abs(count($context) - $context_looks_good) < max(2, count($context) / 4)) AND ($fainter > 0) AND ($black < count($context) / 3)){
                        $kill[] = array($x, $y);
                    }
                };
             }//pixel is grey enough
         }//pixel not black    
      }
   }
   kill_starmap($i, $kill, $radius);
   imagejpeg($i, "input_6_no_smallerstars.jpg", 100);   
}
//========================================================================
function try_loading_a_cache_if_exists(&$i){
    $icv = 'input_9_carved.jpg';
    
    //remove this cache, start from scratch
    if (file_exists($icv)){
       unlink($icv);
    }
    
    if (!file_exists($icv)){
       report('image input loaded, no cache<br />');
       
       
       kill_off_the_highlights($i);
       
       // kill_off_remaining_starlike_areas($i);

       report('saving the carved version<br />');
       imagejpeg($i, $icv, 100);
       
       foreach (glob("*.jpg") as $filename){
          echo '<img src="'.$filename.'?r='.mt_rand().'" />';
       }       
       //die();
    }else{
       report('image outpIt loaded from cache<br />');
       imagedestroy($i);
       $i = imagecreatefromstring(file_get_contents($icv));
    }
}
function filter_lighterthanoriginal(&$i){
    //apply a "lighter than" filter, assuming the "to be replaced" areas were black
    $o = imagecreatefromstring(file_get_contents("input.jpg"));
    for ($x = 0; $x< imagesx($o); $x++){
       for ($y =0; $y< imagesy($o); $y++){
          $pixo  = imagecolorat($o, $x, $y);
          $pixi  = imagecolorat($i, $x, $y);
          if (luminanceof($pixi) > luminanceof($pixo)){
            imagesetpixel($o, $x, $y, $pixi);
          } 
       };      
    };
    imagejpeg($o, 'output_z_final.jpg', 100);
}

//================================== main ==================================



try_loading_a_cache_if_exists($i);


//count the pixels to heal
report('starting to count the pixels that need to be healed');

$pixels_to_heal = 0;

$skip_y = array();

for ($y = imagesy($i)-1; $y > -1; $y--){
   $th = 0;
   for ($x = imagesx($i)-1; $x > -1; $x--){       
          $pix = imagecolorat($i, $x, $y);
          if (maxcolororluminance($pix) < $GLOBALS['minlumi']){
             $pixels_to_heal++;
             $th++;
          }         
   }
   if ($th==0){
      $skip_y[] = $y;
      report('nothing to heal in line '.$y.', adding to ignore list');
   }
}

report('counted '.$pixels_to_heal.' pixels to heal, and save in '.$stages_to_save.' stages.<br />');
report(count($skip_y).' lines already being ignored');

$save_stage_at_x_heals = floor($pixels_to_heal / max($stages_to_save, 1))+1;

//iterate and heal

$healed = 0;


$maxx = imagesx($i);
$maxy = imagesy($i);

//$maxx = 300;
//$maxy = 300;
//$iterations = 10;

$st = 0;



$lumi_delta = 255 - $GLOBALS["maxlumi"];


$healed_in_this_round = 0;

$GLOBALS['maxlumi_orig'] = $GLOBALS["maxlumi"];

$skip_y_orig = $skip_y;
$st2 = 0;

$range_sx = array();
for ($sx = 0; $sx < imagesx($i); $sx = $sx + $GLOBALS['max_damage']){
   $range_sx[] = $sx;
};
shuffle($range_sx);
$range_sy = array();
for ($sy = 0; $sy < imagesy($i); $sy = $sy + $GLOBALS['max_damage']){
   $range_sy[] = $sy;
};
shuffle($range_sy);

//for random pairs of areas

$pairs = array();
foreach ($range_sy as $sy){
   foreach ($range_sx as $sx){
      $pairs[] = array($sx, $sy);
   }
};

shuffle($pairs);

for ($main_iteration = 0; $main_iteration < 2; $main_iteration++){
 foreach ($pairs as $pair){
    $sx = $pair[0];
    $sy = $pair[1];
    $GLOBALS['mincontext'] = 5; //start high then fall down
    $GLOBALS["maxlumi"] = $GLOBALS["maxlumi_orig"];
    $startx = $sx;
    $starty = $sy;
    $endy = min(imagesy($i), $starty + $GLOBALS['max_damage']);
    $endx = min(imagesx($i), $startx + $GLOBALS['max_damage']);
    $skip_y = $skip_y_orig;
    $to_heal_in_iteration = 1;//dummy
    $max_iterations = 50; // 10 still leaves holes
    $old_healed = $healed;
    for ($iterations = 0; $iterations < $max_iterations; $iterations++) if ($to_heal_in_iteration > 0){
        $h1 = $healed;   
        report("starting to heal [$startx, $starty, $endx, $endy]. Already healed = ".$healed.' of '.$pixels_to_heal.', '.
               $healed_in_this_round.' in the last round. Mincontext now = '.$GLOBALS['mincontext'].'<br />'); 
        
        //checking lines that are ready    
        $to_heal_in_iteration = 0;
        for ($y =$starty; $y<$endy; $y++) if (!in_array($y, $skip_y)){  
            $to_heal = 0;
            for ($x = $startx; $x<$endx; $x++){
                $pix = imagecolorat($i, $x, $y);
                if (maxcolororluminance($pix) < $GLOBALS['minlumi']){
                    $to_heal++;
                    $to_heal_in_iteration++;
                }            
            }
            if ($to_heal == 0){
               $skip_y[] = $y;
               report('added line '.$y.' to those that are ready and can be ignored<br />');
           }
        };
   
  
        $xrange = range($startx, $endx-1);
        shuffle($xrange);   
        $yrange = range($starty, $endy-1);
        shuffle($yrange);   
   
        $method = $iterations; //egyszer igy, egyszer ugy

        // don't rerun empty lines   
        report('starting to handle pixels');
        if ($method % 2 == 1){
             report('method 1 running<br />');
             foreach($yrange as $y) if (!in_array($y, $skip_y)){
                foreach($xrange as $x){
                   handle_this_pixel($i, $x, $y, $healed, $save_stage_at_x_heals, $st);
                }
             }
        }else{
             report('method 2 running<br />');
             foreach($xrange as $x){
                foreach($yrange as $y) if (!in_array($y, $skip_y)){
                   handle_this_pixel($i, $x, $y, $healed, $save_stage_at_x_heals, $st);
               }
             }
        }
        report('done');   
   
        //checking if the process' tolerance is low enough
        $healed_in_this_round = $healed - $h1;   
        if ($healed_in_this_round < 100){
           if ($GLOBALS['mincontext'] == 3){
               report('few pixels healed, mincontext already small so increasing lumi tolerance<br />');
               $GLOBALS['maxlumi'] = min(255, floor($GLOBALS['maxlumi']+$lumi_delta / 10));
           };          
           if ($GLOBALS['mincontext'] > 3){
              report('few pixels healed so decreased mincontext from '.$GLOBALS['mincontext'].'<br />');
              $GLOBALS['mincontext']--;
           };   
        }
   
        /* if ($iterations > 1){
           ////not doing this with the smaller areas 
           if ($GLOBALS['mincontext'] > 4){
              report('decreased mincontext from '.$GLOBALS['mincontext'].'<br />');
              $GLOBALS['mincontext']--;
           }
        }
        */

        if ($iterations > 3){
           if ($GLOBALS['mincontext'] > 3){
              report('decreased mincontext from '.$GLOBALS['mincontext'].'<br />');
              $GLOBALS['mincontext']--;
           }
        }
    };
    $st2++;
    $basename_oh = 'outpot_2_asregion_'.numberaskilo($st2);
    if ($old_healed != $healed){
        imagejpeg($i, $basename_oh.'.jpg');  
    }else{
        file_put_contents($basename_oh.'.txt', 'placeholder');
    }        
  };//for pairs
};//main_iteration
    
imagejpeg($i, 'output_8_before_lighterthan.jpg', 100);    
    
filter_lighterthanoriginal($i);    

// output a blinking interface to reveal the differences

echo '<img id = "img1" src="input.jpg" width="800" />';
echo '<img id = "img2" src="output_z_final.jpg?'.mt_rand().'" width = "800" style="display:none"/>';

?>
</div>
<script>
setInterval(function (){
   var a = document.getElementById("img1");
   var b = document.getElementById("img2");
   if (a.style.display!="none"){
      a.style.display = "none";
      b.style.display = "";
   }else{
      a.style.display = "";
      b.style.display = "none";
   }
}, 1000);
</script>

 

facebooktwittergoogle_plusredditpinterestlinkedinmail