This notebook shows how the complex imaging workflow associated with astrophotography can be automated using the Wolfram Language’s intrinsic image processing and data manipulation functionality. The galaxy images shown below were taking in May of 2023 from my somewhat light polluted backyard in Urbana, Illinois. The image of M101 also shows a supernova SN 2023ixf that was discovered just the prior week.
Implementation
Implementation
Compute the master bias frame from a series of FITS files. Bias frames are used to estimate the noise associated with reading the image out of the CCD chip.
Bias frames need to be taken with the sensor covered and at the shortest possible exposure time.
Bias frames need to be taken with the sensor covered and at the shortest possible exposure time.
ComputeMasterBias[theFiles_,outputName_]:=Module[{biasFrames,wmData,id,img},Print["Loading bias frames"];biasFrames=Import[#][1]&/@theFiles;Print["Finding the median of "<>ToString[Length[biasFrames]]<>" bias frames"];wmData=Median[Flatten/@(ImageData/@biasFrames)];id=ImageDimensions[biasFrames[[1]]];img=Image[Partition[wmData,id[[1]]],ColorSpace->"Grayscale"];Print["Exporting master bias frame: "<>outputName];Export[outputName,img];outputName];
Compute the master dark frame from a series of FITS files and the master bias frame. Dark frames are used to estimate the thermal noise associated with the imaging sensor at a given temperature. Dark frames need to be taken with the sensor covered at the same temperature and exposure time as the normal imaging frames.
ComputeMasterDark[theFiles_,biasFile_,outputName_]:=Module[{mbFrame,darkFrames,wmData,id,img,exposureTime},Print["Loading master bias frame"];mbFrame=Import[biasFile][1];Print["Loading dark frames"];darkFrames=Import[#][1]&/@theFiles;(*Hereweassumethatallthedarkframeshavethesameexposuretime*)exposureTime=Import[theFiles[[1]],{"FITS","ExposureTime"}][1];Print["Subtracting master bias noise from each dark frame"];darkFrames=ImageSubtract[#,mbFrame]&/@darkFrames;Print["Finding the median of "<>ToString[Length[darkFrames]]<>" dark frames"];wmData=Median[Flatten/@(ImageData/@darkFrames)];Print["Constructing master dark frame"];id=ImageDimensions[darkFrames[[1]]];img=Image[Partition[wmData,id[[1]]],ColorSpace->"Grayscale"];Print["Exporting master dark frame: "<>outputName];Export[outputName,img,{"FITS","ExposureTime"->exposureTime}];outputName];
Get a list of the best stars for alignment from a given image and max detection cutoff.
getAlignmentTargets[img_,max_]:=Module[{targets},targets=MaxDetect[img,max];targets=ComponentMeasurements[targets,{"Centroid","TotalIntensity"}];targets=ReverseSortBy[targets,(#[[2]])[[2]]&];#[[1]]&/@targets[[All,2]]];
Given a reference set of star positions and a corresponding target set of star positions, find the rigid transform between the two sets of coordinates within a given error tolerance.
getAlignmentTransform0[ref_List,target_List,maxError_]:=Module[{reflen,targlen,refs,targ,i,xform,failed},reflen=Length[ref];If[reflen<2,Print["found too few ("<>ToString[reflen]<>") reference points"];Return[$Failed]];reflen=Min[3,reflen];refs=Take[ref,reflen];targlen=Length[target];If[targlen<reflen,Print["found too few ("<>ToString[targlen]<>") target points"];Return[$Failed]];targ={};For[i=1,i<=reflen,i++,AppendTo[targ,MinimalBy[target,EuclideanDistance[refs[[i]],#]&][[1]]]];xform=Quiet[FindGeometricTransform[refs,targ,TransformationClass->"Rigid"]];(*Print["xform error = "<>ToString[AccountingForm@xform[[1]]]];Print["max error allowed = "<>ToString[AccountingForm@maxError]];*)failed=(xform[[1]]>maxError)||!RealValuedNumberQ[xform[[1]]];If[failed,Return[$Failed]];xform[[2]]];getAlignmentTransform[ref_,img_,cutoff_,maxError_]:=Module[{refTargets,refTriangle,imgTargets,imgTriangle,a,b,c,reflen,imglen},refTargets=getAlignmentTargets[ref,cutoff];imgTargets=getAlignmentTargets[img,cutoff];getAlignmentTransform0[refTargets,imgTargets,maxError]];
Align an input astro image with a reference astro image based on the best transform between corresponding stars.
alignAstroImageWL[ref_,inp_,cutoff_,maxError_]:=Module[{reflocs,inplocs,xform},xform=getAlignmentTransform[ref,inp,cutoff,maxError];If[xform==$Failed,Print["Failed to find transform"];Return[$Failed]];ImagePerspectiveTransformation[inp,xform,DataRange->Full]];
Given a set of image FITS files, a master bias FITS file, a master dark FITS file, create a set of aligned FITS files in the specified output directory.
Each image file is calibrated by subtracting the master bias and dark noise and then each image is normalized before alignment to convert from raw electron counts to a [0,1] value range.
Each image file is calibrated by subtracting the master bias and dark noise and then each image is normalized before alignment to convert from raw electron counts to a [0,1] value range.
AlignFileSequence[theFiles_,biasFile_,darkFile_,outputDir_,maxAlignmentError_,drizzle_]:=Module[{darkframe,biasframe,targetframe,thisframe,name,fullname,f},(*Loadthemasterdarkframe*)darkframe=Import[darkFile][1];Print["dark frame: "<>darkFile];(*Loadthemasterbiasframe*)biasframe=Import[biasFile][1];Print["bias frame: "<>biasFile];(*Takethetargetforalignmenttobethefirstofthelightframes*)targetframe=Import[theFiles[[1]]][1];Print["target frame: "<>theFiles[[1]]];(*Calibratethetargetframe*)targetframe=ImageSubtract[targetframe,biasframe];targetframe=ImageSubtract[targetframe,darkframe];targetframe=ImageAdjust[targetframe];If[drizzle>1.0,targetframe=ImageResize[targetframe,Scaled[drizzle],Resampling->Automatic]];(*Writeoutthetargetframe*)name=Last[FileNameSplit[theFiles[[1]]]];fullname=FileNameJoin[{outputDir,name}];Print["writing out reference image to "<>fullname];Export[fullname,targetframe];For[f=2,f<=Length[theFiles],f++,thisframe=Import[theFiles[[f]]][1];thisframe=ImageSubtract[thisframe,biasframe];thisframe=ImageSubtract[thisframe,darkframe];thisframe=ImageAdjust[thisframe];If[drizzle>1.0,thisframe=ImageResize[thisframe,Scaled[drizzle],Resampling->Automatic]];thisframe=alignAstroImageWL[targetframe,thisframe,.05,maxAlignmentError];If[thisframe==$Failed,Print["Skipping image"];Continue[]];name=Last[FileNameSplit[theFiles[[f]]]];fullname=FileNameJoin[{outputDir,name}];Print["writing out transformed image to "<>fullname];Export[fullname,thisframe];];outputDir];
Estimate and remove sky glow from an image, using npts sample points
RemoveSkyglowGS[im_,npts_]:=Module[{pts,skyGlow},pts=Table[{RandomInteger[{1,ImageDimensions[im][[1]]}],RandomInteger[{1,ImageDimensions[im][[2]]}]},{npts}];skyGlow=Mean[Map[ImageValue[im,#]&,pts]];ImageSubtract[im,skyGlow]];RemoveSkyglowRGB[im_,npts_]:=ColorCombine[RemoveSkyglowGS[#,npts]&/@ColorSeparate[im,"RGB"],"RGB"];RemoveSkyglow[im_,npts_:500]:=If[ImageColorSpace[im]=="Grayscale",RemoveSkyglowGS[im,npts],RemoveSkyglowRGB[im,npts]];
Implementation of the IRAS ZSCALE algorithm for stretching image values.
ZScaleImageC[im_,npoints_Integer,contrast_Real,values_List]:=Module[{fit,midpoint,inten,z1,z2},midpoint=values[[npoints/2]];fit=Association[FindFit[values,intercept+slope*(v-midpoint),{intercept,slope},v]];inten[v_]:=fit[intercept]+fit[slope]*(v-midpoint);z1=inten[midpoint]+(fit[slope]/contrast)(1.0-midpoint);z2=inten[midpoint]+(fit[slope]/contrast)(npoints-midpoint);ImageAdjust[im,{0.0,0.0},{z1,z2},{0.0,1.0}]];ZScaleStretchImage[im_,contrast_Real:0.5,npoints_Integer:1000]:=Module[{width,height,points,values,idim=ImageDimensions[im]},width=idim[[1]];height=idim[[2]];points=Table[{RandomInteger[{0,width}],RandomInteger[{0,height}]},npoints];values=Sort[ImageValue[im,points]];If[contrast==0.0,ImageAdjust[im,{0.0,0.0},{values[[1]],values[[npoints]]},{0.0,1.0}],ZScaleImageC[im,npoints,contrast,values]]];
Given a set of (light) image FITS files, a master bias FITS file, a master dark FITS file, this routine will calibrate, align and stack the image files to increase the signal-to-noise ratio (SNR) of the observed object. This will discard a the fract fraction of the outlying pixels and replace them with the next most popular values using the default Winsorized mean combineMethod. Other combine methods can be specified. Before combining the images, their histograms are brought into alignment so a meaningful average can be determined.
alignAndStack0[imageFiles_List,mbFile_String,mdFile_String,fract_,drizzle_,combineMethod_,alignmentError_]:=Module[{tempDir,alignedImages,refImage,wmData,id},Print["Aligning "<>ToString[Length[imageFiles]]<>" raw files"];tempDir=AlignFileSequence[imageFiles,mbFile,mdFile,CreateDirectory[],alignmentError,drizzle];Print["Importing aligned files"];alignedImages=Import/@FileNames[FileNameJoin[{tempDir,"*.fit"}]];(*Cleanupalignedfiles*)DeleteDirectory[tempDir,DeleteContents->True];If[Length[alignedImages]==0,Print["Failed to align any images."];Return[Null]];alignedImages=#[1]&/@alignedImages;(*NormalizeimportedFITSdatatocreateviewableimages*)(*alignedImages=ImageAdjust/@alignedImages;*)(*AdjustallthealignedimagessotheyhavesimilarhistogramstothefirstimageThisisnecessarybeforewecancombinetheimages*)Print["Applying histogram transforms"];refImage=alignedImages[[1]];alignedImages=HistogramTransform[#,refImage]&/@alignedImages;Print["Computing mean; replacing the outlying "<>ToString[100*fract]<>"% pixels."];(*Computeawinsorizedmeantoeliminateoutlyingbadpixels*)wmData=combineMethod[Flatten/@(ImageData/@alignedImages),fract];Print["Constructing final image"];id=ImageDimensions[alignedImages[[1]]];Image[Partition[wmData,id[[1]]],ColorSpace->"Grayscale"]];Options[AlignAndStack]={TrimFraction->0.05,CombineMethod->WinsorizedMean,Drizzle->1.0,MaxTransformError->0.75};Options[ImagePipeline]={MaxNoise->0.01,Contrast->0.5,TrimFraction->0.05,CombineMethod->WinsorizedMean,Drizzle->1.0,MaxTransformError->0.75};
Main alignment and stacking entry point
Main imaging pipeline. This aligns and stacks the monochrome image (light) FITS files after calibrating them using the given master bias and master dark FITS files.
After alignment and stacking, noise is reduced using the TotalVariationFilter and then the image is automatically stretched using the IRAS ZScale algorithm.
Finally, sky glow is estimated and subtracted, before a final image adjust.
After alignment and stacking, noise is reduced using the TotalVariationFilter and then the image is automatically stretched using the IRAS ZScale algorithm.
Finally, sky glow is estimated and subtracted, before a final image adjust.
Imaging Pipeline Examples
Imaging Pipeline Examples
Bias Frame Computation
Bias Frame Computation
Bias frames record the noise associated with reading the imaging data out of the CCD array.
Dark Frame Computation
Dark Frame Computation
Dark frames record the noise associated with the temperature (0° C) of the CCD array.
Galaxy M51 (The Whirlpool Galaxy)
Galaxy M51 (The Whirlpool Galaxy)
Final stretching is usually done in GIMP using the Levels and Curves tools.
Galaxy M101 (Plus Supernova SN 2023ixf near the top of frame)
Galaxy M101 (Plus Supernova SN 2023ixf near the top of frame)
The arrow indicates the position of SN 2023ixf
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Deep space object imaging pipeline examples: M51 & M101 galaxies
by Tom Sherlock
Wolfram Community, STAFF PICKS, August 3, 2023
https://community.wolfram.com/groups/-/m/t/2982964
by Tom Sherlock
Wolfram Community, STAFF PICKS, August 3, 2023
https://community.wolfram.com/groups/-/m/t/2982964