%include "inc:util.imp"
%include "iff:iffinc.imp"
%include "level1:graphinc.imp"
!Program recombines three monochrome images into an RGB image.

%begin
%record (iffhdr fm) iffin, iffout
%recordformat hf(%halfinteger rgbval, count)
%record (hf) %array h(0:65535)
%string (255) rf, gf, bf, outfile
%halfarray c(0:65535), map(0:255)
%owninteger i,k,rc,r,g,b,ri,gi,bi,x

   %integer entries
   %routine sort(%record (hf) %arrayname h)
      !Sort array h throughout its length in order determined by h_count
      %integer i,j,k,l,n1
      %record (hf) temp

      n1 = 1
      l = entries - 2
      n1 = 2*n1 %while n1<l
      l = 0
       %while n1 > 1 %cycle
         n1 = n1 // 2
         i = l
          %cycle
            i = i+1
            j = i+n1
             %exitif j > entries
            k = i
             %cycle
               %exitif h(j)_count < h(k)_count
               temp = h(j); h(j) = h(k); h(k) = temp
               j = k
               k = k - n1
             %repeatuntil k < 1
          %repeat
       %repeat
    %end                 ; ! of routine SORT

%routine countfile(%integer image, %integerarrayname hist, %integername density)

   %integer i,c
   density=0 ;!32-bit integer should be OK.
   %for i=0, 1, 255 %cycle; hist(i)=0; %repeat
   %for i=0,1,iffin_ht*iffin_wid-1 %cycle
      c = byteinteger(image+i)
      hist(c) = hist(c)+1; density=density+c
   %repeat
   density = int(density/(iffin_ht*iffin_wid))
%end

!***** 1: read in the images *****
printline("Parameters?") %and %stop %unless cliparam -> rf.(",").gf %c
%and gf -> gf.(",").bf
outfile="" %unless bf -> bf.("/").outfile

r=0
rc = iff readin(rf, iffin, r)
iff show header(iffin, 0)
%stop %if rc#0
g=0
rc = iff readin(gf, iffin, g)
iff show header(iffin, 0)
%stop %if rc#0
b=0
rc = iff readin(bf, iffin, b)
iff show header(iffin, 0)
%stop %if rc#0
printline("Images read")

!***** 2: Get the average density of each image *****

!!countfile(r, rc, rd)
!!stretch(r, rc)
!!countfile(g, gc, gd)
!!countfile(b, bc, bd)

!***** 3: Scale and combine the values *****

%for i=0,1,iffin_ht*iffin_wid-1 %cycle
    ri = {int}(byteinteger(r+i)>>3)&31 ;!scales it to range 0-31
    gi = {int}(byteinteger(g+i)>>3)&31
    bi = {int}(byteinteger(b+i)>>3)&31
    c(i) = (bi<<5+gi)<<5+ri ;!15-bit value
%repeat
printline("15-bit Merged image created. ")

!***** 4: Compute a histogram of the scaled values *****

%for i=0,1,65535 %cycle; h(i)=0 ; %repeat

%for i=0, 1, iffin_ht*iffin_wid-1 %cycle
   h(c(i))_count = h(c(i))_count + 1
%repeat
printline("H computed")

!***** 5: Sort the counts into order *****

!Compress the table down first by eliminating "zero" entries
entries=0
%for i=0,1,65535 %cycle
   %if h(i)_count#0 %start
      h(entries)_count=h(i)_count; h(entries)_rgbval = i; h(i)_count=0
      entries=entries+1
   %finish
%repeat
printline("H compacted: ".itos(entries,-1)." non-zero entries")
sort(h)
printline("H sorted")

!***** 6: take the 256 most common colours and put these in the LUT *****

%for i=0,1,255 %cycle; map(i) = h(i)_rgbval; %repeat
printline("Map created")

! also put them back in the COUNT field for us to do the image translation

%for i=0,1,255 %cycle; h(h(i)_rgbval)_count = i; %repeat

! Now translate the image.  odd colours will be left black.

x = heapget(iffin_ht*iffin_wid)
%for i=0,1,iffin_ht*iffin_wid-1 %cycle
   byteinteger(x+i) = h(c(i))_count
%repeat
printline("Image translated")

clear
update colour map(map(0))
col fill(0, 0, iffin_wid-1, iffin_ht-1, byteinteger(x))

!***** 7: Now write the modified picture back to file "rgb.iff" *****

%if outfile # "" %start
iffout=iffin
iffout_mapaddr = addr(map(0)); iffout_maplen=256; iffout_mapwid=16
iffout_hlen=0

rc = iff writeout(outfile, iffout, x)
write(rc, 3); newline
%finish

%endofprogram