% function qmz = quantizeMS(scans,header,keepScans,scaleFac)% % Quantize mass spec data, and convert it into% a 2D matrix (scans,mz)%% if keepScans is present, then only those scans are % retained  (ie the others are culled)%% scaleFac is from preprocessing, and all abundance values% should be multiplied by scaleFacfunction qmz = quantizeMS(scans,header,keepScans,scaleFac)if (exist('keepScans')==1)    scans=scans(keepScans);    header=header(keepScans);end%convert it to 2D matrix of dimensions: (numScans x 1200)LOW=400; HIGH=1600;% note that numDig represents to what extent we are% rounding off the nearest decimal place (ie how we are% quantizing the mz data).% For example numDigits=2 means that every mass is of % the form XYZ.5  or XYZ.0numDigits=2;%numQuantizedMasses = (HIGH-LOW+1)*numDigits; %% SLIGHTLY WRONGnumQuantizedMasses = mzToCoord(HIGH)-mzToCoord(LOW)+1;data = zeros(length(header),numQuantizedMasses);for s=1:length(scans)%     if (mod(s,100)==0)%         display(strcat('Working on:',num2str(s)));%     end    temp=scans{s};    abundance=temp(:,3);    mz=temp(:,4);    %round mz to nearest 0.5      mzr=round(mz*numDigits)/numDigits;    %check for mz collisions:    collisions=[diff(mzr)];    colInd = find(collisions==0);    if (length(colInd)>0)        %merge collisions        s        error('still getting collisions');        %         currentCol = 1;        %         %are there every more than 2 collisions?        %         bigCollisions = collisions(1:end-1)+collisions(2:end);        %         areBigCol = find(bigCollisions==0);        %         if (length(areBigCol>0))        %             error('must change code to handle larger collisions');        %         end        %        %         newAbund = abundance(1:colInd(1)-1);        %         newmzr = mzr(1:colInd(1)-1);        %         for col=1:length(colInd)        %             tempInd = colInd(col);        %             mergeAbun = sum(abundance(tempInd:(tempInd+1)));        %             newAbund = cat(1,newAbund,mergeAbun);        %             newmzr   = cat(1,newmzr,mzr(tempInd));        %         end        %         data(newmzr-LOW+1,s)=newAbund;    else        %need to map mzr into correct matrix index        %mzInd = (mzr-LOW+1)*numDigits; %% OLD, and *slightly* wrong:        mzInd = mzToCoord(mzr,LOW,numDigits); %% NEW                        if (min(mzr)<LOW) | (max(mzr)>HIGH)            imstats(mzr);            error('out of range, must remove high or low values');        end        data(s,mzInd)=abundance';    endendqmz=data;qmz=qmz*scaleFac;