; file: fitscube2mpeg.pro = make mpeg movie from an image cube in a fits file ; init: Oct 27 2002 from /home/rutten/rr/wrk/SST/2002-09-xx makempeg.pro ; last: Nov 23 2015 Rob Rutten Oslo pro fitscube2mpeg,cubefile,mpegfile,$ xrange=xrange,yrange=yrange,trange=trange,$ rebin=rebin,cutxga=cutxga,sharpen=sharpen,$ zeromid=zeromid,sqrtint=sqrtint,logint=logint,$ bytscale=bytscale,cutmin=cutmin,cutmax=cutmax,$ wrappercent=wrappercent,wraprms=wraprms,$ clock=clock,starttime=starttime,cadence=cadence,$ datebanner=datebanner,diagname=diagname,jpeg=jpeg,$ gap=gap,overlay=overlay,verbose=verbose ;+ ; NAME: ; fitscube2mpeg.pro ; PURPOSE: ; make an mpeg movie from an image sequence in a fitsfile [x,y,t] cube ; CALL: ; fitscube2mpeg,cubefile,mpegfile,$ ; xrange=xrange,yrange=yrange,trange=trange,$ ; rebin=rebin,cutxga=cutxga,sharpen=sharpen,$ ; zeromid=zeromid,sqrtint=sqrtint,logint=logint,$ ; bytscale=bytscale, cutmin=cutmin,cutmax=cutmax,$ ; wrappercent=wrappercent,wraprms=wraprms,$ ; clock=clock,starttime=starttime,cadence=cadence,$ ; datebanner=datebanner,diagname=diagname,jpeg=jpeg,$ ; gap=gap,overlay=overlay,verbose=verbose ; EXAMPLE CALL: ; at the bottom of the program file ; INPUTS: ; cubefile = string wth path/filename including .fits of input fitscube ; KEYWORDS: ; xrange, yrange, trange = partial cutout specification (2-element arrays) ; cutxga = 1/0: cut to 1024x768 XGA size for beamers ; rebin = multiplier to nx and ny for rebinning (larger or smaller image) ; cutxga = cut frame size to 1024x768 (beamer size XGA) ; sharpen = 1/0: sharpen with fixed ImageMagick parameters (for talks) ; zeromid = +1/0: keep zero at byte 127 (magnetogram, Dopplergram) ; sqrtint = 1/0 display sqrt(int) ; logint = 1/0: display alog10(int) ; bytscale: options for bytscaling the movie ; bytscale = 0 don't bytscale the data (should be bytes) ; bytscale = 1 use (min,max) of first image for whole movie ; bytscale = 2 use (min,max) of the whole image sequence (default) ; bytscale = 3 bytscale every movie frame individually (not for /zeromid) ; bytscale = 4 use cutmin and cutmax values for this movie ; cutmin = minimum value for bytscale=4 ; negative = - percentage fraction of min of the whole sequence ; for mag and dop cutmin is set to -cutmax to maintain zero ; cutmax = maximum value for bytscale=4 ; negative = - percentage fraction of max of the whole sequence ; wrappercent = wraparound at mean + (value/100)*mean ; wraprms = wraparound at mean + value*rms ; clock = 1/0: insert running clock ; starttime: specify if no header-keyword starttim ; standard anytim format, eg '2014-01-30 12:11' ; cadence: specify (in s) if no header-keyword cadence ; datebanner = 1/0: insert running date-time specifier ; diagname: insert when specified, e.g. 'AIA 171' ; jpeg = 1/0: write first image as jpeg (without clock etc) ; gap = nr frames of grey at end (to mark start of replay) ; overlay = bytarr (nx,ny) frame or bytarr (nx,ny,nt) cube with overlay ; verbose = 1/0: print progress ; OUTPUTS: ; mpegfile = string with path+filename including .mpg for output file ; RESTRICTIONS: ; requires a private tmp directory: ~/tmp ; cuts spatial extent to multiples of 16 pixels ; METHOD: ; uses assoc to accommodate large fitscube files ; suited for call by doallfiles.pro to process multiple cubes ; HISTORY: ; Nov 18 2013 RR: from cube2mpeg ; Jan 21 2014 RR: option to sharpen ; Jan 31 2014 RR: corrections, add datebanner, diagname, cutxga ; Feb 4 2014 RR: add zeromid, wrappercent, wraprms ; Feb 8 2014 RR: RGB > YUV conversion (apparently necessary) ; Aug 22 2014 RR: options cutmin, cutmax ; Sep 5 2014 RR: options gap, overlay ;- ; answer no-parameter query if (n_params() lt 2) then begin print,'fitscube2mpeg,cubefile,mpegfile,$' print,' [xrange=xrange,yrange=yrange,trange=trange,$' print,' rebin=rebin,cutxga=cutxga,$' print,' zeromid=zeromid,sqrtint=sqrtint,logint=logint,bytscale=bytscale,$' print,' wrappercent=wrappercent,wraprms=wraprms,$' print,' clock=clock,starttime=starttime,cadence=cadence,$' print,' datebanner=datebanner,diagname=diagname,jpeg=jpeg,$' print,' gap=gap,overlay=overlay,verbose=verbose]' return endif ; keyword defaults if (n_elements(xrange) eq 0) then xrange=[0,-1] if (n_elements(yrange) eq 0) then yrange=[0,-1] if (n_elements(trange) eq 0) then trange=[0,-1] if (n_elements(rebin) eq 0) then rebin=1 if (n_elements(cutxga) eq 0) then cutxga=0 if (n_elements(sharpen) eq 0) then sharpen=0 if (n_elements(zeromid) eq 0) then zeromid=0 if (n_elements(sqrtint) eq 0) then sqrtint=0 if (n_elements(logint) eq 0) then logint=0 if (n_elements(bytscale) eq 0) then bytscale=2 else bytscale=bytscale if (n_elements(wrappercent) eq 0) then wrappercent=0 if (n_elements(wraprms) eq 0) then wraprms=0 if (n_elements(clock) eq 0) then clock=0 if (n_elements(starttime) eq 0) then starttime=0 if (n_elements(cadence) eq 0) then cadence=0 if (n_elements(datebanner) eq 0) then datebanner=0 if (n_elements(diagname) eq 0) then diagname='' if (n_elements(jpeg) eq 0) then jpeg=0 if (n_elements(gap) eq 0) then gap=0 if (n_elements(overlay) eq 0) then overlay=0 if (n_elements(verbose) eq 0) then verbose=0 ; print pacifier print,' ' print,' --- fitscube2mpeg starts on '+mpegfile ; set endian bigendian=1 ; get unique temporary file identifier rannr=fix(abs(randomn(seed))*10000) strrannr=strmid(string(rannr+1E5,format='(i6)'),2) ; define temporary files for optional sharpening pnginfile='~/tmp/raw'+strrannr+'.png' pngoutfile='~/tmp/sharp'+strrannr+'.png' ; get cube geometry and file datatype from the fits header inheader=headfits_rr(cubefile) if (n_elements(inheader) eq 1) then begin print,' ##### fitscube file not found: ',cubefile wait,5 return endif inheadersize=(1+fix(n_elements(inheader)/36.))*2880 bitpix=fxpar(inheader,'bitpix') nxfile=fxpar(inheader,'naxis1') nyfile=fxpar(inheader,'naxis2') ntfile=fxpar(inheader,'naxis3') ; set dimension ranges as requested, rename to not muck input params xmin=xrange[0] xmax=xrange[1] ymin=yrange[0] ymax=yrange[1] itstart=trange[0] itend=trange[1] if (xmax eq -1) then xmax=nxfile-1 if (ymax eq -1) then ymax=nyfile-1 if (itend eq -1) then itend=ntfile-1 nx=xmax-xmin+1 ny=ymax-ymin+1 nt=itend-itstart+1 blank=bytarr(nx,ny) ; check the requested cutout size if (nx gt nxfile or ny gt nyfile or nt gt ntfile $ or xmax gt nxfile-1 or ymax gt nyfile-1 or itend gt ntfile-1) then begin print,' ### ERROR: xrange or yrange or trange excess' print,' ### nx, ny, nt of fitscube file = ',$ ntostr([nxfile,nyfile,ntfile]) return endif ; get and check size overlay sizeoverlay=size(overlay) if (sizeoverlay[0] gt 0) then begin if (sizeoverlay[0] eq 2 and (sizeoverlay[1] ne nxfile or $ sizeoverlay[2] ne nyfile)) then begin print,' ##### overlay dimensions differ from image' return endif if (sizeoverlay[0] eq 3 and (sizeoverlay[1] ne nxfile or $ sizeoverlay[2] ne nyfile or sizeoverlay[3] ne ntfile)) then begin print,' ##### overlay dimensions differ from movie' return endif endif ; set date time values if (clock ne 0 or datebanner ne 0) then begin headstart=fxpar(inheader,'starttim') ;RR fxpar returns 0 when not present if (headstart ne 0) then starttime=headstart if (headstart eq 0 and starttime eq 0) then message,$ ' ##### oops, clock and/or banner require starttime info' headcad=fxpar(inheader,'cadence') if (headcad ne 0) then cadence=headcad if (headcad eq 0 and cadence eq 0) then message,$ ' ##### oops, clock and/or banner require cadence value' t0time=anytim2tai(starttime) segstart=anytim2utc(t0time+itstart*cadence,/ccsds) date=strmid(segstart,0,10) hh=strmid(segstart,11,2) mm=strmid(segstart,14,2) ss=strmid(segstart,17,2) segstarttai=anytim2tai(segstart) endif ; define time sequence for the clock if (clock eq 1) then begin tt=findgen(nt)*cadence+hh*3600.+mm*60.+ss hr=fix(tt/3600.) mod 24 mn=fix((tt mod 3600)/60) sc=fix(tt mod 60) tt=nnumber(hr,2)+":"+nnumber(mn,2)+":"+nnumber(sc,2) endif ; open input file for assoc get_lun, unit_in if (bigendian) then openr,unit_in,cubefile,/swap_if_little_endian $ else openr,unit_in,cubefile if (bitpix eq -32) then $ inassoc=assoc(unit_in,fltarr(nxfile,nyfile),inheadersize) if (bitpix eq 16) then $ inassoc=assoc(unit_in,intarr(nxfile,nyfile),inheadersize) if (bitpix eq 8) then $ inassoc=assoc(unit_in,bytarr(nxfile,nyfile),inheadersize) ; get min, max of first image minint=min(inassoc[0]) maxint=max(inassoc[0]) insertcolor=255*(avg(inassoc[0])-minint lt (maxint-minint)/2.) ; bytscale > 1: get min and max of whole sequence whole field if (bytscale gt 1) then begin for it=itstart,itend do begin mima=minmax(inassoc[it]) if (mima[0] lt minint) then minint=mima[0] if (mima[1] gt maxint) then maxint=mima[1] endfor if ((bytscale eq 4 and not keyword_set(cutmin)) or $ (bytscale eq 4 and not keyword_set(cutmax))) then begin print,' ##### bytscale = 4 requires cutmin and cutmax' return endif if (verbose eq 1) then print,' ===== minint = '+ntostr(minint)+$ ' maxint = '+ntostr(maxint) endif ; if wrap requested find mean and rms intensities whole field whole sequence if (wrappercent ne 0 or wraprms ne 0) then begin bytscale=0 summean=0. sumrms=0. for it=itstart,itend do begin momim=moment(inassoc[it]) summean=summean+momim[0] sumrms=sumrms+sqrt(momim[1]) endfor avermean=summean/(itend-itstart+1) averrms=sumrms/(itend-itstart+1) if (wrappercent ne 0) then maxint=avermean+wrappercent/100.*avermean if (wraprms ne 0) then maxint=avermean+wraprms*averrms endif if (zeromid ne 0) then begin zeromax=max([abs(minint),abs(maxint)]) if (minint lt 0) then minint=-zeromax if (maxint gt 0) then maxint=zeromax endif if (sqrtint) then begin minint=sqrt(minint) maxint=sqrt(maxint) endif if (logint) then begin minint=alog10(minint) maxint=alog10(maxint) endif ; initalize mpeg ;; mpegid=mpeg_open([nxframe,nyframe], filename=mpegfile,quality=100) ;RR discarded because the IDL mpeg writing commands give bad quality ; start big loop over timesteps it = frame by frame processing ; -------------------------------------------------- for it=itstart,itend do begin ; print progress if (verbose ne 0) then print,$ ' --- timestep ', ntostr(it)+' of ',ntostr(itend) ; get image image=inassoc[it] ; cut as requested image=image[xmin:xmax,ymin:ymax] ; rebin if requested if (rebin ne 1) then image=rebin(image,nx*rebin,ny*rebin) ; cut size to xga if requested if (cutxga ne 0) then begin excessx=max([(nx*rebin-1024)/2,0]) excessy=max([(ny*rebin-768)/2,0]) image=image[excessx:nx*rebin-excessx-1,excessy:ny*rebin-excessy-1] endif ; cut to multiple of 16 px sizeim=size(image) nxframe=(sizeim[1]/16)*16 nyframe=(sizeim[2]/16)*16 image=image[0:nxframe-1,0:nyframe-1] ; check for size excess if (nxframe gt 2048 or nyframe gt 2048) then $ message, ' ### fatal errror: requested mpeg frame exceeds 2048x2048' ; rescale intensity if requested if (sqrtint) then image=sqrt(image) if (logint) then image=alog10(image) ; wrap if requested if (wrappercent ne 0 or wraprms ne 0) then $ image=byte((float(image)-minint)/(maxint-minint)*255) ; bytscale as requested if (bytscale eq 1 or bytscale eq 2) then $ image=bytscl(image,min=minint,max=maxint) if (bytscale eq 3) then image=bytscl(image) if (bytscale eq 4) then begin if (cutmax lt 0) then uppercut=float(-cutmax)/100.*maxint $ else uppercut=cutmax if (cutmin lt 0) then lowercut=float(-cutmin)/100.*abs(minint) $ else lowercut=cutmin if (zeromid) then lowercut=-uppercut image=bytscl(image,min=lowercut,max=uppercut) endif ; sharpen if requested (fixed ImageMagick parameter choice); ; won't work before bytscale when negative values (DO EUV mag, dop); ; doing after bytscale may give black wraps for sharpened brightest, ; but I like that, makes overexposure clearer than blanked cutoff if (sharpen) then begin write_png,pnginfile,image spawn,'convert -unsharp 3.0x1.0+1.5+0 '+pnginfile+' '+pngoutfile image=read_png(pngoutfile) endif ; write the first image as jpeg for usage as presentation clicker if (it eq itstart and jpeg eq 1) then write_jpeg,$ str_replace(mpegfile,'.mpg','_im0.jpg'),image ; insert the clock (use DOT code by Pit Suetterlin) if (clock eq 1) then begin while !d.window ne -1 do wdelete,!d.window window,xsize=nxframe,ysize=nyframe,/pixmap,retain=2 tv, blank clocksize=(nyframe/8)/16*16 ; multiple of 16 px for mpeg clocksize=max([clocksize,32]) clock,tt[it-itstart],/dev,size=clocksize,pos=[nyframe/50,nyframe/15] imclock=tvrd() image[where(imclock gt 1)]=insertcolor endif ; insert running date-time banner if (datebanner ne 0) then begin while !d.window ne -1 do wdelete,!d.window window,xsize=nxframe,ysize=nyframe,/pixmap,retain=2 tv, blank tittai=segstarttai+(it-itstart)*cadence banner=strmid(anytim2utc(tittai,/ccsds),0,19) banner=str_replace(banner,'T',' ') xyouts,nyframe/50,nyframe/50,diagname+' '+banner,color=255,$ /device,charsize=1.*sqrt(nyframe/200.) imbanner=tvrd() image[where(imbanner gt 1)]=insertcolor endif ; insert diagname if no banner if (diagname ne '' and datebanner eq 0) then begin while !d.window ne -1 do wdelete,!d.window window,xsize=nxframe,ysize=nyframe,/pixmap,retain=2 tv,blank xyouts,nyframe/50,nyframe/50,diagname,color=255,$ /device,charsize=1.5 imdiagname=tvrd() image[where(imdiagname gt 1)]=insertcolor endif ; do the overlay insert if requested if (sizeoverlay[0] ne 0) then begin if (sizeoverlay[0] eq 2) then overim=overlay $ else overim=reform(overlay[*,*,it]) ; manipulate overlay to same size as image overim=overim[xmin:xmax,ymin:ymax] if (rebin ne 1) then overim=rebin(overim,nx*rebin,ny*rebin) if (cutxga ne 0) then $ overim=overim[excessx:nx*rebin-excessx-1,excessy:ny*rebin-excessy-1] overim=overim[0:nxframe-1,0:nyframe-1] while !d.window ne -1 do wdelete,!d.window window,xsize=nxframe,ysize=nyframe,/pixmap,retain=2 tv,blank tv,overim overimpx=tvrd() ; apply the overlay (use all non black pixels, so set to 1 for near black) image[where(overimpx ge 1)]=insertcolor endif ;RR need RGB 0-255 > YUV 16-235 format conversion, yak image=byte(fix(image*((235.-16.)/255.)+16.)) ; write movie frame as ~/tmp/png with unique sequence name filename='mpegim'+strrannr+strmid(string(it-itstart+1e6,format='(i7)'),1,6) write_png,'~/tmp/'+filename+'.png',image ;; ; write mpeg movie frame with IDL recipe ;RR NO! too low quality ;; mpeg_put,mpegid,image=image,frame=it,/order endfor ; end of the big loop over images > movie frames ; free the input file free_lun,unit_in ; add gap to mark repeat if requested if (gap ne 0) then begin gapframe=byte(image*0+50) ;RR nice restful dark grey for it=itend+1,itend+gap do begin gapfile='mpegim'+strrannr+strmid(string(it-itstart+1e6,format='(i7)'),1,6) write_png,'~/tmp/'+gapfile+'.png',gapframe endfor endif ; complete the mpeg movie with IDL recipe ; RR NO! too low quality ;; mpeg_save,mpegid ;; mpeg_close,mpegid ; construct the mpeg if (verbose ne 0) then print, ' ---- producing mpeg' ;; spawn,'avconv -f image2 -i ~/tmp/mpegim'+strrannr+$ ;; '%06d.png -v quiet -r 24 -b 65536k -y '+mpegfile ;RR avconv doesn't work in UU science.staff computer so back to ffmpeg spawn,'ffmpeg -f image2 -i ~/tmp/mpegim'+strrannr+$ '%06d.png -r 24 -b 65536k -y -v quiet '+mpegfile ;RR trials to avoid YUV 16-235 format conversion = 'yuv420p', no success ;; '%06d.png -r 24 -pix_fmt rgb24 -b 65536k -y '+mpegfile ;RR autoselecting 'yuv420p'; for .mp2 no output; for .mov clip at YUV values ;; '%06d.png -r 24 -f m24 -b 65536k -y '+mpegfile ; delete all the /tmp images for this movie spawn,'rm -f ~/tmp/mpegim'+strrannr+'* ' spawn,'rm -f '+pnginfile spawn,'rm -f '+pngoutfile ; conclusion if (verbose ne 0) then print,' === wrote movie file ',mpegfile end ; ================= main = test per IDLWAVE H-c ========================== ; this SDO test dir has 304, 171, 1600, 1700, mag, cont, dop fitscube files ; make moving overlay circle for test nx=241 ny=241 nt=11 overlay=bytarr(nx,ny,nt) for it=0,10 do begin window,xsize=nx,ysize=ny,/pixmap,retain=2 tv,bytarr(nx,ny) disksize=(ny/8)/16*16 ; same as clock radius=max([disksize,48])/2 xc=ny/2+it*2 yc=ny/2+it*2 draw_circle,xc,yc,radius,npts=200,/device overlay[*,*,it]=tvrd() endfor cubefile='~/data/SDO/2014-01-10-test/cubes/aia1700.fits' mpegfile='~/tmp/test.mpg' ; [241,241,11] fitscube2mpeg,cubefile,mpegfile,/cutxga,/clock,/datebanner,$ diagname='SDO',/verbose,trange=[5,10],wraprms=5,gap=20,overlay=overlay spawn,'playfs '+mpegfile ;; ;RR blink: should be same, apart from clock + banner and sharpening, wrap ;; incube=readfits(cubefile) ;; movcube=movie2cube(mpegfile) ;; sizemc=size(movcube) ;; inbcube=bytscl(incube[0:sizemc[1]-1,0:sizemc[2]-1,0:sizemc[3]-1]) ;; showex,inbcube,movcube,mag=2 ;; print,'' ;RR needed to get IDLWAVE shell prompt back end