!****************************************************************
!*                                                              *
!*      STRETCH:    Program performs contrast enhancement       *
!*                                                              *
!*                  Version 1.1    8 Apr 1988                   *
!*                                                              *
!****************************************************************

%begin
%external%real%fn%spec LOG(%real x)
%include "iff:iffinc.imp"
%include "level1:graphinc.imp"
%include "inc:util.imp"
%include "inc:random.imp"
%bytearrayname a(0:*)
%bytearrayname b(0:*)
%integerarrayname p(0:*)
%integerarrayname q(0:*)
%record (iffhdr fm) iffin, iffout
%real factor,rinc
%integer ad, bd, pd, qd, rc, max, i, j, k, c, pp, qq, nopixels, mode, inc
%constinteger bands=256
%string (255) infile, outfile, param
%halfname cm0
%integerarray tot(0:255)

%routine countfile(%integer from, size, datatype, %integername mode)
   %integer p,q,r
   %recordformat rf(%real r)
   %record (rf) %name rp
   %for p=0,1,255 %cycle; tot(p)=0; %repeat

   %if datatype=0 %start ;!Bytes
      %for q=from, 1, from+size-1 %cycle
         p = byteinteger(q); tot(p) = tot(p) + 1
      %repeat
   %elseif datatype=1 ;!16-bit word
      %for q=from, 2, from+size-2 %cycle
         p = halfinteger(q); tot(p>>8) = tot(p>>8) + 1
      %repeat         
   %elseif datatype=4 ;!32-bit word
      %for q=from, 4, from+size-4 %cycle
         p = integer(q); tot(p>>24) = tot(p>>24) + 1
      %repeat
   %elseif datatype=5 ;!32-bit real
      %for q=from, 4, from+size-4 %cycle
         rp == record(q); p = int(rp_r); tot(p>>24) = tot(p>>24) + 1
      %repeat
   %finish
   q=-1
   %for i=1, 1, 254 %cycle; q=tot(i) %and mode=i %if tot(i)>q; %repeat
%end

%routine stretch(%record (iffhdr fm) %name iffh,
 %bytearrayname a,b(0:*), %real factor)
 %integer i, x
 %for i=0, 1, iffh_ht*iffh_wid-1 %cycle
    x = int(a(i) * factor)
    x=255 %if x>255
    b(i)=x
 %repeat
 col fill(0, 0, iffh_wid-1, iffh_ht-1, b(0))
%end

%routine grey map
   %integer i, c
   %half %array CM (0:255)
   %for i=0,1,255 %cycle; c = i>>3; CM(i) = (c<<5 + c)<<5 + c; %repeat
   Update Colour Map (cm(0))
%end

%routine do count(%bytearrayname a(0:*), %integer margin)
   %integer i, k, mx
   countfile(addr(a(0)), max, iffin_datatype & 7, mode)
   !Grey level histogram.  Bars for values 0 and 255 may be 512
   !pixels long, rest may be 64.
   !Scale the values so we actually see something
!!   %for i=0,1,255 %cycle
!!      %if tot(i)#0 %then tot(i) = int(log(tot(i))*1024)
!!   %repeat
   mx=-1
   mx=tot(0)//8 %if tot(0)//8>mx
   %for i=1,1,254 %cycle; mx=tot(i) %if tot(i)>mx; %repeat
   mx=tot(255)//8 %if tot(255)//8>mx

   %for i=0,1,255 %cycle
      k=margin-(tot(i)*64//mx)
      colour(0)
      fill(margin-64, i*2, k, i*2+1)
      colour(255)
      fill(k,i*2, margin, i*2+1)
   %repeat
   
%end

   %routine sort(%integerarrayname key, data(0:*))
      !Sort array data throughout its length by KEY
      %integer i,j,k,l,n1,temp

      n1 = 1
      l = max - 2
      n1 = n1+n1 %while n1<l
      l = 0
       %while n1 > 1 %cycle
         n1 = n1 >>1
         i = l
          %cycle
            i = i+1
            j = i+n1
             %exitif j > max
            k = i
             %cycle
                %exitif key(j) > key(k)
               temp = key(j); key(j) = key(k); key(k) = temp
               temp = data(j);   data(j) = data(k);   data(k) = temp
               j = k
               k = k - n1
             %repeatuntil k < 1
          %repeat
       %repeat
    %end

%shortfn Random 16
   Seed = (23501*Seed + 6923) & 32767
   %result = Seed
%end

clear

initialise random
param = cli param
outfile="" %and infile=param %unless param -> infile.("/").outfile
ad=0
rc = iff readin(infile, iffin, ad)
%if iffin_maplen=0 %start ;!No colour map - grey levels
   grey map
%else
   cm0 == halfinteger(iffin_mapaddr); update colour map(cm0)
%finish

%if rc=0 %start   ;!If that went OK..
   a == array(ad)
   col fill(0, 0, iffin_wid-1, iffin_ht-1, a(0))
   max = iffin_ht * iffin_wid
   prompt("(M)ultiplier or (H)istogram equalisation:")
   %cycle; readsymbol(c); %repeatuntil c='M' %or c='H'
   bd = heapget(max); b == array(bd)
   do count(a, 688)
   %if c='M' %start

      %cycle
         prompt("Factor:"); read(factor)
         %exit %if factor <=0
         stretch(iffin, a,b,  factor)
         do count(b, 604)
         col fill(0, 0, iffin_wid-1, iffin_ht-1, b(0))
      %repeat
   %else
      pd = heapget(max*4); p == array(pd)
      qd = heapget(max*4); q == array(qd)
      %for i=0,1,max-1 %cycle
         p(i) = byteinteger(ad+i)<<16!random16
         q(i) = i
      %repeat
      printstring("sorting...")
      sort(p, q)
      newline
      nopixels = max//bands
     
      %cycle
         prompt("New modal value (-1 for existing mode, 999 to exit)")
         read(i)
         %exit %if i=999
         mode=i %if 0<=i<=255

         pp=0
         %if mode#0 %start
            rinc = 2*max/(256*mode)
            %for i=0, 1, mode %cycle
               %for j=1, 1, int((i+1)*rinc) %cycle
                  b(q(pp)) = i; pp=pp+1
               %repeat
            %repeat
            newline
         %finish
         qq=max-1
         %if mode#255 %start
            rinc = 2*max/(256*(256-mode))
            %for i=0, 1, 255-mode %cycle
               %for j=1, 1, int((i+1)*rinc) %cycle
                  b(q(qq)) = 255-i; qq=qq-1
               %repeat
            %repeat
         %finish
         %if pp<qq %start
           %for i=pp, 1, qq %cycle
              b(q(i)) = mode
           %repeat
         %finish
         col fill(0, 0, iffin_wid-1, iffin_ht-1, b(0))
         do count(b, 604)
      %repeat
      heapput(pd); heapput(qd)
   %finish
   heapput(ad)

   %if outfile#"" %start
      iffout=iffin; iffin_mapaddr=0
      rc = iff writeout(outfile, iffout, addr(b(0)))
   %finish
%else
   phex(rc); newline
%finish

%endofprogram
