// Copyright ©2025-2026 Karel Kulhavý
//
// blacknoise.ajla is free software: you can redistribute it and/or modify it under
// the terms of the GNU General Public License as published by the Free
// Software Foundation, either version 3 of the License, or (at your option)
// any later version.
//
// blacknoise.ajla is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with
// blacknoise.ajla. If not, see .
// How to Use GNU Licenses for Your Own Software:
// https://www.gnu.org/licenses/gpl-howto.html
const lf:=bytes.[10];
type ctype:=real; // ctype like Calculation TYPE for large pixel calculation in arrays.
const blur_stdev_for_white_level:=0.5; // Standard deviation of the Gaussian (single-sided).
const gaussian_doublecycles_spread:=2; // 3=classical Gaussian approximation, 1 is not separable (horizontal and vertical smears but not much diagonal), 2 still
// suffers from a bit separability problem
const gaussian_doublecycles_white_level:=3; // 3=classical Gaussian approximation, 1 is not separable (horizontal and vertical smears but not much diagonal), 2 still
// suffers from a bit separability problem
uses spawn;
uses charset;
record image[
w:int;
h:int;
d:list(ctype);
]
fn emit_raw_sRGB_16bpc(implicit w: world, h:handle, i:image):world
[
var xc y c:int;
var i_d:=i.d;
var row_subpixels:=i.w*3;
var blob:bytes:=fill(cast(byte,0),row_subpixels*2);
for y:=0 to i.h do [
for xc:=0 to row_subpixels do [
var val:ctype;
val:=i_d[y*row_subpixels+xc];
val:=clip(val,0,1);
{ sRGB encoding. https://en.wikipedia.org/wiki/SRGB inverse OETF }
if val <= 0.0031308 then val := 12.92 * val;
else val := 1.055*power(val,1.0/2.4)-0.055;
val *= 65535;
val += 0.5; { For rounding to nearest }
{ MSB first: https://en.wikipedia.org/wiki/Netpbm#File_formats}
blob[xc*2+0]:=cast(int,val) shr 8;
blob[xc*2+1]:=cast(int,val) and 255;
]
write(h,blob);
]
]
fn emit_ppm(implicit w: world, d:dhandle, fn_ppm:bytes, i:image):world
[
var h:=wopen(d,fn_ppm, open_flag_create, open_mode_default);
write(h,"P6"+lf+
ntos(i.w)+" "+ntos(i.h)+lf+
"65535"+lf);
emit_raw_sRGB_16bpc(w,h,i);
]
fn runcmd(implicit w: world, cwd:dhandle, handles: list(handle), shell_command:bytes, env : treemap(bytes, bytes)):world
[
var shell, cmdline := spawn_command(shell_command,env);
var ph:= spawn(cwd, handles, shell, cmdline, env);
if is_exception ph then abort forced_exit(1,"Error spawning the following program. Could be because it is not installed, not visible in the PATH"
+" environment variable, it has problem with libraries, disk is full, RAM is full, cannot find a data file, "
+"a data file has wrong format or other reasons. Command executed: "
+shell_command+" . Error: "+exception_string ph);
var exit_code := wait(ph);
if exit_code < 0 then abort forced_exit(1,"Command "+shell_command+" started, but was killed by signal " + ntos(-exit_code)+". Could be because it ran "
+"out of RAM and was killed by the system or crashed due to its bug.");
]
fn stdout_from_cmd(implicit w: world, d:dhandle, h: list(handle), shell_command:bytes, env : treemap(bytes, bytes)):(world,bytes)
[
var shell, cmdline := spawn_command(shell_command,env);
var data := get_process_output(d, h, shell, cmdline, env);
if is_exception data then abort forced_exit(1,"Error spawning the following program. Could be because it is not installed, not visible in the PATH"
+" environment variable, it has problem with libraries, disk is full, RAM is full, cannot find a data file, "
+"a data file has wrong format or other reasons. Command executed: "
+shell_command+" . Error: "+exception_string data +nl);
return w,data;
]
fn eat_whitespaces(s: bytes):bytes
[
var c: byte;
again:
if not len_at_least(s,1) then return s;
c:=s[0];
if c=9 or c=10 or c=13 or c=32 then [
s:=s[1..];
goto again;
]
return s;
]
fn eat_number(s: bytes):(bytes,int)
[
var c: byte;
var n:=0;
again:
if not len_at_least(s,1) then return s,n;
c:=s[0];
if c>='0' and c<='9' then [
n:=n*10+c-'0';
s:=s[1..];
goto again;
]
return s,n;
]
fn eat_ws_number(s: bytes):(bytes,int)
[
s:=eat_whitespaces(s);
var n:int;
s,n:=eat_number(s);
return s,n;
]
// Returns the rest of data as 2nd return value.
fn PPM_image_data2image(img_w img_h maxval: int,data:bytes):(image, bytes)
[
var samples:=img_w*img_h*3;
var img_d:=fill(cast(ctype,0.0),samples);
var x y c:int;
var bytes_consumed:int;
{ data := list_flatten(data); This is counterproductive, increases both RAM consumption and run time. }
var one_over_maxval:=1.0/maxval;
if maxval>=256 then [
bytes_consumed:=samples*2;
var xdata := data[ .. bytes_consumed];
xdata := list_flatten(xdata);
for sample:=0 to samples do [
img_d[sample]:=(cast(ctype,xdata[sample*2])*256.0+xdata[sample*2+1])*one_over_maxval;
]
]else[
bytes_consumed:=samples;
var xdata := data[ .. bytes_consumed];
xdata := list_flatten(xdata);
for sample:=0 to samples do [
img_d[sample]:=cast(ctype, xdata[sample])*one_over_maxval;
]
]
return image.[w:img_w, h:img_h,d:img_d], data[bytes_consumed..];
]
// If images have at least one dimension different, whole Ajla crashes with an assertion failed.
fn add_images (i j: image):image
[
xeval assert(i.h=j.h,"Different widths of images when trying to add them");
xeval assert(i.w=j.w,"Different heights of images when trying to add them");
var n_subpixels:=i.h*i.w*3;
for subpixel:=0 to n_subpixels do i.d[subpixel]+=j.d[subpixel];
return i;
]
// returns i-j
// If images have at least one dimension different, whole Ajla crashes with an assertion failed.
fn sub_images (i j: image):image
[
xeval assert(i.h=j.h,"Different widths of images when trying to add them");
xeval assert(i.w=j.w,"Different heights of images when trying to add them");
var n_subpixels:=i.h*i.w*3;
for subpixel:=0 to n_subpixels do i.d[subpixel]-=j.d[subpixel];
return i;
]
// If images have at least one dimension different, whole Ajla crashes with an assertion failed.
fn mul_images (i: image, mul: ctype):image
[
var n_subpixels:=i.h*i.w*3;
for subpixel:=0 to n_subpixels do i.d[subpixel]*=mul;
return i;
]
fn stdev2decay_multiplier(stdev: ctype, stages: int): ctype
[
var decay_time:= stdev/sqrt(cast(ctype,stages));
return 1.0-exp(-1.0/decay_time);
]
// s is strength, 0=no change, 0.2=maybe reasonable
// 0=infinite blur, small number >0=large blur, large number <1=small blur, 1=no blur, >1...<2=ringing, 2=ringing never stops, >2 divergent ringing
fn blur (i: image, stdev: ctype,gaussian_double_cycles:int): image
[
var x xr y yr c cycle:int;
var capacitor:ctype;
var i_d:=i.d;
var i_h:=i.h;
var i_w:=i.w;
if stdev <=0 then return i;
var mult:=stdev2decay_multiplier(stdev,2*gaussian_double_cycles);
for c:=0 to 3 do [
for y:=0 to i_h do [
x:=0;
capacitor:=i_d[(i_w*y+x)*3+c];
for cycle:=0 to gaussian_double_cycles do [
for x:=0 to i_w do [
var dst:=i_d[(i_w*y+x)*3+c];
capacitor+=(dst-capacitor)*mult;
i_d[(i_w*y+x)*3+c]:=capacitor;
]
for xr:=0 to i_w do [
var x:=i_w-xr-1;
var dst:=i_d[(i_w*y+x)*3+c];
capacitor+=(dst-capacitor)*mult;
i_d[(i_w*y+x)*3+c]:=capacitor;
]
]
]
for x:=0 to i_w do [
y:=0;
capacitor:=i_d[(i_w*y+x)*3+c];
for cycle:=0 to gaussian_double_cycles do [
for y:=0 to i_h do [
var dst:=i_d[(i_w*y+x)*3+c];
capacitor+=(dst-capacitor)*mult;
i_d[(i_w*y+x)*3+c]:=capacitor;
]
for yr:=0 to i_h do [
var y:=i_h-yr-1;
var dst:=i_d[(i_w*y+x)*3+c];
capacitor+=(dst-capacitor)*mult;
i_d[(i_w*y+x)*3+c]:=capacitor;
]
]
]
]
i.d:=i_d;
return i;
]
fn print_img_stats (implicit w: world, h: list(handle), i: image, ident:bytes): world
[
var r_min:=i.d[0];
var r_max:=i.d[0];
var g_min:=i.d[0];
var g_max:=i.d[0];
var b_min:=i.d[0];
var b_max:=i.d[0];
var pixels:=i.h*i.w;
for pixel:=1 to pixels do [
var r:=i.d[pixel*3+0];
var g:=i.d[pixel*3+1];
var b:=i.d[pixel*3+2];
if rr_max then r_max:=r;
if g>g_max then g_max:=g;
if b>b_max then b_max:=b;
]
var rgb_min:=min(r_min, min(g_min, b_min));
var rgb_max:=max(r_max, max(g_max, b_max));
write(h[2],ident+
": "+ntos(i.w)+"x"+ntos(i.h)+" RGB "+ntos_base_precision(rgb_min,10,1)+" to "+ntos_base_precision(rgb_max,10,1)+", "+
"R "+ntos_base_precision( r_min,10,1)+" to "+ntos_base_precision( r_max,10,1)+", "+
"G "+ntos_base_precision( g_min,10,1)+" to "+ntos_base_precision( g_max,10,1)+", "+
"B "+ntos_base_precision( b_min,10,1)+" to "+ntos_base_precision( b_max,10,1)+"." +nl);
]
// Returns peak if nonzero, otherwise 1 as peak.
// s is strength, 0=no change, 0.2=maybe reasonable
fn spread_negatives (implicit w: world, h:list(handle), i: image, blur_end:ctype, steps: int): (world, image)
[
var subpixels:=i.h*i.w*3;
for step := 0 to steps do [
var ratio:ctype;
if steps <= 1 then ratio:=1; else ratio:=cast(real, step)/(steps-1);
var current_blur:=power(blur_end,ratio);
var neg_img:=i;
neg_img.d:=fill(cast(ctype,0),subpixels);
for subpixel:=0 to subpixels do [
if i.d[subpixel]<0 then [
neg_img.d[subpixel]:=i.d[subpixel];
i.d[subpixel]:=0;
]
]
neg_img:=blur(neg_img, current_blur, gaussian_doublecycles_spread);
i := add_images(i, neg_img);
print_img_stats(h,i,"After spreading with stdev "+ntos_base_precision(current_blur,10,1)+" px");
]
print_img_stats(h,i,"After spreading of negatives");
return w, i;
]
// Normalizes to 0-1 range
fn normalize (i: image, peak: ctype):image
[
var inv_peak:=cast(ctype,1)/peak;
var i_d:=i.d;
var length:=i.h*i.w*3;
for i:=0 to length do i_d[i]*=inv_peak;
i.d:=i_d;
return i;
]
fn nonzero_peak_min_from_RGB(i: image):ctype
[
var i_d:=i.d;
var pixels:=i.h*i.w;
var peak_r:ctype:=0;
var peak_g:ctype:=0;
var peak_b:ctype:=0;
for pixel:=0 to pixels do [
var r:=i_d[pixel*3+0];
var g:=i_d[pixel*3+1];
var b:=i_d[pixel*3+2];
if r>peak_r then peak_r:=r;
if g>peak_g then peak_g:=g;
if b>peak_b then peak_b:=b;
]
var peak:=min(peak_r,min(peak_g,peak_b));
if peak <=0 then peak:=1; { Prevent division by 0. This can still happen because float numbers are unreliable and can change value a bit }
return peak;
]
fn dng2white_balance(implicit w:world, d:dhandle, h:list(handle), env:treemap(bytes,bytes),dng_filename: bytes): (world, bytes)
[
var prefix_:="Camera multipliers: ";
var cmdline:="dcraw -i -v -w "+dng_filename+"|LC_ALL=C grep '^"+prefix_+"'";
var data:bytes;
write(h[2],"Running this dcraw commandline to extract white balance multipliers from the photo DNG:"+nl+cmdline+nl);
w,data:=stdout_from_cmd(d,h,cmdline,env);
data :=data[len(prefix_)..];
data:=list_filter(data,lambda (b:byte):bool[return b>=32;]);
return data;
]
// Returns the rest of data as 2nd return value.
fn PPM2image(implicit w: world, h: list(handle), data:bytes, ident: bytes ): (world, image, bytes)
[
var img_w img_h maxval:int;
data,img_w:=eat_ws_number(data[2..]);
data,img_h:=eat_ws_number(data);
data,maxval:=eat_ws_number(data);
write(h[2],"Loading "+ident+" PPM image "+ntos(img_w)+"x"+ntos(img_h)+" ("+ntos_base_precision(cast(real,img_w)*img_h/1e6,10,1)+
" Mpix), "+ntos(maxval+1)+" levels from dcraw."+nl);
data:=data[1..];
var img, rest:=PPM_image_data2image(img_w, img_h, maxval, data);
return w,img,rest;
]
// Like PPM2image, but eats up all images, averages them and doesn't return any rest of data.
// Each input PPM is normalized to white level 1 according to its own maxval, then they are averaged, at the end they are normalized to white level 65535.
fn PPMs2image(implicit w: world, h: list(handle), data:bytes, ident:bytes): (world, image)
[
var i j:image;
w, i, data :=PPM2image(w, h, data, ident);
var n_images:=1;
more:
if not len_greater_than(data, 0) then
goto exit_;
//if len (data) <= 0 then goto exit_;
w, j, data :=PPM2image(w, h, data, ident);
i := add_images(i, j);
n_images += 1;
goto more;
exit_:
i:=mul_images(i,cast(ctype,65535)/n_images);
return w,i;
]
fn reduce_resolution(i:image, divisor: int):image
[
var o:=i;
var counts:list(ctype);
if divisor <=1 then return i;
o.h:=(i.h-1) div divisor+1;
o.w:=(i.w-1) div divisor+1;
var out_pixels:=o.h*o.w;
var out_subpixels:=out_pixels*3;
counts:=fill(cast(ctype,0),out_pixels);
o.d:=fill(cast(ctype,0),out_subpixels);
for y:=0 to i.h do [
var oy:=y div divisor;
for x:=0 to i.w do [
var ox:=x div divisor;
counts[o.w*oy+ox]+=1;
for c:=0 to 3 do o.d[(o.w*oy+ox)*3+c]+=i.d[(i.w*y+x)*3+c];
]
]
for oy:=0 to o.h do
for ox:=0 to o.w do
for c:=0 to 3 do
o.d[(o.w*oy+ox)*3+c]/=counts[o.w*oy+ox];
return o;
]
// flag will be either "-q 3" (full resolution, 66.666...% or more) or "-h" (half resolution, <66.666...%).
fn percentage2h_divisor(percentage: real):(flag: bytes, divisor : int, effective_percentage:real)
[
if percentage >= 100.0/1.5 then return "-q 3",1,100;
if percentage <= 0 then abort forced_exit(1,"Error: Cannot have zero or less resolution percentage.");
var divisor:=round(50.0/percentage);
return "-h",divisor,50.0/divisor;
]
fn main
[
if len (args)<6 then [
write(h[2],""
+nl+"Usage: ajla blacknoise.ajla rot_flip resolution_percent blur_end steps dark_frames_dir filename.dng"
+nl+""
+nl+"Example: portrait document: ajla blacknoise.ajla 90 17 64 25 portrait.dng"
+nl+" landscape document: ajla blacknoise.ajla 0 17 64 25 landscape.dng"
+nl+" lower resolution : ajla blacknoise.ajla 90 17 64 25 portrait.dng"
+nl+""
+nl+"rot_flip is """",0-7,90,180,270. It is unreliable, can be turned upside down in"
+nl+" various phones because it is passed to dcraw which lacks specification:"
+nl+" """" landscape"
+nl+" 0 landscape"
+nl+" 1 mirrored landscape (flipped left and right)"
+nl+" 2 landscape that was first mirrored (flipped left and right) and then"
+nl+" turned upside down"
+nl+" 3 upside down landscape"
+nl+" 4 portrait that was first mirrored (flipped left and right) and then"
+nl+" turned upside down"
+nl+" 5 upside down portrait"
+nl+" 6 portrait"
+nl+" 7 mirrored portrait (flipped left and right) - often used for selfies"
+nl+" 90 portrait"
+nl+" 180 upside down landscape"
+nl+" 270 upside down portrait"
+nl+""
+nl+"resolution_percent can be >0 to 100. For 66.7% up, full resolution will be"
+nl+" used. For less than 66.7%, the resolution will first be"
+nl+" halved in the dcraw program to 50% using the -h flag, which"
+nl+" turns off the demosaicing algorithm, which I believe is ben-"
+nl+" eficial because it won't be confused by the strong noise."
+nl+" For 33.3% and less, there will be additional division of the"
+nl+" Resolution by an integer (2,3,4...) so the resulting exact"
+nl+" percentages achievable will be 100%, 50%, 25%, 16.7%, 12.5%,"
+nl+" 10%, 8.3%, 7.1%, 6.3%, 5.6%, 5%,..."
+nl+""
+nl+"blur_end Standard deviation of the end Gaussian blur. Recommended 64. Start blur"
+nl+" is 1. Proceeds logarithmically."
+nl+""
+nl+"steps number of blur steps to distribute negative values. Recommended 25. First"
+nl+" step has standard deviation of Gaussian blur 1, last step blur_end."
+nl+" Proceeds logarithmically."
+nl+""
+nl+"Output is a 16-bit per pixel PPM that uses sRGB color space and transfer function,"
+nl+" .dng in the filename is replaced with .ppm."
+nl+""
+nl);
return w;
]
var rot_flag:bytes;
[
rot_flag:=args[0];
if len(rot_flag)<=0 then rot_flag:="0"; {Prevent empty flag causing inconsistent rotation of picture and dark frame when they were shot
with different orientation}
if len(rot_flag)>0 then rot_flag:="-t "+rot_flag;
]
var blur_end:ctype:=instance_real_number_real64.from_bytes(args[2]);
var steps:=ston(args[3]);
var dark_dir:=args[4];
var fn_dng:=args[5];
var basefn:bytes;
[
var suffix:=".dng";
if list_ends_with(bytes_locase(fn_dng),bytes_locase(suffix)) then [
basefn:=fn_dng[..len(fn_dng)-len(suffix)];
]else[
basefn:=fn_dng;
]
]
var fn_ppm:=basefn+".ppm";
var white_balance:bytes;
w,white_balance:=dng2white_balance(w,d,h,env,fn_dng);
write(h[2],"White balance multipliers which will be used both for image and dark frames: "+white_balance+nl);
var resampling_flag, resolution_divisor, effective_percentage:=percentage2h_divisor(ston(args[1]));
write(h[2],"Requested resolution percentage: "+args[1]+", rounded to resolution percentage "+ntos_base_precision(effective_percentage,10,2)
+"%, which is implemented by running dcraw with flag "+resampling_flag+", which first "
+select(resampling_flag="-h","keeps full resolution and turns AHD debayering algorithm on","halves resolution and turns debayering off")
+select(resolution_divisor<=1," and then the resolution is further reduced by a factor of "+ntos(resolution_divisor),"and the resolution is not reduced anymore")
+"."+nl);
var data:bytes;
{ -W is necessary to prevent different scaling of image and dark data.
-4 is necessary because we need photolinear.
-H 1 is necessary because -H 0 wouldn't work anyway, subtracting the dark frame would cause the white highlights to turn green.
-w must not be present because camera white balance can be different for image and dark frame data.
-j I don't include because I don't want the document to be stretched on cameras with non-square pixels.}
var dcraw_opts:="-c -k 0 "+rot_flag+" -r "+white_balance+" -H 1 -c -v -4 -W "+resampling_flag;
var cmdline:="dcraw "+dcraw_opts+" "+fn_dng;
write(h[2],"Trying to run this commandline to convert DNG into 16-bit PPM and load it via a pipe from:"+nl+cmdline+nl);
w,data:=stdout_from_cmd(d,h,cmdline,env);
var img:=PPMs2image(h, data, "photograph");
print_img_stats(h,img,"Image loaded from DNG");
cmdline:="parallel ""dcraw "+dcraw_opts+""" ::: """+dark_dir+"/""*.[Dd][Nn][Gg]";
write(h[2],"Trying to run this commandline to convert dark frame DNG's into 16-bit PPM's and load them via a pipe from:"+nl+cmdline+nl);
w,data:=stdout_from_cmd(d,h,cmdline,env);
var black_frame:=PPMs2image(h, data, "dark frame");
print_img_stats(h,black_frame,"Dark frame loaded from DNGs");
img := sub_images(img, black_frame);
print_img_stats(h,img,"After dark frame subtraction");
if resolution_divisor <=1 then write(h[2],"Will not reduce resolution further."+nl);
else[
write(h[2],"Reducing resolution "+ntos(resolution_divisor)+" times."+nl);
img:=reduce_resolution(img, resolution_divisor);
print_img_stats(h,img,"After reducing resolution "+ntos(resolution_divisor)+" times");
]
write(h[2],"Starting spreading negatives with Gaussian standard deviation from 1 to "+ntos_base_precision(blur_end,10,1)+
" pixels in "+ntos(steps)+" steps."+nl);
img:=spread_negatives(w, h, img, blur_end, steps);
var blurred_image:=blur(img, blur_stdev_for_white_level,gaussian_doublecycles_white_level);
var white_level:=nonzero_peak_min_from_RGB(blurred_image);
write(h[2],"Blurring finished. Starting normalizing the white level to "+ntos_base_precision(white_level,10,1)+"."+nl);
img:=normalize(img, white_level);
write(h[2],"Normalizing fininshed. Writing a "+ntos(img.w)+"x"+ntos(img.h)+" 16-bit PPM into "+fn_ppm+"."+nl);
emit_ppm(d,fn_ppm,img);
return w;
]